coeffquant.m

IIR Filter Coefficient Quantization Example ％ written by DRB on 6-Dec-2004 for ECE3703 ％ ％ Looks

% IIR Filter Coefficient Quantization Example
% written by DRB on 6-Dec-2004 for ECE3703
%
% Looks at unquantized filter (ideal), direct form, and cascaded second
% order section (SOS).

% general parameters
Fs=44100;                       % sampling frequency
Wp=[0.3 0.4];                   % edges of passband (frequencies normalized to Fs/2)
Ws=[0.29 0.41];                 % edges of stopbands (frequencies normalized to Fs/2)
Rp=0.2;                         % dB spec for passband
Rs=40;                          % dB spec for stopband
nfft=2^12;                      % number of fft points in magnitude response

% direct form filter parameters
Nbits=15;                       % number of bits for fixed-point quantization of numerator
Dbits=15;                       % number of bits for fixed-point quantization of denominator
Nfrac=15;                       % number of bits in fractional part of numerator
Dfrac=8;                        % number of bits in fractional part of denominator

% SOS filter parameters
Nbits_sos=15;                   % number of bits for fixed-point quantization of numerator
Dbits_sos=15;                   % number of bits for fixed-point quantization of denominator
Nfrac_sos=15;                   % number of bits in fractional part of numerator
Dfrac_sos=14;                   % number of bits in fractional part of denominator

% ------------------------------------------------------------------
% UNQUANTIZED IDEAL FILTER
% ------------------------------------------------------------------

    % The following code designs the unquantized filter
    n=ellipord(Wp,Ws,Rp,Rs);        % get the minimum order of the elliptical filter
    [b,a]=ellip(n,Rp,Rs,Wp,'bandpass');     % design the elliptical filter

    % Plot the magnitude response of the unquantized filter
    figure(1)
    [h,f]=freqz(b,a,nfft,Fs);
    plot(f/1E3,10*log10(abs(h).^2),'LineWidth',2);
    hold on

    % Plot the poles and zeros of the unquantized filter
    figure(2)
    p=roots(a);
    z=roots(b);
    plot(real(p),imag(p),'bx',real(z),imag(z),'bo');
    hold on

% ------------------------------------------------------------------
% QUANTIZED DIRECT FORM FILTER
% ------------------------------------------------------------------
    bq=round(b*2^Nfrac)/2^Nfrac;        % quantized numerator
    aq=round(a*2^Dfrac)/2^Dfrac;        % quantized denominator
    Nmax=sum(2.^[0:Nbits-1])/2^Nfrac;   % upper limit on numerator coefficients
    Nmin=-2^(Nbits-Nfrac);              % lower limit on numberator coefficients
    Dmax=sum(2.^[0:Dbits-1])/2^Dfrac;   % upper limit on denominator coefficients
    Dmin=-2^(Dbits-Dfrac);              % lower limit on denominator coefficients
    for i=1:length(bq),                 % Now correct for any overflow
        if bq(i)>Nmax,                  % (assume saturation)
            bq(i)=Nmax;
        elseif bq(i)<Nmin,
            bq(i)=Nmin;
        end
        if aq(i)>Dmax,
            aq(i)=Dmax;
        elseif aq(i)<Dmin,
            aq(i)=Dmin;
        end
    end

    % Plot the magnitude response of the direct form filter
    figure(1)
    [h,f]=freqz(bq,aq,nfft,Fs);
    plot(f/1E3,10*log10(abs(h).^2),'r','LineWidth',2);
    grid on
    axis([0 22.050 -60 10]);
    xlabel('Frequency (kHz)');
    ylabel('|H(e^{j\omega})|^2 (dB)')
    title('Magnitude response');
    legend('ideal','direct form');

    % Plot the poles and zeros of the direct form filter
    figure(2)
    p=roots(aq);
    z=roots(bq);
    plot(real(p),imag(p),'rx',real(z),imag(z),'ro',...
        cos(2*pi*[0:0.01:1]),sin(2*pi*[0:0.01:1]),'--');
    axis([-1.75 1.75 -1.75 1.75]);
    axis square
    grid on
    xlabel('Real part');
    ylabel('Imaginary part');
    title('Poles and zeros');
    legend('ideal poles','ideal zeros','direct form poles','direct form zeros');
      
    pause

% ------------------------------------------------------------------
% QUANTIZED SOS FILTER
% ------------------------------------------------------------------

    % first make unquantized SOS filters
    p=roots(a);                     % poles of unquantized filter
    z=roots(b);                     % zeros of unquantized filter
    A=zeros(length(p)/2,3);         % matrix for SOS feedback
    B=zeros(length(z)/2,3);         % matrix for SOS feedforward
    index=1;
    for i=1:size(A,1),
        A(i,:)=poly([p(index) p(index+1)]);
        B(i,:)=poly([z(index) z(index+1)]);
        index=index+2;
    end
    B=B*b(1)^(1/size(B,1));         % scale numerator coefficients equally

    % now quantize
    Bq=round(B*2^Nfrac_sos)/2^Nfrac_sos;        % quantized numerators
    Aq=round(A*2^Dfrac_sos)/2^Dfrac_sos;        % quantized denominator
    Nmax=sum(2.^[0:Nbits_sos-1])/2^Nfrac_sos;   % upper limit on numerator coefficients
    Nmin=-2^(Nbits_sos-Nfrac_sos);              % lower limit on numberator coefficients
    Dmax=sum(2.^[0:Dbits_sos-1])/2^Dfrac_sos;   % upper limit on denominator coefficients
    Dmin=-2^(Dbits_sos-Dfrac_sos);              % lower limit on denominator coefficients
    for i=1:size(A,1),                          % Now correct for any overflow
        for j=1:size(A,2),
        if Bq(i,j)>Nmax,                        % (assume saturation)
            Bq(i,j)=Nmax;
        elseif Bq(i,j)<Nmin,
            Bq(i,j)=Nmin;
        end
        if Aq(i,j)>Dmax,
            Aq(i,j)=Dmax;
        elseif Aq(i,j)<Dmin,
            Aq(i,j)=Dmin;
        end
        end
    end
    
    % now reassemble actual quantized overall transfer function
    p=[];
    z=[];
    for i=1:size(A,1),
        p=[p ; roots(Aq(i,:))];
        z=[z ; roots(Bq(i,:))];    
    end
    aq=poly(p);
    bq=poly(z)*b(1);

    % Plot the magnitude response of the SOS filter
    figure(1)
    [h,f]=freqz(bq,aq,nfft,Fs);
    plot(f/1E3,10*log10(abs(h).^2),'g--','LineWidth',2);
    grid on
    axis([0 22.050 -60 10]);
    xlabel('Frequency (kHz)');
    ylabel('|H(e^{j\omega})|^2 (dB)')
    title('Magnitude response');
    legend('ideal','direct form','SOS');
    hold off

    % Plot the poles and zeros of the SOS filter
    figure(2)
    p=roots(aq);
    z=roots(bq);
    plot(real(p),imag(p),'gx',real(z),imag(z),'go',cos(2*pi*[0:0.01:1]),sin(2*pi*[0:0.01:1]),'--');
    axis([-1.75 1.75 -1.75 1.75]);
    axis square
    grid on
    xlabel('Real part');
    ylabel('Imaginary part');
    title('Poles and zeros');
    legend('ideal poles','ideal zeros','direct form poles','direct form zeros','', 'SOS poles','SOS zeros');
    hold off