[SciPy-User] frequency components of a signal buried in a noisy time domain signal
Nils Wagner
nwagner@iam.uni-stuttgart...
Sat Feb 27 04:06:38 CST 2010
On Fri, 26 Feb 2010 22:03:08 -0500
Ivo Maljevic <ivo.maljevic@gmail.com> wrote:
> Nils,
> I modified a little bit your code. Try this version,
>maybe you will find it
> helpful. I wouldn't bother much
> with the next power of 2 unless you are going to
>implement something in real
> time - FFT calculation
> with non-power of 2 on a PC is just fine. Also, I would
>consider
> periodograms only if the signal changes over time, and
>you want
> to be able to track these changes over time. Otherwise,
>taking the whole
> signal and doing one FFT may be OK.
>
> Take a look at the code below (blue) , which is a
>modified version of your
> code (and read the comments).
>
> Ivo
>
>
> import numpy as np
> import matplotlib.pyplot as plt
> from scipy.fftpack import fft, fftshift
>
> def signal_spectrum(x, Fs):
> n = len(x)
> f = np.linspace(-Fs/2, Fs/2, n)
> X=abs(fftshift(fft(x))) # move the upper spectrum to
>negative freqs
> return [f, X]
>
> pi = np.pi
>Fs = 1000. # Sampling frequency
> T = 1./Fs # Sample time
> L = 1000 # length of signal
> t = np.arange(0,L)*T
> x = 0.7*np.sin(2*pi*50*t)+np.sin(2*pi*120*t)
> y = x + 2*np.random.randn(len(t))
> plt.figure()
> plt.plot(Fs*t[:50],y[:50])
> plt.title('Signal corrupted with zero-mean random
>noise')
> plt.xlabel('Time (milliseconds)')
>
> NFFT = 2**np.ceil(np.log2(L)) # only if you really want
>this
>
> plt.figure()
> Y = fft(y,NFFT)/L
> f = Fs/2*np.linspace(0,1,NFFT/2+1)
> plt.plot(f,2*abs(Y[:NFFT/2+1]))
> plt.title('Single-sided amplitude spectrum of y(t)')
> plt.xlabel('Frequency (Hz)')
> plt.ylabel('|Y(f)|')
>
> # Another way of calculating and plotting the spectrum
> f,Yd = signal_spectrum(y, Fs)
> plt.figure() # you want to add this to be able to see
>multiple plots
> plt.plot(f,Yd)
> plt.title('Spectrum of the the signal')
> plt.xlabel('Frequency (Hz)')
> plt.ylabel('|Y(f)|')
> plt.grid(True)
>
> plt.show()
Hi Ivo,
Thank you very much for your hints.
My experience is, that spectrograms work well with
synthetic signals.
What is needed to use spectrograms for real-life signals ?
I mean do I need a filter etc. ?
Nils
from scipy import linspace, vectorize, cos, pi
from pylab import plot, show, imshow, subplot, specgram,
xlabel, ylabel, savefig
#
# Example taken from
http://en.wikipedia.org/wiki/Short-time_Fourier_transform
#
def x(t):
" Synthetic signal "
if t < 5:
return cos(2*pi*10*t)
if t >= 5. and t < 10:
return cos(2*pi*25*t)
if t >=10. and t< 15:
return cos(2*pi*50*t)
if t >=15. and t<= 20:
return cos(2*pi*100*t)
t = linspace(0.,20.,8001) # sampled at 400 Hz
dt = t[1]-t[0]
NFFT = 512 # the length of the
windowing segments
Fs = int(1.0/dt) # the sampling frequency
x_vec = vectorize(x)
signal = x_vec(t)
ax1=subplot(211)
plot(t,signal)
subplot(212,sharex=ax1)
Pxx, freqs, bins, im = specgram(signal, NFFT=NFFT, Fs=Fs,
noverlap=5)
xlabel(r'Time $t$ [s]')
ylabel('Frequency $f$ [Hz]')
savefig('spectrogram',dpi=60)
show()
More information about the SciPy-User
mailing list