%spectral leakage demonstration
clear
N = 128;

%assign normalized frequency range, f = k/N
%where f=1 is sampling rate

df=1/N;
%df=1/1024;
df=1/2048;
f=0:df:1/16;

figure(1)
y = ones(N,1);

for k=1:length(f)
s(k)=0;
for l=1:N
s(k) = s(k) + (y(l)*exp(i*2*pi*(l-1)*f(k)));
end
end

spec=abs(s);
subplot(2,1,1)
plot(1:N,y)
title('Rectangular Window; N=128')
axis([0 N-1 0 1.5]);
subplot(2,1,2)
semilogy(f,spec/N);
xlabel('f/r')
axis([0 1/16 1e-6 1]);
%print -adobecset -depsc -r300  rect.ps

pause
figure(2)

y = bartlett(N);
for k=1:length(f)
s(k)=0;
for l=1:N
s(k) = s(k) + (y(l)*exp(i*2*pi*(l-1)*f(k)));
end
end

spec=abs(s);
subplot(2,1,1)
plot(1:N,y)
xlabel('n')
title('Bartlett Window; N=128')
axis([0 N-1 0 1.5]);
subplot(2,1,2)
semilogy(f,spec/sum(y));
xlabel('f/r')
axis([0 1/16 1e-6 1]);
%print -adobecset -depsc -r300  bart.ps

pause
figure(3)

y = kaiser(N,4);
for k=1:length(f)
s(k)=0;
for l=1:N
s(k) = s(k) + (y(l)*exp(i*2*pi*(l-1)*f(k)));
end
end
spec=abs(s);
subplot(2,1,1)
plot(1:N,y)
xlabel('n')
title('Kaiser Window; N=128')
axis([0 N-1 0 1.5]);
subplot(2,1,2)
semilogy(f,spec/sum(y));
axis([0 1/16 1e-6 1]);
xlabel('f/r')
%print -adobecset -depsc -r300  kaiser.ps

pause
figure(4)

y = hamming(N);
for k=1:length(f)
s(k)=0;
for l=1:N
s(k) = s(k) + (y(l)*exp(i*2*pi*(l-1)*f(k)));
end
end
spec=abs(s);
subplot(2,1,1)
plot(1:N,y)
xlabel('n')
title('Hamming Window; N=128')
axis([0 N-1 0 1.5]);
subplot(2,1,2)
semilogy(f,spec/sum(y));
axis([0 1/16 1e-6 1]);
xlabel('f/r')
%print -adobecset -depsc -r300  ham.ps

pause
figure(5)

y = hanning(N);
for k=1:length(f)
s(k)=0;
for l=1:N
s(k) = s(k) + (y(l)*exp(i*2*pi*(l-1)*f(k)));
end
end
spec=abs(s);
subplot(2,1,1)
plot(1:N,y)
xlabel('n')
title('Hanning Window; N=128')
axis([0 N-1 0 1.5]);
subplot(2,1,2)
semilogy(f,spec/sum(y));
axis([0 1/16 1e-6 1]);
xlabel('f/r')
%print -adobecset -depsc -r300  han.ps

load tapers128.out
nt=8;
for j=2:nt+1
taper(:,j-1)=tapers128(:,j);
end
for t=1:nt
y=taper(:,t);
for k=1:length(f)
ms(k,t)=0;
for l=1:N
ms(k,t) = ms(k,t) + (y(l)*exp(i*2*pi*(l-1)*f(k)));
end
end
end
spec=abs(ms);
subplot(2,1,1)
plot(1:N,taper)
xlabel('n')
title('Multitaper Windows; N=128; BW = 1/32')
axis([0 N-1 -2.5 2.5]);
subplot(2,1,2)
semilogy(f,spec/N);
xlabel('f/r')
axis([0 1/16 1e-6 1]);
%print -adobecset -depsc -r300  multi.ps