%% --------------------------------------------------------
% Single-section DF-II IIR filter simulator
% Floating-point / fixed-point 
% DRB November 2012
% --------------------------------------------------------
% User parameters below
% --------------------------------------------------------
Fs = 16000;       % Sampling Frequency
Fstop1 = 3200;    % First Stopband Frequency
Fpass1 = 4000;    % First Passband Frequency
Fpass2 = 5000;    % Second Passband Frequency
Fstop2 = 6000;    % Second Stopband Frequency
Astop1 = 40;      % First Stopband Attenuation (dB)
Apass  = 0.5;     % Passband Ripple (dB)
Astop2 = 40;      % Second Stopband Attenuation (dB)
match  = 'both';  % Band to match exactly
T = 10;           % duration of white noise
Mb = 19;          % number of frac bits in numerator
Ma = 12;          % number of frac bits in denominator
Mu = 20;          % number of fractional bits in intermediate values 
maxu = 0;         % track maximum intermediate value (float)
maxui = 0;        % track maximum intermediate value (int)
maxy = 0;         % track maximum output (float)
maxyi = 0;        % track maximum output (int)
Px = -5;          % shift plan: input right pre-shift
Pa = 0;           % shift plan: feedback coefficient right pre-shift
Pufb = 11;        % shift plan: intermediate value right pre-shift for feedback
Pua = 1;          % shift plan: post feedback product right shift
Pb = 0;           % shift plan: feedforward coefficient right pre-shift
Puff = 17;        % shift plan: intermediate value right pre-shift for feedforward
Pub = 2;          % shift plan: post feedforward product right shift
Pfinal = 5;       % shift plan: pre-shift output (if necessary)
% --------------------------------------------------------

%% generate 99% full-scale white noise (single channel, unquantized)
x = 0.99*(2*rand(1,T*Fs)-1);

%% design floating point filter
h  = fdesign.bandpass(Fstop1, Fpass1, Fpass2, Fstop2, Astop1, Apass, Astop2, Fs);
Hd = design(h, 'ellip', 'MatchExactly', match);
[b, a] = tf(Hd);
[H,ff] = freqz(b,a,512,Fs);

%% do DF-II single-section floating point filter
u = zeros(1,length(b));     % intermediate values
y = zeros(1,length(x));     % outputs
for n=1:length(x),
    u(2:end) = u(1:end-1);  % move u-buffer down
    u(1) = x(n)-sum(u(2:end).*a(2:end));
    if abs(u(1)>maxu),
        maxu = u(1);
    end
    y(n) = sum(u.*b);
end

disp(['Maximum float intermediate value: ' num2str(maxu)]);
disp(['Maximum float output value: ' num2str(max(abs(y)))]);

%% Quantize coefficients
bi = round(b*2^Mb);         % integer values
ai = round(a*2^Ma);         % integer values
bfp = bi/2^Mb;             % fixed-point values
afp = ai/2^Ma;             % fixed-point values

%% do simplified DF-II single-section fixed point filter
% this code just shows the effect of coefficient quantization
% there is no shifting here
ufp = zeros(1,length(b));   % fixed-point intermediate values
yfp = zeros(1,length(x));   % fixed-point outputs
for n=1:length(x),
    ufp(2:end) = ufp(1:end-1);  % move u-buffer down
    ufp(1) = x(n)-sum(ufp(2:end).*afp(2:end));
    yfp(n) = sum(ufp.*bfp);
end

%% do actual DF-II single-section fixed point filter with shifts
ui = zeros(1,length(b));   % integer intermediate values
yi = zeros(1,length(x));   % integer outputs
xi = round(x*2^15);        % integer inputs (between -32768 and 32767)
feedbackprodi = zeros(1,length(ai));
feedforwardprodi = zeros(1,length(bi));
for n=1:length(xi),
    ui(2:end) = ui(1:end-1);  % move u-buffer down
    tmp = fix(xi(n)/2^Px);    % pre-shift x
    for k=2:length(ai),
        feedbackprodi(k) = fix(ui(k)/2^Pufb)*fix(ai(k)/2^Pa);
        tmp = tmp-fix(feedbackprodi(k)/2^Pua);
    end
    ui(1) = tmp;  
    if abs(ui(1)>maxui),
        maxui = ui(1);
    end
    for k=1:length(bi),
        feedforwardprodi(k) = fix(ui(k)/2^Puff)*fix(bi(k)/2^Pb);
        yi(n) = yi(n) + fix(feedforwardprodi(k)/2^Pub);
    end
    yi(n) = fix(yi(n)/2^Pfinal);
end

disp(['Maximum integer intermediate value: ' num2str(maxui)]);
disp(['Maximum integer output value: ' num2str(max(abs(yi)))]);

%% plot PSDs
[Px,fx] = pwelch(x,1024,512,1024,Fs); % estimate psd of input
[Py,fy] = pwelch(y,1024,512,1024,Fs); % estimate psd of output
[Pyfp,fyfp] = pwelch(yfp,1024,512,1024,Fs); % estimate psd of fixed-point output
[Pyi,fyi] = pwelch(yi/2^15,1024,512,1024,Fs); % estimate psd of integer output
plot(fx/1e3,10*log10(Px)+44,...
     fy/1e3,10*log10(Py)+44,...
     fyfp/1e3,10*log10(Pyfp)+44,...  
     fyi/1e3,10*log10(Pyi)+44,...
     ff/1e3,10*log10(abs(H).^2)); % plot psd in dB
grid on
legend('input','floating-point filter output','fixed-point filter output',...
       'integer filter output','theoretical');
xlabel('frequency (kHz)');
ylabel('magnitude response (dB)'); 