%Matlab file used for class demoS
clf
clear all
set(0,'defaultaxesfontsize',20);
Fs=20; Ts=1/Fs;
t=[-Fs/2:Ts:Fs/2];
std=0.0975; %std dev for showing linear phase
%
t0=0;
%t0=1/20;
eps=10^(-7);
a=1;
%
%
% Tsa=6; 
% Tr=0.05;
% xst=zeros(size(t));
% ts=-10:Tsa:10;
% for k=1:length(ts),
% xst=xst+(1/Tr)*rect((t-ts(k))/Tr);
% end
% x=exp(-pi*((a*(t-t0)).^2));
% x=x.*xst;
%
% h=Tsa*sin(pi*(1/Tsa)*(t-eps))./(pi*(t-eps)); df=0.05;
% h=(1/(2*df))*h.*(sin(2*pi*df*(t-eps))./(pi*(t-eps)));
%
%
%x=(t-t0).*exp(-pi*((a*(t-t0)).^2));
x=cos(2*pi*2*t);
%x=cos(2*pi*5*(t-t0)).*exp(-pi*((a*(t-t0)).^2));
%x=cos(2*pi*5*(t-t0)).*x;
%x=exp(-j*2*pi*5*(t-t0)).*x;
%x=cos(2*pi*4*(t-t0))+0.75*cos(2*pi*7*(t-t0));
%
hf=sin(2*5*pi*(t-eps))./(pi*(t-eps));
df=0.1; hf=(1/(2*df))*hf.*(sin(2*pi*df*(t-eps))./(pi*(t-eps)));
hbt=2*hf.*sin(2*5*pi*t);
h=hbt;
xt=Ts*conv(x,hbt); xt=xt(201:end-200);
xd=x;
%x=x+j*xt;
x=xd.*cos(2*pi*5*t)-xt.*sin(2*pi*5*t);
%
%x=sin(2*pi*5*(t-eps))./(pi*(t-eps));
%df=0.1; x=(1/(2*df))*x.*(sin(2*pi*df*(t-eps))./(pi*(t-eps)));
%x=x.*cos(2*pi*5*t);
%
%x=sin(2*pi*5*(t-t0-eps))./(pi*(t-t0-eps));
%df=0.1; x=(2*df)*x.*(sin(2*pi*df*(t-eps))./(pi*(t-eps)));
%y=(1/Ts)*diff(x); yt=y;
%x=cos(2*pi*4*t)+0.75*cos(2*pi*7*t);
%x=rect((t-t0)/1);
%
% h=sin(5*pi*1*(t-eps))./(pi*(t-eps)); df=0.1;
% h=(1/(2*df))*h.*(sin(2*pi*df*(t-eps))./(pi*(t-eps)));
% h=2*h.*cos(2*pi*5*t);
%
%x=cos(2*pi*2*(t-t0))+cos(2*pi*5*(t-t0))+cos(2*pi*8*(t-t0));
%h=2*h.*cos(2*pi*5*t);
%h=3*2*h.*cos(2*pi*1.5*t);
%h=2*2*h.*cos(2*pi*4.5*t);
%h=h+2*h.*cos(2*pi*7.5*t);
%h=exp(j*2*pi*3*t).*sin(2*pi*3*(t-eps))./(pi*(t-eps));
%yz=Ts*conv([y 0],h); yt=y; y=yz;
y=Ts*conv(x,h);
%
%
N=length(t);
df=Fs/(N-1);
freqomega=-Fs/2:df:Fs/2;
%x=h;
Ns=length(freqomega);
for in=1:Ns,
    Xfft(1,in)=Ts*(x*exp(-j*2*pi*freqomega(in)*t).');
    %Xsfft(1,in)=Ts*(xst*exp(-j*2*pi*freqomega(in)*t).');
    Hfft(1,in)=Ts*(h*exp(-j*2*pi*freqomega(in)*t).');
    %Hfft(1,in)=Ts*(hbt*exp(-j*2*pi*freqomega(in)*t).');
    %Xfft(1,in)=Ts*sum(x.*exp(-j*2*pi*freqomega(in)*t));
    %Xfft(1,in)=Ts*[yt 0]*exp(-j*2*pi*freqomega(in)*t).';
    %Yfft(1,in)=Ts*y(201:end-200)*exp(-j*2*pi*freqomega(in)*t).';
    Yfft(1,in)=Ts*xd*exp(-j*2*pi*freqomega(in)*t).';
end
EnergyTime=sum((abs(x).^2)')
EnergyFreq=sum((abs(Xfft).^2)')
Ratio=EnergyFreq/EnergyTime
%phoff=exp(j*200*2*pi*freqomega); Xph=angle(Xfft.*phoff);
%Xph=angle(Xfft);
Xph_p=unwrap(angle(Xfft(201:end)));
Xph=[-Xph_p(end:-1:2) Xph_p];
Hph=(angle(Hfft(1:end)));
figure(1)
subplot(211)
x=real(x);
plot(t,x,'LineWidth',3);
axis([-Fs/2,Fs/2 min(x) max(x)]);
set(gca,'XTick',-Fs/2:1:Fs/2);
%plot(t,[yt 0],'LineWidth',3);
%axis([-Fs/2,Fs/2 min(yt) max(yt)]);
title('Time-Domain Signal'),...
xlabel('Time (s)'),ylabel('Amplitude'), grid
subplot(212)
title('Magnitude of Fourier Transform')
plot(freqomega,abs(Xfft),'-r','LineWidth',3);
axis([-Fs/2 Fs/2 0 max(abs(Xfft))]);
set(gca,'XTick',-Fs/2:1:Fs/2);
%plot(freqomega,imag(Xfft),'-r','LineWidth',3);
%axis([-Fs/2 Fs/2 -max(abs(Xfft)) max(abs(Xfft))]);
title('Magnitude of Fourier Transform'),...
xlabel('Frequency (Hz)'), ylabel('Magnitude'), grid
% 
figure(2)
subplot(211)
z=(y(201:end-200));
plot(t,y(201:end-200),'LineWidth',3);
axis([-Fs/2,Fs/2 min(z) max(z)]);
%plot(t,[yt 0],'LineWidth',3);
%axis([-Fs/2,Fs/2 min(yt) max(yt)]);
title('Filtered Time-Domain Signal'),...
xlabel('Time (s)'),ylabel('Amplitude'), grid
subplot(212)
title('Magnitude of Fourier Transform')
plot(freqomega,abs(Yfft),'-r','LineWidth',3);
axis([-Fs/2 Fs/2 0 max(abs(Yfft))]);
set(gca,'XTick',-Fs/2:1:Fs/2);
%plot(freqomega,imag(Xfft),'-r','LineWidth',3);
%axis([-Fs/2 Fs/2 -max(abs(Xfft)) max(abs(Xfft))]);
title('Magnitude of Fourier Transform'),...
xlabel('Frequency (Hz)'), ylabel('Magnitude'), grid
%
figure(3)
subplot(211)
plot(t,h,'LineWidth',3);
axis([-Fs/2,Fs/2 min(h) max(h)]);
set(gca,'XTick',-Fs/2:1:Fs/2);
%plot(t,[yt 0],'LineWidth',3);
%axis([-Fs/2,Fs/2 min(yt) max(yt)]);
title('Time-Domain Impulse'),...
xlabel('Time (s)'),ylabel('Amplitude'), grid
subplot(212)
title('Magnitude of Fourier Transform')
plot(freqomega,abs(Hfft),'-r','LineWidth',3);
axis([-Fs/2 Fs/2 0 max(abs(Hfft))]);
set(gca,'XTick',-Fs/2:1:Fs/2);
%plot(freqomega,imag(Xfft),'-r','LineWidth',3);
%axis([-Fs/2 Fs/2 -max(abs(Xfft)) max(abs(Xfft))]);
title('Magnitude of Fourier Transform'),...
xlabel('Frequency (Hz)'), ylabel('Magnitude'), grid
% figure(4)
% subplot(211)
% plot(t,xst,'LineWidth',3);
% axis([-Fs/2,Fs/2 min(xst) max(xst)]);
% set(gca,'XTick',-Fs/2:1:Fs/2);
% %plot(t,[yt 0],'LineWidth',3);
% %axis([-Fs/2,Fs/2 min(yt) max(yt)]);
% title('Time-Domain Sampling Signal'),...
% xlabel('Time (s)'),ylabel('Amplitude'), grid
% subplot(212)
% title('Magnitude of Fourier Transform')
% plot(freqomega,abs(Xsfft),'-r','LineWidth',3);
% axis([-Fs/2 Fs/2 0 max(abs(Xsfft))]);
% set(gca,'XTick',-Fs/2:1:Fs/2);
% %plot(freqomega,imag(Xfft),'-r','LineWidth',3);
% %axis([-Fs/2 Fs/2 -max(abs(Xfft)) max(abs(Xfft))]);
% title('Magnitude of Fourier Transform'),...
% xlabel('Frequency (Hz)'), ylabel('Magnitude'), grid
% 
%
figure(4)
title('Phase of Fourier Transform')
plot(freqomega,Hph,'LineWidth',3);
%axis([-Fs/2 Fs/2 -120 120]);
axis([-Fs/2 Fs/2 -pi pi]);
set(gca,'XTick',-Fs/2:1:Fs/2);
%axis([-pi pi -pi pi]);
xlabel('Frequency (Hz)'), ylabel('Phase'), grid
slope=(1/(2*pi))*Ts*sum(diff(Xph))