
% A simple analysis tool for windowing functions
%
%http://forums.parallax.com/discussion/169298/adc-sampling-breakthrough
% SaucySoliton at forums.parallax.com
% 11/23/2018
% Made for GNU Octave

n=4096;
t=n*3;


rec=zeros(t,1);
rec(1:n)=1;

ham=zeros(t,1);
ham(1:n)=hamming(n);

trap32=rec;
trap32(1:32)=(1:32)/32;
trap32(n-31:n)=1-((1:32)/32);

trap64=rec;
trap64(1:64)=(1:64)/64;
trap64(n-63:n)=1-((1:64)/64);

trap128=rec;
trap128(1:128)=(1:128)/128;
trap128(n-127:n)=1-((1:128)/128);

trap256=rec;
trap256(1:256)=(1:256)/256;
trap256(n-255:n)=1-((1:256)/256);




seven=rec;
seven(1:586*7)=repmat([1 0 1 1 0 0 0]',586,1);
%seven(1:586*8)=repmat([1 0 1 0 1 0 1 0]',586,1);

%http://forums.parallax.com/discussion/comment/1454753/#Comment_1454753
tukeybase2464=[01 02 03 05 08 11 14 18 22 26 30 34 38 42 46 50 53 56 59 61 62 63 64]'/64; 
tukeybase1664=[1 4 7 12 17 23 29 35 41 47 52 57 60 63 64]'/64;
%http://forums.parallax.com/discussion/comment/1454977/#Comment_1454977
tukeybase1732=[1 3 5  7 10 13 16 19 22 25 27 29 31 32]'/32;
tukey16=rec;
tukey16(1:length(tukeybase1664))=tukeybase1664;
tukey16(n-length(tukeybase1664)+1:n)=flipud(tukeybase1664);
tukey24=rec;
tukey24(1:length(tukeybase2464))=tukeybase2464;
tukey24(n-length(tukeybase2464)+1:n)=flipud(tukeybase2464);

tukey17d=rec;
tukeybase1732d=kron(tukeybase1732,[1 1]');
tukey17d(1:length(tukeybase1732d))=tukeybase1732d;
tukey17d(n-length(tukeybase1732d)+1:n)=flipud(tukeybase1732d);



%tukeyfbase=(1+sin((linspace(-7,7,16)')/16*pi))/2
%tukeyfbase=round(((1+sin((linspace(-14,14,32)')/32*pi))/2)*64)/64
tukeyfbase=((1+sin((linspace(-4,4,8)'+.1)/10*pi))*32)/64
tukeyf=rec;
tukeyf(1:length(tukeyfbase))=tukeyfbase;
tukeyf(n-length(tukeyfbase)+1:n)=flipud(tukeyfbase);


vor=rec;
vor(1:n)=sin(0.5*pi*sin((0:n-1 + 0.5)/n*pi).^2 );


sinc2base=conv(ones(n/2,1),ones(n/2,1),'full');
sinc2base=sinc2base/max(sinc2base);
size(sinc2base);
sinc2=rec*0;
sinc2(1:length(sinc2base))=sinc2base;

len3=round(n/3)
sinc3base=conv(conv(ones(len3,1),ones(len3,1),'full'),ones(len3,1));
sinc3base=sinc3base/max(sinc3base);

size(sinc3base);
sinc3=rec*0;
sinc3(1:length(sinc3base))=sinc3base;


signals=[ rec trap32 trap256 sinc2 sinc3 vor tukey17 tukey17d];

traces = [];
for c=1:size(signals,2)
	psd=abs(fft(signals(:,c)));
	psd=psd/max(psd); % normalize amplitude

	psd=max([ psd [0;psd(1:end-1)] [0;0;psd(1:end-2)] ],[],2); % neighborhood 3 maximum
	%psd=conv(psd,[1 1 1]'/3,'same'); % smoothing
	% without maximum or smoothing, plot oscillates between a value and near zero
	traces = [traces 20*log10(abs(psd)+1e-7) ];
end

M=1
R=1024
N=2
f=(1:t)'/t*n;
cicresponse=abs(sin(pi*M*f)./sin(pi*f/R)).^N ;
lfgain=max(cicresponse)
cicresponse=cicresponse/max(cicresponse);
%traces = [traces  20*log10( cicresponse + 1e-7) ];
plotrange=floor(t/2);

plot((1:plotrange')/(plotrange*2),traces(1:plotrange,:),'-')
ylim([-120 0]);

title('Window Filter Frequency Response, Lenth 4096')
legend('Rectangular/Sinc^1','Trapezoid32','Trapezoid256','Triangular/Sinc^2','AlmostGaussian/Sinc^3','Vorbis4096','Tukey17/32','Tukey17rep2/32')
ylabel('dB');
xlabel('relative frequency')
print('plotdf.png');


plot(signals(1:n,:))
title('Window Filter Impulse Response, Lenth 4096')
legend('Rectangular/Sinc^1','Trapezoid32','Trapezoid256','Triangular/Sinc^2','AlmostGaussian/Sinc^3','Vorbis4096','Tukey17/32','Tukey17rep2/32')
print('plotdt.png');
