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)

Bases: object

__init__(filename=None)

Initialize the Signal object. Inputs:

Parameters:filename (str) – the signal’s (future) filename

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 (e.g., “Empty signal object created”)

If no filename is given, then create an object containing just a report. Note that you can’t do anything with this empty object except to call the open function. If you intend to use the object but not save it to disk, then you should create it with a dummy filename.

load_nparray(s, s_name, s_unit, dt)

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]

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

close()

Update report; write the file to disk; close the file.

open(filename)

Open the file for reading and writing. The report comes back as a np.ndarray. Need to convert back to a 1D array by flattening, then convert to a list so we can continue appending.

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)

Create a bandpass filter with

Parameters:

• bw (float) – filter bandwidth. \(\Delta f\) [kHz]
• order (int) – filter order, \(n\) (defaults to 50)

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

workup/freq/filter/bp  

The associated filtering function 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\).

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.print_hdf5_item_structure(g, offset=' ')

Prints the input file/group/dataset (g) name and begin iterations on its content

demodulate.testsignal_sine()

demodulate.testsignal_sine_fm()

demodulate.testsignal_sine_exp()