demodulate module

We demodulate the signal in the following steps:

1. Apply a window to the signal \(S(t)\) in order to make it smoothly ramp up from zero and ramp down to zero.
2. Fast Fourier Transform the windowed signal to obtain \(\hat{S}(f)\).
3. Identify the primary oscillation frequency \(f_0\). Apply a bandpass filter centered at \(f_0\) to reject signals at other frequencies.
4. Apply a second filter which zeros out the negative-frequency components of \(\hat{S}(f)\).
5. Apply an Inverse Fast Fourier Transform to the resulting data to obtain a complex signal \(z(t) = x(t) + \imath \: y(t)\).
6. Compute the instantaneous phase \(\phi\) and amplitude \(a(t)\) using the following equations. Unwrap the phase.

(1)\[\begin{equation} \phi(t) = \arctan{[\frac{y(t)}{x(t)}]} \end{equation}\]

(2)\[\begin{equation} a(t) = \sqrt{x(t)^2 + y(t)^2} \end{equation}\]

7. Calculate the “instantaneous” frequency \(f(t)\) by dividing the instantaneous phase data into equal-time segments and fitting each segment to a line. The average frequency \(f(t)\) during each time segment is the slope of the respective line.

class demodulate.Signal(filename=None, mode='w-', driver='core', backing_store=False)

Bases: object

__init__(filename=None, mode='w-', driver='core', backing_store=False)

Initialize the Signal’s hdf5 data structure. Calling with no arguments results in an in-memory only object. To save the data to disk, provide both a filename and backing_store=True.

Parameters:

• filename (str) – the signal’s filename.
• mode (str) – file open mode (see h5py.File)
• driver (str) – hdf5 driver (see h5py.File)
• backing_store (bool) – If True, save the file to disk.

Add the following objects to the Signal object

Parameters:

• f – an h5py object stored in core memory
• report (str) – a string summarizing in words what has been done to the signal (e.g., “Empty signal object created”)

Examples:

s = Signal()                                      # In-memory only
s = Signal('not-saved.h5')                        # Still in-memory only
s = Signal('save-to-disk.h5', backing_store=True) # Save to disk

load_nparray(s, s_name, s_unit, dt, s_help='cantilever displacement')

Create a Signal object from the following inputs.

Parameters:

• s (np.array or list) – the signal vs time
• s_name (str) – the signal’s name
• s_name – the signal’s units
• dt (float) – the time per point [s]
• s_help (str) – the signal’s help string

Add the following objects to the Signal object

Parameters:

• f (h5py.File) – an h5py object stored in core memory
• report (str) – a string summarizing in words what has been done to the signal

load_hdf5(f, s_dataset='y', t_dataset='x', infer_dt=True, infer_attrs=True)

Load an hdf5 file saved with default freqdemod attributes.

Parameters:

• f – A filename, or h5py file or group object, which contains t_dataset and s_dataset
• s_dataset (str) – signal dataset name (relative to h5object)
• t_dataset (str) – time dataset name (relative to h5object)
• infer_dt (bool) – If True, infer the time step dt from the contents of the t_dataset
• infer_attrs (bool) – If True, fill in any missing attributes used by freqdemod.

load_hdf5_general(f, s_dataset, s_name, s_unit, t_dataset=None, dt=None, s_help='cantilever displacement')

Load data from an arbitrarily formatted hdf5 file.

Parameters:

• f – A filename, or h5py file or group object, which contains s_dataset
• s_dataset (str) – signal dataset name (relative to h5object)
• s_name (str) – the signal’s name
• s_unit (str) – the signal’s units
• t_dataset (str) – time dataset name (optional; or specify dt)
• dt (float) – the time per point [s]
• s_help (str) – the signal’s help string

close()

Update report, close the file. This will write the file to disk if backing_store=True was used to create the file.

save(dest, save='time_workup', overwrite=False)

Save the current signal object to a new hdf5 file, with control over what datasets to save.

Parameters:

• dest – Copy destination. A filename, or h5py file or group object.
• save – A string describing the datasets to be saved, or a list of groups / datasets to save. all: All datasets and groups input: x, y input_no_t: y time_workup: x, y, workup/time time_workup_no_t: y, workup/time time_workup_no_s: workup/time fit_phase: x, y, workup/fit fit_phase_no_t: y, workup/fit fit_phase_no_s: workup/fit
• overwrite – If true, overwrite an existing destination file.

open(filename)

Open the file for reading and writing.

plot(ordinate, LaTeX=False, component='abs')

Plot a component of the Signal object.

param str ordinate:  the name the y-axis data key param str component:  abs (default), real, imag, or both; if the dataset is complex, which component do we plot param boolean LaTeX:  use LaTeX axis labels; True or False (default)

Plot self.f[ordinate] versus self.f[y.attrs[‘abscissa’]]

If the data is complex, plot the absolute value.

time_mask_binarate(mode)

Create a masking array of True/False values that can be used to mask an array so that an array is a power of two in length, as required to perform a Fast Fourier Transform.

Parameters:mode (str) – “start”, “middle”, or “end”

With “start”, the beginning of the array will be left intact and the end truncated; with “middle”, the array will be shortened symmetically from both ends; and with “end” the end of the array will be left intact while the beginning of the array will be chopped away.

Create the mask in the following np.array:

Parameters:self.f['workup/time/mask/binarate'] (bool) – the mask; a np.array of boolean values

This boolean array of True and False values can be plotted directly – True is converted to 1.0 and False is converted to 0 by the matplotlib plotting function.

time_window_cyclicize(tw)

Create a windowing function with

Parameters:tw (float) – the window’s target rise/fall time [s]

The windowing function is a concatenation of

1. the rising half of a Blackman filter;
2. a constant 1.0; and
3. the falling half of a Blackman filter.

The actual rise/fall time may not exactly equal the target rise/fall time if the requested time is not an integer multiple of the signal’s time per point.

If the masking array workup/time/mask/binarate is defined then design the windowing function workup/time/window/cyclicize so that is is the right length to be applied to the masked signal array. If on the other hand the masking array is not already defined, then densign the window to be applied to the full signal array instead.

If workup/time/mask/binarate is defined, then also create a masked version of the x-axis for plotting, x_binarated. Make the x_binarated array the new abscissa associated with the new workup/time/window/cyclicize array.

fft()

Take a Fast Fourier transform of the windowed signal. If the signal has units of nm, then the FT will have units of nm/Hz.

freq_filter_Hilbert_complex()

Generate the complex Hilbert transform filter (\(Hc\) in the attached notes). Store the filter in:

workup/freq/filter/Hc

The associated filtering function is

\[\begin{split}\begin{equation} \mathrm{rh}(f) = \begin{cases} 0 & \text{if $f < 0$} \\ 1 & \text{if $f = 0$} \\ 2 & \text{if $f > 0$} \end{cases} \end{equation} \end{split}\]

freq_filter_bp(bw, order=50, style='brick wall')

Create a bandpass filter with

Parameters:

• bw (float) – filter bandwidth, \(\Delta f\) [kHz]
• order (int) – filter order, \(n\) (defaults to 50)
• style (string) – “brick wall” (default) or “cosine”

Note that the filter width bw should be specified in kHz and not Hz. Store the filter in:

workup/freq/filter/bp  

The center frequency, \(f_0\), is determined automatically from the positive-frequency peak in the Fourier-transform of the cantilever signal. The filtering function for the brick wall filter is

\[\begin{equation} \mathrm{bp}(f) = \frac{1}{1 + (\frac{|f - f_0|}{\Delta f})^n} \end{equation}\]

The absolute value makes the filter symmetric, even for odd values of \(n\). With the bw parameter set to 1 kHz, the filter will pass a 2 kHz band of frequencies, from 1 kHz below to 1 kHz above the center frequency. The power spectrum will show noise falling away starting 1 kHz away from the carrier.

The filtering function for the cosine filter is

\[\begin{split}\begin{equation} \mathrm{bp}(f) = \begin{cases} 0 & f < f_0 - \Delta f \\ \cos{\left( \frac{\pi}{2} \frac{f-f_0}{\Delta f} \right)} & f_0 - \Delta f < f < f_0 + \Delta f \\ 0 & f_0 + \Delta f < f \end{cases} \end{equation}\end{split}\]

time_mask_rippleless(td)

Defined using

Parameters:td (float) – the dead time [s]

that will remove the leading and trailing rippled from the phase (and amplitude) versus time data.

Programming notes. As a result of the applying cyclicizing window in the time domain and the bandpass filter in the frequency domain, the phase (and amplitude) will have a ripple at its leading and trailing edge. We want to define a mask that will let us focus on the uncorrupted data in the middle of the phase (and amplitude) array

Defining the required mask takes a little thought.We need to first think about what the relevant time xis is. If:

self.f.__contains__('workup/time/mask/binarate') == True:

then the data has been trimmed to a factor of two and the relevant time axis is:

self.f['/workup/time/x_binarated']

Otherwise, the relevant time axis is:

self.f['x'] 

We will call this trimming mask:

self.f['workup/time/mask/rippleless']

We will apply this mask to either self.f['/workup/time/x_binarated'] or self.f['x'] to yield a new time array for plotting phase (and amplitude) data:

self.f['/workup/time/x_rippleless']

If no bandpass filter was defined, then the relevant time axis for the phase (and amplitude) data is either:

self.f['/workup/time/x_binarated']
(if self.f.__contains__('workup/time/mask/binarate') == True)    

or:

self.f['x']
(if self.f.__contains__('workup/time/mask/binarate') == False) 

ifft()

If they are defined,

• apply the complex Hilbert transform filter,
• apply the bandpass filter,

then

• compute the inverse Fourier transform,
• if a trimming window is defined then trim the result

fit_phase(dt_chunk_target)

Fit the phase vs time data to a line. The slope of the line is the (instantaneous) frequency. The phase data is broken into “chunks”, with

Parameters:dt_chunk_target (float) – the target chunk duration [s]

If the chosen duration is not an integer multiple of the digitization time, then find the nearest chunk duration which is.

Calculate the slope \(m\) of the phase vs time line using the linear-least squares formula

\[\begin{equation} m = \frac{n \: S_{xy} - S_x S_y}{n \: S_{xx} - (S_x)^2} \end{equation}\]

with \(x\) representing time, \(y\) representing phase, and \(n\) the number of data points contributing to the fit. The sums involving the \(x\) (e.g., time) data can be computed analytically because the time data here are equally spaced. With the time per point \(\Delta t\),

\[\begin{equation} S_x = \sum_{k = 0}^{n-1} x_k = \sum_{k = 0}^{n-1} k \: \Delta t = \frac{1}{2} \Delta t \: n (n-1) \end{equation}\]\[\begin{equation} S_{xx} = \sum_{k = 0}^{n-1} x_k^2 = \sum_{k = 0}^{n-1} k^2 \: {\Delta t}^2 = \frac{1}{6} \Delta t \: n (n-1) (2n -1) \end{equation}\]

The sums involving \(y\) (e.g., phase) can not be similarly precomputed. These sums are

\[\begin{equation} S_y = \sum_{k = 0}^{n-1} y_k = \sum_{k = 0}^{n-1} \phi_k \end{equation} \]\[\begin{equation} S_{xy} = \sum_{k = 0}^{n-1} x_k y_k = \sum_{k = 0}^{n-1} (k \Delta t) \: \phi_k \end{equation} \]

To avoid problems with round-off error, a constant is subtracted from the time and phase arrays in each chuck so that the time array and phase array passed to the least-square formula each start at zero.

fit_amplitude()

Fit the data stored at:

workup/time/a

To a decaying exponential. Store the resulting fit at:

workup/fit/exp

Programming Notes

• From the lmfit documentation [link]: “Note that the calculation of chi-square and reduced chi-square assume that the returned residual function is scaled properly to the uncertainties in the data. For these statistics to be meaningful, the person writing the function to be minimized must scale them properly.”

plot_fit(fit_group, LaTeX=False)

Plot the fit stored in fit_group.

__repr__()

Print out the report.

list(offset='', indent=' ')

List all file/group/dataset elements in the hdf5 file by iterating over the file contents.

Source:

https://confluence.slac.stanford.edu/display/PSDM
/How+to+access+HDF5+data+from+Python
#HowtoaccessHDF5datafromPython-HDF5filestructure

demodulate.testsignal_sine()

demodulate.testsignal_sine_fm()

demodulate.testsignal_sine_exp()