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


n=34;
nfft=4096;
sysclock=250

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

% data file source:
% http://forums.parallax.com/discussion/comment/1455002/#Comment_1455002
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];


% other filters previously tested 

delta70taps=[0 2 5 -13 -16 -18 -51 -53 -56 64 67 69];
delta70=zeros(71,1);
delta70( 1+delta70taps(find(delta70taps >=0)) )=1;
delta70( 1-delta70taps(find(delta70taps <0 )) )=-1; % it is assumed that 0 is positive

delta47taps=[0 2 5 -11 -14 -16 -30 -32 -35 41 44 46];
delta47=zeros(71,1);
delta47( 1+delta47taps(find(delta47taps >=0)) )=1;
delta47( 1-delta47taps(find(delta47taps <0 )) )=-1; % it is assumed that 0 is positive

delta30taps=[0 1 3 -11 -13 -14 -15 -16 -18 26 28 29];
delta30=zeros(71,1);
delta30( 1+delta30taps(find(delta30taps >=0)) )=1;
delta30( 1-delta30taps(find(delta30taps <0 )) )=-1; % it is assumed that 0 is positive


hann30=cumsum(cumsum( delta30));
tukey47=cumsum(cumsum( delta47));
tukey70=cumsum(cumsum( delta70));

%http://forums.parallax.com/discussion/comment/1457896/#Comment_1457896
blackman22=[   1, 2, 4, 6, 8,11,14,17,20,22,23	...
   23,22,20,17,14,11, 8, 6, 4, 2, 1 ]';


cma67=[  1, 3, 6,10,15,21,28,36,45,53,60,66,71,75,78  ...		[15] Ramp up   =  568 \
         80*ones(1,37)   ...							[37] Plateau   = 2960  > 4096
         78,75,71,66,60,53,45,36,28,21,15,10, 6, 3, 1]'; 



a=8 
b=47
m=1

t0=1;
t1 = t0+a;
t2 = t1+a;
t3 = t2+b;
t4 = t3+a;
t5 = t4+a;

deltatest=zeros(80,1);
deltatest([t0 t2])=1;
deltatest([t1])=-2;
deltatest([t3 t5])=-1;
deltatest([t4])=2;
deltataps=(cumsum(cumsum(cumsum(deltatest ))));
deltasum=sum( deltataps )

		

% Find factor combinations for the desired sum
ncma=3;
facts=[];

for s=[255 510 511 1020:1023]
  facts=[facts ; factorgroups(s,ncma) ];	
end
[fs fi]=sort(sum(facts,2));
facts=facts( fi,:); % sort by length
siglist=[];
for r=1:size(facts,1)
	siglist=[siglist ; 'cmafilter(' num2str(facts(r,1)) ',' num2str(facts(r,2)) ','  num2str(facts(r,3)) ')' ];
end


traces = [];
highbits = [];
adcmean = [];
adcqmean = [];
adcstd = [];
adcqstd = [];
sums = [];
adcscale = [];
hfnoise = [];
vpp = [];
lengths = [];
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)) ];
	sums = [sums ; sum(impulse) ];
	hfnoise = [hfnoise ; sum( psd(sumstart:sumend).^2) ];
	lengths = [lengths ; length(find(impulse)) ];
	thissort = sort(thisadc);
	highs=thissort(floor(length(thissort)*0.99),:);
	lows=thissort(floor(length(thissort)*0.01),:);
	vpp = [vpp ; (highs-lows)/thisscale];
end

f3db=max((traces(1:nfft/2,:)>-3).*((1:nfft/2)' /nfft*sysclock ))';
f6db=max((traces(1:nfft/2,:)>-6).*((1:nfft/2)' /nfft*sysclock ))';


% Assemble the filter comparison chart
chartdata=[mean(vpp,2) mean(adcstd,2) mean(hfnoise,2) lengths sums f3db f6db];

%disp(['WindowFunction  UnquantizedVpp98%  UnquantizedStdDev  HFNoisePower  HighBitCount' ]);
disp(['WindowFn                Vpp98%         StdDev   HFNoisePower         Length            Sum      freq -3dB      freq -6dB' ]);
disp([ siglist repmat(' ',size(siglist,1),1) num2str(chartdata) ]);

% analytical CIC decimator spectrum response
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
