Friday, January 8, 2010

fft and ifft for continuous signal




% for fft look under fft.pdf under engineering docs
clear
clc

% xxxxx
% note:
% also see page 75 computer exploration book
%
% where N samples were taken from x which was non zero for only T secs
% then timeperiod for each sample tau=T/N
%
% it also means that max amount of freq that can detected in this
% sampling is 2pi/tau
%
% which implies w is taken from 2*pi*k/(T=N*tau) where k varies from 0:N
% i.e it will be 0,2pi/N*tau=2pi/T,2pi2/N*tau = 2pi2/T,
% 2piN/T=2piN/N*tau=2pi/tau
%
% where 2pi/tau is sample freq
%
% wk varies as 2pik/(T=N*tau) for  0<=k<=N
%
% for descrete times we do not consider the effect of intermediat time values
% we take tau to be equal to 1
%
% hence wk varies from 0,2pi/N,2pi2/N...2pi
%
% w=2pik/N for 0<=k<=N

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector

% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
y = 0.7*cos(2*pi*50*t) + cos(2*pi*120*t);

% Next power of 2 from length of y
nfft = 2^nextpow2(L);

% Normalized fft => fft/N returns ak see page 47 on exploration book
Y = fft(y,nfft)/L;

% shifted y
Ys = fftshift(fft(y,nfft))/L;

% The FFT contains information between 0 and fs from the fft.pdf in
% engineering docs folder on desktop
f = (Fs)*linspace(0,1,nfft);

% to get only one side of the spectrum
fhalf=(Fs/2)*linspace(0,1,nfft/2+1);

% shifted frequencies
fs_half=(Fs/2)*linspace(-1,1,nfft);

figure(1)
clf

subplot(411)
plot(t(1:50),y(1:50))
title('Sum of sinusoid signals')
xlabel('time (seconds)')

subplot(412)
plot(f,abs(Y))
title('Full Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

% Plot single-sided amplitude spectrum.
subplot(413)
plot(fhalf,abs(Y(1:nfft/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

subplot(414)
plot(fs_half,abs(Ys))
title('Double-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

return

No comments:

Post a Comment