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


n=34;
nfft=4096;
sysclock=250

sumstart=floor(nfft/20);
sumend=floor(nfft/2);

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),[nfft 1]);

ham=zeros(nfft,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);

%http://forums.parallax.com/discussion/comment/1455621/#Comment_1455621
tukey1732=[01,03,05,07,10,13,16,19,22,25,27,29,31 ...		'13 up, up/down sum = 208 * 2
		32*ones(9,1)',31,32*ones(9,1)' ...		'19 top, top sum = 607
		31,29,27,25,22,19,16,13,10,07,05,03,01]';	%'13 down, total sum = 1023 (>>2 = 255)

tukey1732t=[01,03,05,07,10,13,16,19,22,25,27,29,31 ...		'13 up, up/down sum = 208 * 2
		32*ones(3,1)' ...		'19 top, top sum = 607
		31,29,27,25,22,19,16,13,10,07,05,03,01]';	%'13 down, total sum = 1023 (>>2 = 255)

tukey1732p=[01,03,05,07,10,13,16,19,22,25,27,29,31 ...		'13 up, up/down sum = 208 * 2
		32,33,32 ...		'19 top, top sum = 607
		31,29,27,25,22,19,16,13,10,07,05,03,01]';	%'13 down, total sum = 1023 (>>2 = 255)

hann1732=[01,03,05,07,10,13,16,19,22,25,27,29,31 ...		'13 up, up/down sum = 208 * 2
		32*ones(1,1)' ...		'19 top, top sum = 607
		31,29,27,25,22,19,16,13,10,07,05,03,01]';	%'13 down, total sum = 1023 (>>2 = 255)

tukey1332=[	01,02,05,08,12,16,20,24,27,30,31  ...		'11 up, up/down sum = 176 * 2
		32*ones(1,10),31,32*ones(1,10)	 ...			'21 top, top sum = 671
		31,30,27,24,20,16,12,08,05,02,01]'; %		'11 down, total sum = 1023 (>>2 = 255)

tukey1633=[	01,02,04,06,09,12,15,18,21,24,27,29,31,32 ...	'14 up, up/down sum = 231 * 2 = 462
	33*ones(1,17)					...	'17 top, top sum = 561
	32,31,29,27,24,21,18,15,12,09,06,04,02,01 ]';%	'14 down, total sum = 1023 (>>2 = 255)
tri=[1:23 22:-1:1]';

%http://forums.parallax.com/discussion/comment/1456266/#Comment_1456266
tukey17_13_32 = [01,03,05,07,10,13,16,19,22,25,27,29,31	...	'13 up, up/down sum = 208 * 2
		32*ones(1,9),31,32*ones(1,9)				...	'19 top, top sum = 607
		31,29,27,25,22,19,16,13,10,07,05,03,01	]';	%'13 down, total sum = 1023 (>>2 = 255)


tukey15_13_64 = [01,03,07,12,18,25,32,39,46,52,57,61,63	...	'13 up, up/down sum = 416 * 2
		64*ones(1,9),63,64*ones(1,9)	...			'19 top, top sum = 1215
		63,61,57,52,46,39,32,25,18,12,07,03,01	]';	%'13 down, total sum = 2047 (>>3 = 255)

tukey17_13_32_x2 = [02,06,10,14,20,26,32,38,44,50,54,58,62  ...		'13 up, up/down sum = 416 * 2
		64*ones(1,9),62,64*ones(1,9)			...	'19 top, top sum = 1214
		62,58,54,50,44,38,32,26,20,14,10,06,02	]';	%'13 down, total sum = 2046 (>>3 = 255)

tukey17_13_64 =[ 02,05,09,14,20,26,32,38,44,50,55,59,62	...	'13 up, up/down sum = 416 * 2
		64*ones(1,9),63,64*ones(1,9)		...		...	'19 top, top sum = 1215
		62,59,55,50,44,38,32,26,20,14,09,05,02	]';	%'13 down, total sum = 2047 (>>3 = 255)


tukey17_15_64 = [01,02,05,09,14,20,26,32,38,44,50,55,59,62,63	... '15 up, up/down sum = 480 * 2
		64*ones(1,8),63,64*ones(1,8)			... 		'17 top, top sum = 1087
		63,62,59,55,50,44,38,32,26,20,14,09,05,02,01 ]'; % '15 down, total sum = 2047 (>>3 = 255)


vn=34
vor=sin(0.5*pi*sin((0:vn-1 + 0.5)/vn*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;


cmabase = conv(conv(conv(ones(19,1),ones(15,1),'full'),ones(12,1)),ones(8,1));
cmatest=round(cmabase/(sum(cmabase)/1023));
cmasum=sum(cmatest);
cmalen=length(find(cmatest));
cma9999 = conv(conv(conv(ones(9,1),ones(9,1),'full'),ones(9,1)),ones(9,1));
cmalen=length(cmatest)
cma121212 = conv(conv(ones(12,1),ones(12,1),'full'),ones(12,1));
addendscma=sum(sum(bsxfun(@bitget,cmatest,1:32)));



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

adcchan=1
if 0
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');
end


siglist=['tukey17_13_32' ; 'tukey15_13_64' ; 'tukey17_13_32_x2' ; 'tukey17_13_64' ; 'tukey17_15_64'];


traces = [];
highbits = [];
adcmean = [];
adcqmean = [];
adcstd = [];
adcqstd = [];
adcscale = [];
hfnoise = [];
vpp = [];
for c=1:size(siglist,1)
 	thisname=siglist(c,:)
	impulse=eval( siglist(c,:) );
	psd=abs(fft(impulse,nfft));
	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) ];

	highbits = [highbits ; sum(sum(bsxfun(@bitget,impulse,1:32))) ];

	thisadc=conv2(abit,impulse,'valid');
	thisscale=round( sum(impulse)/255 );  % try to calculate the scale to get max=255
	adcscale=[adcscale ; thisscale];
	adcmean = [adcmean ; mean(thisadc/thisscale) ];
	adcqmean = [adcqmean ; mean(floor(thisadc/thisscale)) ];
	adcstd = [adcstd ; std(thisadc/thisscale) ];
	adcqstd = [adcqstd ; std(floor(thisadc/thisscale)) ];
	
	hfnoise = [hfnoise ; sum( psd(sumstart:sumend).^2) ];

	thissort = sort(thisadc);
	highs=thissort(floor(length(thissort)*0.99),:);
	lows=thissort(floor(length(thissort)*0.01),:);
	vpp = [vpp ; (highs-lows)/thisscale];
end

chartdata=[mean(vpp,2) mean(adcstd,2) mean(hfnoise,2) highbits];

disp(['WindowFunction  UnquantizedVpp98%  UnquantizedStdDev  HFNoisePower  HighBitCount' ]);
disp([ siglist repmat(' ',size(siglist,1),1) num2str(chartdata) ]);

M=10
R=1
N=4
f=(0:nfft-1)'/nfft;
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 1
plotrange=floor(nfft/2);
plot((1:plotrange')/(plotrange*2)*sysclock,traces(1:plotrange,:),'-')
ylim([-120 0]);

title('Window Filter Frequency Response')
legend(siglist);
ylabel('dB');

if sysclock==1
  xlabel('relative frequency')
else
  xlabel(['MHz, sysclock=' num2str(sysclock)])
end
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 0
	hold off
  plot(round(cmatest(3:end-2)/29.61),'x')
  hold on
  y=[0 20];
  confstr='r'
  plot( [1 1]*(0.5),y,confstr)
  plot( [1 1]*(8.5),y,confstr)
  plot( [1 1]*(16.5),y,confstr)
  plot( [1 1]*(24.5),y,confstr)
  plot( [1 1]*(32.5),y,confstr)
  hold off

end
