%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Code for DSP Lecture of December 29, 2007
%
%  Different parts are separated by "pause"
%   You need to press a key to proceed
%   to following part
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all

%%%%%%%%%%% Basic signal and its DFT
N= 256;
sig= zeros(N,1);
sig(N/4:N/2-1-40)= 1.;
figure(1);
subplot(4,1,1); plot(sig)
fsig0= fft(sig);
fsig= fftshift(fft((sig)));
subplot(4,1,2); plot(abs(fsig0));
subplot(4,1,3); plot(abs(fsig))
subplot(4,1,4); plot((phase(fsig)))
pause


%%%%%%%% Delay of the signal by m samples
N= 2048;
sig= zeros(N,1);
sig(1:N/4)= 1.;
fsig= fftshift(fft(fftshift(sig)));
m=N/13.;                    %%%% ripples appear when noninteger delay is used
ind= -N/2:N/2-1;
fsig1= fsig.*exp(-j*2*pi*ind'/N*m);
sig1= fftshift(ifft(fftshift(fsig1)));
figure (4)
subplot(3,1,1); plot((sig))
subplot(3,1,2); plot(abs(sig1))
% example of a circular shift
m= N/16.;
fsig1= (fsig).*exp(-j*2*pi*ind'/N*m);
sig1= fftshift(ifft(fftshift(fsig1)));
subplot(3,1,3); plot(abs(sig1))
pause


%%%%%%% basic linear/circular convolution properties
% case when both are similar
N1= 16;
s1= zeros(N1,1);
s2= zeros(N1,1);
s1(1:8)= 1;
s2(2:5)=[1 2 3 4];
s3= conv(s1,s2);            % compute using the conv routine -> linear
figure(1)
subplot(4,1,1); plot(s1,'o')
subplot(4,1,2); plot(s2,'o')
subplot(4,1,3); plot(s3,'o')

fs1= (fft((s1)));
fs2= (fft((s2)));
fs4= fs1.*fs2;
s4= abs((ifft((fs4))));
subplot(4,1,4); plot(s4,'o');
pause

% case when both are different
N1= 16;
s1= zeros(N1,1);
s2= zeros(N1,1);
s1(2:15)= 1;
s2(2:5)=[1 2 3 4];
s3= conv(s1,s2);            % compute using the conv routine -> linear
figure(2)
subplot(5,1,1); plot(s1,'o')
subplot(5,1,2); plot(s2,'o')
subplot(5,1,3); plot(s3,'o')

fs1= (fft((s1)));
fs2= (fft((s2)));
fs4= fs1.*fs2;
s4= abs((ifft((fs4))));
subplot(5,1,4); plot(s4,'o');


%%%% Compute the correct linear conv using fft
ss1= zeros(N1*2,1);
ss2= zeros(N1*2,1);
ss1(1:N1)= s1;
ss2(1:N1)= s2;
fss1= (fft((ss1)));
fss2= (fft((ss2)));
fss4= fss1.*fss2;
ss4= abs((ifft((fss4))));
subplot(5,1,5); plot(ss4,'o');
pause

%%%%%%%  Sampling related issues
%%% Leakage
N2= 32;
ind1= 1:N2;
ss1= sin(2*pi*ind1/8);
ss2= sin(2*pi*ind1/9);
fss1= fftshift(fft(fftshift(ss1)));
fss2= fftshift(fft(fftshift(ss2)));
figure(4)
title('Spectral Leakage Example');
subplot(2,2,1); plot(ss1);
subplot(2,2,2); plot(abs(fss1),'o');
subplot(2,2,3); plot(ss2);
subplot(2,2,4); plot(abs(fss2),'o');
pause

%%%% Gibbs ringing
N3= 64;
sss1= zeros(N3,1);
sss1(1:N3/8)= 1;
fsss1= fft(sss1);
fsss2= fsss1;
fsss2(N3*3./8:N3*5/8.)=0;
sss2= ifft(fsss2);
figure(3);
title('Gibbs Ringing Example');
subplot(2,2,1); plot(sss1);
subplot(2,2,2); plot(fftshift(abs(fsss1)));
subplot(2,2,3); plot(real(sss2));
subplot(2,2,4); plot(fftshift(abs(fsss2)));




