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

n=34;
t=4096;

gtext=textread('gio_bitstream_8ch_64000clocks.txt', '%s');
vtext=textread('vio_bitstream_8ch_64000clocks.txt', '%s');
gnum=hex2dec (gtext);
vnum=hex2dec (vtext);
gbit=[];
vbit=[];
for j=1:8
 gbit=[ gbit bitget(gnum,j)];
 vbit=[ vbit bitget(vnum,j)];
end
abit=[gbit vbit];


rec=resize(ones(n,1),[t 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);


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);



vn=34
vor=resize(sin(0.5*pi*sin((0:vn-1 + 0.5)/vn*pi).^2 )',[t 1]);



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=10 %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;
for offs=1:10:91
 %offs
  sinc3(1+offs:length(sinc3base)+offs)=sinc3(1+offs:length(sinc3base)+offs)+sinc3base;
end
sinc3=sinc3/max(sinc3);

sinc888base=conv(conv(ones(8,1),ones(8,1)),ones(8,1));
%sinc888base=sinc888base/sum(sinc888base);
sinc888=resize(sinc888base,[t,1]);

sinc1688base=conv(conv(ones(16,1),ones(8,1)),ones(8,1));
%sinc1688base=sinc1688base/sum(sinc1688base);
sinc1688=resize(sinc1688base,[t,1]);


cmatest = conv(conv(conv(ones(12,1),ones(10,1),'full'),ones(9,1)),ones(7,1));
cma9999 = conv(conv(conv(ones(9,1),ones(9,1),'full'),ones(9,1)),ones(9,1));
cmalen=length(cma161208)
cma121212 = conv(conv(ones(12,1),ones(12,1),'full'),ones(12,1));


testfilter=cmatest;
testlength=length(testfilter)
adccma=conv2(abit,testfilter/sum(testfilter),'valid');
adcrect=conv2(abit,ones(testlength,1)/testlength,'valid');

adcchan=1
plot([adcrect(:,[0 8]+adcchan) adccma(:,[0 8]+adcchan)]*256)
title(['Filtered ADC values, filter length=' num2str(testlength) ])
legend('Gio Rectangular','Vio Rectangular','Gio Sinc','Vio Sinc')
print('plotadc1.png');

signals=[ rec resize(cmatest,[t 1]) resize(cma121212,[t 1]) ];%resize(cma9999,[t 1])];

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=10
R=1
N=4
f=(0:t-1)'/t;
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) ];

if 0
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')
legend('Rectangular/Sinc^1','Cascaded Moving Average 12*10*8*6','Cascaded Moving Average 12*12*12','Cascaded Moving Average 9*9*9*9','AlmostGaussian/Sinc^3','Vorbis4096','Tukey17/32','Tukey17rep2/32')
ylabel('dB');
xlabel('relative frequency')
print('plotmf.png');
end

if 0
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('plotmt.png');
end

if 1
  
end
