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

fft and ifft in time and freq domains for descrete signals

% first run prob7_3a b c d.m
clear
clc

nr=[-52:124];

syms n

% see table 5.2 page 392 oppenheim book
% to get X(jw) as 1 for 0<|w|<W and 0 else where
W=pi/2;
x1_n=(W/pi)*sinc(W/pi*n);

% note that x1 is defined for -52<=n<=124
x1=subs(x1_n,n,nr);

% to cover us for nan values
nan=find(isnan(x1));
x1(nan)=1;

N=2048;
wk=[0:N-1]*2*pi/N;

% note here fft was calculated for N=2048 hence it is defines
% N*ak for -52<=k<=-52+N-1
X1=fft(x1,N);

% test of fft and ifft
% since X1 define N*ak for -52<=k<=-52+2048 hence ifft should define x(n)
% for -52<=n<=-52+N-1
x1i=ifft(X1);

% note that X1i should be same as X1
X1i=fft(x1i,N);

figure(1)
clf

subplot(711)
plot(nr,x1);
title('nr against x1');

subplot(712)
plot(wk,abs(X1));
title('abs of X1 over wk for 0 to 2*pi but N*ak for -52<=k<=-52+N-1');

subplot(713)
plot(wk,unwrap(angle(X1)));
title('unwrapped angle for X1 over wk for 0 to 2*pi but N*ak for -52<=k<=-52+N-1');

subplot(714)
plot(wk,abs(X1i));
title('abs of X1i');

subplot(715)
plot(wk,unwrap(angle(X1i)));
title('unwrapped angle for X1i');

subplot(716)
plot([-52:-52+N-1],x1i); %careful here
title('ifft of X1 but over 2048 values starting from -52');

subplot(717)
plot(nr,x1i(1:length(x1))); %careful here
title('ifft of X1 but over only length of x1');

return