In [49]:
%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

In [57]:
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 [58]:
phi = deg2rad(phase) # phase shift to radians
# 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[58]:
<matplotlib.text.Text at 0x7ff9349c82d0>

In [59]:
windowed = sig * signal.hamming(sig.size, False)
pN = sig.size * 2
p = (pN - sig.size) / 2
padded = np.lib.pad(sig, (p, p), 'constant', constant_values=(0., 0.))
X = fft.rfft(padded) / float(N) # pN-point complex DFT
print sig.size, p, pN, N, X.size


512 256 1024 256 513

In [62]:
frq = fft.rfftfreq(padded.size, d=ts) # get frequency bins
grid(True)
stem(frq, abs(X)) # magnitudes vs frequencies
xlabel('f (Hz)')
ylabel('|X(k)|')
show()
threshold = max(abs(X))/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[62]:
<matplotlib.text.Text at 0x7ff92dcf9650>

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


Out[61]:
<matplotlib.text.Text at 0x7ff92f2e2910>

In [54]: