%Matlab file used for class demo conducted during Session 12
%to illustrate digital subbanding -- transmultiplexers.
%An utterance from a male speaker and one from a
%female speaker are transmultiplexed into different subbands.
%--polyphase implementation
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
%std=1;
a=1/std;
t0=0;
%t0=1/20;
eps=10^(-7);
a=10;
%
x=exp(-pi*((a*(t-t0)).^2));
%x=t.*exp(-pi*((sc*(t-t0)).^2));
%x=cos(2*pi*5*(t-t0)).*exp(-pi*((a*(t-t0)).^2));
%x=exp(j*2*pi*5*(t-t0)).*x;
%
%x=sin(2*pi*6*(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)/2);
%
h=sin(2*pi*7*(t-eps))./(pi*(t-eps));
df=0.1;
h=(1/(2*df))*h.*(sin(2*pi*df*(t-eps))./(pi*(t-eps)));
%h=4*h.*cos(2*pi*4*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);
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).');
    %Xfft(1,in)=Ts*(h*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).';
end
%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];
figure(1)
subplot(211)
plot(t,x,'LineWidth',3);
axis([-Fs/2,Fs/2 min(x) max(x)]);
%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))]);
%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)
plot(t,y(201:end-200),'LineWidth',3);
axis([-Fs/2,Fs/2 min(y) max(y)]);
%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))]);
%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)
title('Phase of Fourier Transform')
plot(freqomega,Xph,'LineWidth',3);
axis([-Fs/2 Fs/2 -120 120]);
%axis([-pi pi -pi pi]);
xlabel('Frequency (Hz)'), ylabel('Phase'), grid
slope=(1/(2*pi))*Ts*sum(diff(Xph))