fft_demo.m

数字信号处理（机械工业出版社）的源码

% First create some data.  Consider data sampled at 1000 Hz.  Start by forming a
% time axis for our data, running from t=0 until t=.25 in steps of 1
% millisecond.  Then form a signal, x, containing sine waves at 50 Hz and 120
% Hz.
t = 0:.001:.25;
% creat signal
x = sin(2*pi*50*t) + sin(2*pi*120*t);
% add noise and show it
y = x + 2*randn(size(t));
plot(y(1:50)), title('Noisy time domain signal')
pause
% Finding the discrete Fourier transform of the noisy signal y
Y = fft(y,256);
% Compute the power spectral density
Pyy = Y.*conj(Y)/256;
f = 1000/256*(0:127);
plot(f,Pyy(1:128)), 
title('Power spectral density'), 
xlabel('Frequency (Hz)'),
pause
% Zoom in and plot only up to 200 Hz
plot(f(1:50),Pyy(1:50)), 
title('Power spectral density'), 
xlabel('Frequency (Hz)'),