clear
% Bresenham fractional divider vs DDS
% First Question: Does Jitter for DDS depend on 
% accumulator width?

%Example: 11 MHz clock= 11 * 5^6 * 2^6
% 310kHz output = 31 * 5^4 *2^4
nn=2*31*11*100;
%nn=100;
% Signal references
% add some phaseshifts to avoid 'exact' zero results

fastsig=    (real(exp(j*(2*pi*(0:nn-1)/(31  )+0.1234))));
slowsig=    (real(exp(j*(2*pi*(0:nn-1)/(1100)+0.1234))));
fastref=sign(real(exp(j*(2*pi*(0:nn-1)/(31  )+0.1234))));
slowref=sign(real(exp(j*(2*pi*(0:nn-1)/(1100)+0.1234))));

spfastsig=fft(fastsig);
spslowsig=fft(slowsig);
spfastref=fft(fastref);
spslowref=fft(slowref);

plot(0:nn-1,fastsig,'b.-',....
     0:nn-1,fastref,'r.-');

plot(0:nn-1,20*log10(abs(spfastsig)),'b.-',....
     0:nn-1,20*log10(abs(spfastref)),'r.-');

plot(0:nn-1,20*log10(abs(spfastref)),'b.-',....
     0:nn-1,20*log10(abs(spslowref)),'r.-');

acc32     =uint64(0);
phaseinc32=uint64(2^32/31);
fastdds32=zeros(size(fastref));

for(k=1:nn)
  acc32=acc32+phaseinc32;
  acc32=bitand(acc32,2^32-1); % too dumb for overflow arith
  fastdds32(k)=2*(bitand(acc32,2^31)>0)-1;
end;
spfastdds32=fft(fastdds32);

plot(0:nn-1,20*log10(abs(spfastref)),'b.-',....
     0:nn-1,20*log10(abs(spfastdds32)),'r.-');
% identical! nn too small to see any difference

BITBIT=20;
acc32     =uint64(0);
phaseinc32=uint64((2^BITBIT)/31);
fastdds16=zeros(size(fastref));

for(k=1:nn)
  acc32=acc32+phaseinc32;
  acc32=bitand(acc32,(2^BITBIT)-1); % too dumb for overflow arith
  fastdds16(k)=2*(bitand(acc32,2^(BITBIT-1))>0)-1;
end;
spfastdds16=fft(fastdds16);

plot(0:nn-1,20*log10(abs(spfastdds16)),'b.-',....
     0:nn-1,20*log10(abs(spfastdds32)),'r.-');
 plot(fastdds16-fastdds32,'.');

 % answer question 1: yes, threshold effect, any width
% broader 19 is ok. Then no difference betwenn fastdds16 and
% fastdds32

%ph=exp(j*2*pi*(0:nn-1)/31)./exp(j*2*pi*(0:nn-1)/1100);

ph=296*31;
fracdiv=zeros(size(fastsig));
for(k=1:nn)
 ph= ph+1*31;
 fracdiv(k)=2*(ph<=((1100*31)/2))-1;
 if(ph>=1100*31) ph=ph-1100*31;end;
end;
spfracdiv=fft(fracdiv);
plot(0:nn-1,20*log10(abs(spslowref)),'b.-',....
     0:nn-1,20*log10(abs(spfracdiv)),'r.-');
plot(1:3000,slowref(1:3000),'.-',1:3000,fracdiv(1:3000),'.-')
sum(abs(fracdiv-slowref))
return

bres1= (0:nn-1)*31


% Bresenham == fractional divider ?
xxx=0:nn-1;
ee=fix((nn/2+(0:nn-1)*deltay)/nn);


return