%
% (C) 2018 Andy Walls <address@hidden>
%

% The X310 DDC chain has 4 Integrator Stages  and 4 single delay Comb stages
% followed by a decimation stage.  We'll model that here, using the
% equations in https://www.altera.com/content/dam/altera-www/global/en_US/pdfs/literature/an/an455.pdf

M = 1; % differential delay
N = 4; % stages
R = 4; % decimation rate to go from 200 Msps down to 50 Msps
d = ones(1, R*M); % FIR filter coefficients before raising to N
                  % Since it is like a cascade of N identical filters
[hd, wd] = freqz(d, [1], R*8*1024);
hd = hd .^ N;
hd = hd/hd(1);
wdturn = [wd(1:(R*8*1024/R))' wd((R*8*1024/R):-1:1)' ...
          wd(1:(R*8*1024/R))' wd((R*8*1024/R):-1:1)'] * 50/pi; % Fs = 200, fold into 0-50 multiple times

% The X310 DDC uses up to 3 x% Nyquist Half-band filters

% Coefficients from
% uhd/fpga-src/usrp3/top/x300/coregen_dsp/hb47.coe
% Generated with round((2^18-2)*halfgen_test(.21,12,1))
%  83% effective BW, 64dB to 85dB rolloff

a = [   -62, 0,    194, 0,  -440, 0,   855, 0,  -1505, 0,  2478, 0, ...
      -3900, 0,   5990, 0, -9187, 0, 14632, 0, -26536, 0, 83009, ...
     131071, ...
      83009, 0, -26536, 0, 14632, 0, -9187, 0,   5990, 0, -3900, 0, ...
       2478, 0,  -1505, 0,   855, 0,  -440, 0,    194, 0, -62];
b = a;

% Coefficients from
% uhd/fpga-src/usrp3/top/x300/coregen_dsp/hb63.coe
% Generated with round((2^18-2)*halfgen_test(.22,16,1))
%  88% effective BW, 64dB to 87dB rolloff
c = [ -35, 0, 95, 0, -195, 0, 352, 0, -582, 0, 907, 0, -1354, 0, 1953, 0, ...
      -2751, 0, 3813, 0, -5249, 0, 7264, 0, -10296, 0, 15494, 0, -27083, 0, ...
      83196, 131071, 83196, ...
      0, -27083, 0, 15494, 0, -10296, 0, 7264, 0, -5249, 0, 3813, 0, -2751, ...
      0, 1953, 0, -1354, 0, 907, 0, -582, 0, 352, 0, -195, 0, 95, 0, -35];

[ha,wa] = freqz(a, [1], 8*1024);
ha = ha/ha(1);
[hb,wb] = freqz(b, [1], 4*1024);
hb = hb/hb(1);
[hc,wc] = freqz(c, [1], 2*1024);
hc = hc/hc(1);

% Fold freq axis to show aliasing into final passband
% Scale to original 50 Msps Fs
waturn = [wa(1:(8*1024/2))' wa((8*1024/2):-1:1)'] * 25/pi; % Fs = 50, fold 25-50 onto 25-0
wbturn = [wb(1:(4*1024/2))' wb((4*1024/2):-1:1)'] * 12.5/pi; % Fs = 25, fold 12.5-25 onto 12.5-0
wcturn = [wc(1:(2*1024/2))' wc((2*1024/2):-1:1)'] * 6.25/pi; % Fs = 12.5, fold 6.25-12.5 onto 6.25-0

clf;
plot( ...
 wdturn, 20*log10(abs(hd)), ';CIC filter 200 Msps folded to 50 Msps;', ...
 waturn, 20*log10(abs(ha)), ';1st HB filter 50 Msps folded to 25 Msps;', ...
 wbturn, 20*log10(abs(hb)), ';2nd HB filter 25 Msps folded to 12.5 Msps;', ...
 wcturn, 20*log10(abs(hc)), ';3rd HB filter 12.5 Msps folded to 6.25 Msps;');
title('X310 DDC CIC & Default Halfband Filters Alias Folding Performance');
ylabel('Gain (dB)');
xlabel('Frequency (MHz)');
grid on;
hold off;

% Passband droop at the final fold
% At Fs/2 = 6.25 MHz
worst_pass = ...
       20*log10(abs(hd(R*8*1024/32))) ...
     + 20*log10(abs(ha(8*1024/8))) ...
     + 20*log10(abs(hb(4*1024/4))) ...
     + 20*log10(abs(hc(2*1024/2)))
ideal_pass = ...
       20*log10(abs(hd(1))) ...
     + 20*log10(abs(ha(1))) ...
     + 20*log10(abs(hb(1))) ...
     + 20*log10(abs(hc(1)))
worst_droop = ideal_pass - worst_pass