In [49]:
%pylab inline
from scipy import signal
import math
pylab.rcParams['figure.figsize'] = (16, 6)
pylab.rcParams["font.size"] = "20"
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]:
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
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]:
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]:
In [54]: