``````

In [62]:

%pylab inline
from scipy import signal
import math
pylab.rcParams['figure.figsize'] = (16, 6)
pylab.rcParams["font.size"] = "20"

``````
``````

Populating the interactive namespace from numpy and matplotlib

``````

Real values signal FFT, signal, magnitude and phase plots

``````

In [91]:

amp = 0.5      # amplitude of the cosine wave
fc1 = 8.0      # frequency of the cosine wave
fc2 = 48.0
phase = 30.0   # desired phase shift of the cosine in degrees
fs = 100.0     # sampling frequency, 100 samples/s
ts = 1.0 / fs  # sample interval
N = 256        # FFT size
tmax = (2.0*N)/fs # length of signal in seconds
t = arange(0.0, tmax, ts)

``````
``````

In [92]:

# time domain signal with phase shift
sig = amp * cos(2.0 * pi * fc1 * t + phi) + \
amp/4 * cos(2.0 * pi * fc2 * t + phi)
# and plot
grid(True)
plot(t,sig)
title('%3.1f cos(2 pi %d t + pi/6)' % (amp, fc1))
xlabel('time (sec)')
ylabel('x(t)')

``````
``````

Out[92]:

<matplotlib.text.Text at 0x7f7d6c0e5650>

``````
``````

In [93]:

X = fft.rfft(sig) / float(N) # N-point real DFT
frq = fft.rfftfreq(sig.size, d=ts) # get frequency bins

``````
``````

In [94]:

# verify things
print X.size, sig.size, t.size, frq.size, N
print "t  :", t[0], t[-1]
print "frq:", frq[0], frq[-1]
peak_ix = argmax(abs(X))
print "peak of %f at ix %d, %f Hz" % (abs(X[peak_ix]), peak_ix, frq[peak_ix])

``````
``````

257 512 512 257 256
t  : 0.0 5.11
frq: 0.0 50.0
peak of 0.498500 at ix 41, 8.007812 Hz

``````
``````

In [95]:

grid(True)
stem(frq, abs(X)) # magnitudes vs frequencies
xlabel('f (Hz)')
ylabel('|X(k)|')

``````
``````

Out[95]:

<matplotlib.text.Text at 0x7f7d6c209bd0>

``````
``````

In [97]:

mx = max(abs(X))
threshold = mx/10  # tolerance threshold
# maskout values that are below the threshold
X2 = [x if abs(x)>threshold else 0+0j for x in X]
phase = rad2deg(arctan2(imag(X2),real(X2))) # phase information
grid(True)
stem(frq,phase) # phase vs frequencies
xlabel('f (Hz)')
ylabel('dX(k)')

``````
``````

Out[97]:

<matplotlib.text.Text at 0x7f7d6b2b85d0>

``````
``````

In [ ]:

sig_recon = N * fft.irfft(X) # reconstructed signal
grid(True)
plot(t,sig)
title('reconstructed')
xlabel('time (sec)')
ylabel('x(t)')

``````
``````

In [ ]:

``````