freqdemod-develop-1

Author information

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

Date

2014/06/28 created
2014/11/23 revised

Abstract

Here we apply freqdemod to a noiseless sine-wave oscillation, obtaining phase vs time data from oscillation amplitude vs time data. We fit the phase vs time data to a line and obtain a slope that agrees exactly with the known frequency of the sine wave.

Preliminaries

We first have to tell Python where to find the freqdemod module. Since we are running this file from the /docs directory of the module diretory, we have to tell sys.path to look one level up to find the module.



In [1]:

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

Now we are ready to load the module:



In [2]:

    
from freqdemod import Signal









    



---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-2-a1d53745a648> in <module>()
----> 1 from freqdemod import Signal

ImportError: cannot import name Signal

Create and plot a sine wave



In [3]:

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

f0 = 7.457E3

fd = 50.0E3
dt = 1/fd
nt = 300

t = dt*np.arange(nt)
s = np.sin(2*np.pi*f0*t)
s_name = "x"
s_unit = "nm"

plt.plot(t*1E6,s)
plt.xlabel("t [us]")
plt.ylabel(s_name + " [" + s_unit + "]")
plt.show()

Create the Signal object



In [4]:

    
S = Signal(s,s_name,s_unit,dt)
print(S)









    



Signal
======
signal name: x
signal unit: nm
signal lenth = 300
time step = 20.000 us
rms = 0.7059
max = 1
min = -1
 
Signal Report
=============
* Add a signal x[nm] of length 300 and time step 20.000 us.

Trim and window the signal

Before proceeding with the Fast Fourier Transform, force the data to be a factor of two in length



In [5]:

    
S.binarate("middle")
print(S)









    



Signal
======
signal name: x
signal unit: nm
signal lenth = 256
time step = 20.000 us
rms = 0.7084
max = 1
min = -1
 
Signal Report
=============
* Add a signal x[nm] of length 300 and time step 20.000 us.

* Truncate the signal to be 256 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 22 up to 278.

Window the data so that is starts and ends at zero. Look at two possible windowing functions. Use the second one. Below we plot, versus time in microseconds, the unwindowed signal (green), the windowed signal (red), and the windowing function (blue).



In [6]:

    
S.window(128*dt)
plt.plot(S.signal['t']*1E6,S.signal['w'])
plt.plot(S.signal['t']*1E6,S.signal['s'])
plt.plot(S.signal['t']*1E6,S.signal['sw'])
plt.show()



In [7]:

    
S.window(32*dt)
plt.plot(S.signal['t']*1E6,S.signal['w'])
plt.plot(S.signal['t']*1E6,S.signal['s'])
plt.plot(S.signal['t']*1E6,S.signal['sw'])
plt.show()

FT, filter, IFFT

Fourier transform (FT) the data. Plot the real and imaginary parts of the FT'ed data versus frequency (upper plot; red and green lines). Separately, plot the magnitude (lower plot; blue line).



In [8]:

    
S.fft()
plt.plot(S.signal['f'],S.signal['swFT'].real)
plt.plot(S.signal['f'],S.signal['swFT'].imag)
plt.show()
plt.plot(S.signal['f'],abs(S.signal['swFT']))
plt.show()

Apply the right-hand and bandpass filters. In the plot below we see plotted versus frequency on a semi-log scale the FT'ed signal (green), the right handed filter function (red), the bandpass filter function (purple), and the filtered FT'ed signal (blue). The filter order is the default -- 50. Note how sharp the edges of the filter window are; this is a "brick wall" filter.



In [9]:

    
S.filter(bw=2000.0)

yf = abs(S.signal['swFTfilt'])
y = abs(S.signal['swFT'])

plt.plot(S.signal['f'],yf/max(yf))
plt.plot(S.signal['f'],y/max(y))
plt.plot(S.signal['f'],S.signal['rh'])
plt.plot(S.signal['f'],S.signal['bp'])
plt.plot(S.signal['f'],S.signal['rh']*S.signal['bp'])
plt.yscale('log')
plt.ylim([1E-7,3.0])
plt.show()

IFFT the data to obtain a plot of phase in-phase and out-of-phase amplitude versis time (see the greed and blue curves in the below).



In [10]:

    
S.ifft()



In [11]:

    
t = 1E6*S.signal['t'][0:50]
y1 = S.signal['z'].real[0:50]
y2 = S.signal['z'].imag[0:50]

plt.plot(t,y1)
plt.plot(t,y2)
plt.xlabel("t [us]")
plt.show()

Trim the complex signal, phase, and amplitude data to remove the initial and final ripple.



In [12]:

    
S.trim()



In [13]:

    
plt.plot(S.signal['t']*1E6,S.signal['theta'])
plt.xlim(0,S.signal['t_original'][-1]*1E6)
plt.ylabel("phase [cycles/sec]")
plt.xlabel("t [us]")
plt.show()

plt.plot(S.signal['t']*1E6,S.signal['a'])
plt.xlim(0,S.signal['t_original'][-1]*1E6)
plt.ylabel("amplitude")
plt.xlabel("t [us]")
plt.show()

Fit the phase vs time data

Fit the phase vs time data to a line. The slope of the line is the frequency. We fit S.signal['theta'] (y-axis) vs S.signal['t'] (x-axis). Use numpy's polyfit function to carry out the fitting.



In [14]:

    
zf = np.polyfit(S.signal['t'], S.signal['theta'], 1)

print "Best-fit frequency = {0:.3f} Hz".format(zf[0])
print "   Known frequency = {0:.3f} Hz".format(f0)









    



Best-fit frequency = 7457.001 Hz
   Known frequency = 7457.000 Hz

The frequency inferred from the slope of the phsse vs frequency line agrees with the known frequency to six decminal places.

Summary

Let's summarize the procedure



In [15]:

    
# Load the neccessary packages

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

# Generate a signal
#
# f0 = signal frequency
# fd = digitization frequency
# nt = number of signal points

f0 = 4.457E3  
fd = 50.0E3 
nt = 512

dt = 1/fd
t = dt*np.arange(nt)
s = np.sin(2*np.pi*f0*t)

R = Signal(s,"x","nm",1/fd)

# Analyze the signal

R.binarate("middle")
R.window(201E-6)
R.fft()
R.filter(bw=4E3)
R.ifft()
R.trim()

# Print the signal

print(R)









    



Signal
======
signal name: x
signal unit: nm
signal lenth = 512
time step = 20.000 us
rms = 0.7062
max = 1
min = -1
 
Signal Report
=============
* Add a signal x[nm] of length 512 and time step 20.000 us.

* Truncate the signal to be 512 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 512.

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

* Fourier transform the windowed signal.

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

* Apply an inverse Fourier transform.

* Remove the leading and trailing ripple from the complex signal. Compute the signal phase and amplitude.



In [16]:

    
# Plot the signal

plt.plot(1E6*R.signal['t'][0:100],R.signal['z'].real[0:100])
plt.plot(1E6*R.signal['t'][0:100],R.signal['z'].imag[0:100])
# plt.xlim(0,R.signal['t_original'][-1]*1E6)
plt.ylabel(s_name + " [" + s_unit + "]")
plt.xlabel("t [us]")
plt.show()

plt.plot(R.signal['t']*1E6,R.signal['theta'])
plt.xlim(0,R.signal['t_original'][-1]*1E6)
plt.ylabel("phase [cycles]")
plt.xlabel("t [us]")
plt.show()

plt.plot(R.signal['t']*1E6,R.signal['a'])
plt.xlim(0,R.signal['t_original'][-1]*1E6)
plt.ylabel("amplitude")
plt.xlabel("t [us]")
plt.show()