[SciPy-User] frequency components of a signal buried in a noisy time domain signal

[Nils Wagner]

On Fri, 26 Feb 2010 22:03:08 -0500
  Ivo Maljevic <ivo.maljevic at 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()