demodulate
module¶
We demodulate the signal in the following steps:
 Apply a window to the signal \(S(t)\) in order to make it smoothly ramp up from zero and ramp down to zero.
 Fast Fourier Transform the windowed signal to obtain \(\hat{S}(f)\).
 Identify the primary oscillation frequency \(f_0\). Apply a bandpass filter centered at \(f_0\) to reject signals at other frequencies.
 Apply a second filter which zeros out the negativefrequency components of \(\hat{S}(f)\).
 Apply an Inverse Fast Fourier Transform to the resulting data to obtain a complex signal \(z(t) = x(t) + \imath \: y(t)\).
 Compute the instantaneous phase \(\phi\) and amplitude \(a(t)\) using the following equations. Unwrap the phase.
 Calculate the “instantaneous” frequency \(f(t)\) by dividing the instantaneous phase data into equaltime 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)[source]¶ Bases:
object

__init__
(filename=None, mode='w', driver='core', backing_store=False)[source]¶ Initialize the Signal’s hdf5 data structure. Calling with no arguments results in an inmemory 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() # Inmemory only s = Signal('notsaved.h5') # Still inmemory only s = Signal('savetodisk.h5', backing_store=True) # Save to disk

load_nparray
(s, s_name, s_unit, dt, s_help='cantilever displacement')[source]¶ 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)[source]¶ 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')[source]¶ 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
()[source]¶ 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)[source]¶ 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.

plot
(ordinate, LaTeX=False, component='abs')[source]¶ Plot a component of the Signal object.
param str ordinate: the name the yaxis 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
orFalse
(default)Plot
self.f[ordinate]
versus self.f[y.attrs[‘abscissa’]]If the data is complex, plot the absolute value.

time_mask_binarate
(mode)[source]¶ 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 valuesThis boolean array of
True
andFalse
values can be plotted directly –True
is converted to 1.0 andFalse
is converted to 0 by thematplotlib
plotting function.

time_window_cyclicize
(tw)[source]¶ Create a windowing function with
Parameters: tw (float) – the window’s target rise/fall time [s] The windowing function is a concatenation of
 the rising half of a Blackman filter;
 a constant 1.0; and
 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 xaxis for plotting, x_binarated. Make the x_binarated array the new abscissa associated with the new workup/time/window/cyclicize array.

fft
()[source]¶ 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
()[source]¶ 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')[source]¶ 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 positivefrequency peak in the Fouriertransform 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{ff_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)[source]¶ 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']
orself.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
()[source]¶ 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)[source]¶ 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 linearleast 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}^{n1} x_k = \sum_{k = 0}^{n1} k \: \Delta t = \frac{1}{2} \Delta t \: n (n1) \end{equation}\]\[\begin{equation} S_{xx} = \sum_{k = 0}^{n1} x_k^2 = \sum_{k = 0}^{n1} k^2 \: {\Delta t}^2 = \frac{1}{6} \Delta t \: n (n1) (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}^{n1} y_k = \sum_{k = 0}^{n1} \phi_k \end{equation} \]\[\begin{equation} S_{xy} = \sum_{k = 0}^{n1} x_k y_k = \sum_{k = 0}^{n1} (k \Delta t) \: \phi_k \end{equation} \]To avoid problems with roundoff 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 leastsquare formula each start at zero.

fit_amplitude
()[source]¶ 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 chisquare and reduced chisquare 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.”
 From the
