freqdemod-quickstart-3

Author information

John A. Marohn (jam99@cornell.edu)
Department of Chemistry and Chemical Biology
Cornell University
Ithaca, NY USA; 14853-1301

Date

2014/07/01 -- 2014/07/02

Abstract

We use freqdemod's fit() function to recover transient frequency chenges in two frequency-modulated sine-wave signals (1) a signal that we generate here and (2) a signal generated experimentally by Sarah Nathan in a scanning Kelvin probe force microscope experiment.

Preliminaries

Set the working directory to be one directory up from the directory where this file is stored; load the module.



In [1]:

    
import os, sys
module_path = os.path.abspath(os.path.join(os.getcwd(), '..'))
sys.path.append(module_path)

from freqdemod import Signal

Simulated data

Create a frequency-modulated sine wave

Set the digitization frequency to 100 kHz, corresponding to a time step of 10 us. We want a signal with quite a few points, say 512k = 512 x 1024 = 524288 points. The duration of the signal will thus be 524288 x 10E-6 = 5.24288 seconds.

The sine wave has an initial frequency of 4 kHz. Between 1.5 s and 3.5 s, the frequency is ramped linearly to 6 kHz. It remains 6 kHz until the end of the signal.



In [2]:

    
import math
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# set values

fd = 100E3
nt = 512*1024

f_start = 4.000E3
f_end = 6.000E3

t_start = 1.5
t_end = 3.5

# time array

dt = 1/fd

t1 = dt*np.arange(int(round(t_start/dt)))
t2 = dt*np.arange(int(round((t_end-t_start)/dt)))
t3 = dt*np.arange(nt-t1.size-t2.size)

t2plus = t2+t1[-1]+dt
t3plus = t3+t2plus[-1]+dt

t = np.append(np.append(t1,t2plus),t3plus)

# frequency array

f1 = f_start*np.ones(t1.size)
f2 = f_start + t2*(f_end-f_start)/(t2[-1]-t2[0])
f3 = f_end*np.ones(t3.size)

f = np.append(np.append(f1,f2),f3)

# phase accumulator

p = np.zeros(t.size)

p[0] = 0.0
for k in np.arange(1,t.size):
    p[k] = p[k-1] + dt*f[k-1]

p = 2*np.pi*p

# the signal

s = np.cos(p)

Create the Signal objec, demodulate it, and fit the phase vs time data

... to recover the oscillation frequency. We have a few choices to make:

We choose window rise/fall time of 1 ms
We choose a filter bandwidth of 5 kHz, larger than the 2 kHz width of the frequency sweep.
Finally we need to choose the "chunk" time. We want to choose this so that the associated Nyquist frequency is larger than the filter bandwidth. This is so that we do not alias any noise into the low-frequency part of the frequency power spectrum. Choose 1/(10 kHz) = 100 us



In [3]:

    
S = Signal(s,"x","nm",1/fd)
S.binarate("middle")
S.window(tw=1.0E-3)
S.fft()
S.filter(bw=5.0E3)
S.ifft()
S.trim()
S.fit(dt_chunk_target=100E-6)
print(S)









    



Signal
======
signal name: x
signal unit: nm
signal lenth = 524288
time step = 10.000 us
rms = 0.7071
max = 1
min = -1
 
Signal Report
=============
* Add a signal x[nm] of length 524288, time step 10.000 us, and duration 5.243 s

* Truncate the signal to be 524288 points long, a power of two. This was done by chopping points off the beginning and end of the signal array, that is, by using points 0 up to 524288.

* Window the signal with a rising/falling blackman filter having a rise/fall time of 1000.000 us (100 points).

* Fourier transform the windowed signal. It took 797.2 ms to compute the FFT.

* Reject negative frequencies, apply a bandpass filter of bandwidth 5000.0 Hz & order 50, and set the delay time to 250.0 us.

* Apply an inverse Fourier transform.

* Remove the leading and trailing ripple (25 points) from the complex signal. Compute the signal phase and amplitude.

* Curve fit the phase data. The target chunk duration is 100.000 us; the actual chunk duration is 100.000 us (10 points). 52423 chunks will be curve fit; 5242.300 ms of data. It took 29.2 ms to perform the curve fit and obtain the frequency.

Plot the resulting frequency vs time



In [4]:

    
S.plot_phase_fit()

Experimental data

Read in the signal

... which is stored, compactly, in a HDF5 type file.



In [5]:

    
import h5py
with h5py.File('JAM_9p90__copy.h5', 'r') as h5f:
    s = h5f['s'][:]

Create the Signal objec, demodulate it, and fit the phase vs time data

... to recover the oscillation frequency. We have a few choices to make:

We choose window rise/fall time of 10 ms.
We choose a filter bandwidth of 1 kHz -- large enough to capture a ms transient but small enough to give us some signal averaging.
Choose the "chunk" time to be 0.5 ms; the associated Nyquist frequency is 2 kHz -- larger than the filter bandwidth.



In [6]:

    
S = Signal(s,"x","nm",1/1.15E6)
S.binarate("middle")
S.window(tw=10.0E-3)
S.fft()
S.filter(bw=1.0E3)
S.ifft()
S.trim()
S.fit(dt_chunk_target=0.50E-3)
print(S)









    



Signal
======
signal name: x
signal unit: nm
signal lenth = 524288
time step = 0.870 us
rms = 0.1486
max = 0.1885
min = -0.2339
 
Signal Report
=============
* Add a signal x[nm] of length 762375, time step 0.870 us, and duration 0.663 s

* Truncate the signal to be 524288 points long, a power of two. This was done by chopping points off the beginning and end of the signal array, that is, by using points 119043 up to 643331.

* Window the signal with a rising/falling blackman filter having a rise/fall time of 10000.000 us (11500 points).

* Fourier transform the windowed signal. It took 153.0 ms to compute the FFT.

* Reject negative frequencies, apply a bandpass filter of bandwidth 1000.0 Hz & order 50, and set the delay time to 1250.0 us.

* Apply an inverse Fourier transform.

* Remove the leading and trailing ripple (1437 points) from the complex signal. Compute the signal phase and amplitude.

* Curve fit the phase data. The target chunk duration is 500.000 us; the actual chunk duration is 500.000 us (575 points). 906 chunks will be curve fit; 453.000 ms of data. It took 30.4 ms to perform the curve fit and obtain the frequency.

Plot the frequency vs time



In [7]:

    
print "number of elements = {}".format(S.signal['fit_freq'].size)
print "number of nan = {}".format(np.count_nonzero(np.isnan(S.signal['fit_freq'])))

S.plot_phase_fit()









    



number of elements = 906
number of nan = 0

Replot the frequency shift

There is an interesting blip in the middle that we want to take a closer look at. To do this, compute a frequency shift; zoom in to the middle section of the data; and plot vs time shift in ms.



In [9]:

    
import matplotlib.pyplot as plt
import copy

# make a copy of the fitten data

fit_freq_temp = copy.deepcopy(S.signal['fit_freq'])
fit_time_temp = copy.deepcopy(S.signal['fit_time'])

# caculate a frequency shift by subtrating off
#  the average over the first third of the data

fit_freq_size = fit_freq_temp.size
fit_freq_size_avg = fit_freq_temp[0:fit_freq_size/3].mean()
fit_freq_temp = fit_freq_temp - fit_freq_size_avg

# define a zoom window -- define as a perecentage
#  the the window still works if we change the chunk time
#  for curve fitting

start = 0.47*fit_freq_size
stop = 0.52*fit_freq_size

# plot the frequency shift [Hz] versus time [ms]
#  reset the time at the start of the plot to zero

x = 1E3*(fit_time_temp[start:stop] - fit_time_temp[start])
y = fit_freq_temp[start:stop]

plt.rcParams['text.usetex'] = True

plt.plot(x,y,'ko-')
plt.ylabel(r"$\Delta f \: \mathrm{[Hz]}$")
plt.xlabel(r"$\Delta t \: \mathrm{[ms]}$")
plt.show()