Contents

% coefquant2.m
% coefficient quantization

% Copyright 2011
% Wayne T. Padgett
% Last Modified - 4/16/2011 WTP
% 6/17/11 Added blocking audio player and msgbox() information for prettier
% user interface WTP
% 7/9/11 Made it publish more nicely WTP

%
% Most of this code is user interface, the actual signals part is pretty
% short
%

SECTION: INITIALIZATION

clear all;  % cleans up first
close all;

[x,fs] = wavread('the pipe began to leak 2.wav');

disp('filter coefficients')
[b,a]=ellip(6,1,40,.2) % create an example filter

wdl_list = [16, 11, 10, 7, 6];

h = msgbox('Playing Original Unfiltered','Info');
player = audioplayer(x,fs,16);
playblocking(player);
close(h)

y1 = filter(b,a,x);  % apply the unquanitzed coefficients
h = msgbox('Playing Unquantized Low Pass Filtered (all outputs should sound like this)','Info','replace');
player = audioplayer(y1,fs,16);
playblocking(player);
close(h)

for j = 1:length(wdl_list)

new wordlength

    wdl = wdl_list(j)  % wordlength
    % both quanitzed systems work at 11 bits or more
    % direct form fails for 10 bits
    % cascaded sos form fails for 6 bits

    disp('ideal sos matrix')
    [sosm,g] = tf2sos(b,a)     % convert it to second order sections
    sections = size(sosm,1);
    sosm(1,1:3) = sosm(1,1:3)*g^(1/sections); % split the gain across the sections
    sosm(2,1:3) = sosm(2,1:3)*g^(1/sections); % not ideal but good enough for now
    sosm(3,1:3) = sosm(3,1:3)*g^(1/sections);

    % for simplicity (not realism), let's assume a max value of 2 (Q2.x) for
    % all values
    if wdl < 2, error('This code assumes a wordlength greater than 2'); end;
    ibits = 2 % assume 2 integer bits for a range of +/- 2
    fbits = wdl-ibits % all the rest of the bits are fraction bits
    qb = 2^-fbits*round(b*2^fbits) % quantize the direct form numerator
    qa = 2^-fbits*round(a*2^fbits) % quantize the direct form denominator
    qsosm = 2^-fbits*round(sosm*2^fbits) % quantize the second order section matrix

    y2 = filter(qb,qa,x); % apply the quantized direct form coefficients

    w3 = filter(qsosm(1,1:3),qsosm(1,4:6),x);  % apply the second order section coefifcients in stages
    ww3 = filter(qsosm(2,1:3),qsosm(2,4:6),w3);
    y3 = filter(qsosm(3,1:3),qsosm(3,4:6),ww3);

    if max(abs(y2))<Inf  % if the output is unstable, don't use it for scaling
        y2i = y2;
    else
        y2i = 0;
    end
    if max(abs(y3))<Inf
        y3i = y3;
    else
        y3i = 0;
    end

    sc = max(abs([y1; y2i; y3i])); % scale by the largest (non-infinite) output value

    % wait to play the sounds until after the corresponding pictures are
    % displayed


    % check for stability and mark unstable cases
        % multiply the polynomials back out to get the zplane plot (conv = poly
        % multiplication)
    sb = conv(qsosm(1,1:3),qsosm(2,1:3));
    sb = conv(sb,qsosm(3,1:3));
    sa = conv(qsosm(1,4:6),qsosm(2,4:6));
    sa = conv(sa,qsosm(3,4:6));

    dunstable = max(abs(roots(qa)))>= 1;
    sunstable = max(abs(roots(sa)))>= 1;

    if dunstable
        dstr = 'r';  % red plot for unstable
    else
        dstr = 'b';  % blue plot for stable
    end
    if sunstable
        sstr = 'r';
    else
        sstr = 'b';
    end

    scrsz = get(0,'ScreenSize'); % left, bottom, width, height
    h1 = figure(1);
    set(h1,'OuterPosition',[1 1 scrsz(3)/3 scrsz(4)])

    subplot(311)
    plot(y1)
    title(['Filter Output, Wordlength = ' num2str(wdl_list(j))])
    hlabel = ylabel('unquantized');
    pos = get(hlabel,'Position'); % get a position to use to workaround label position bug

    subplot(312)
    plot(y2,dstr)
    hlabel = ylabel('direct');  % ylabel seems to fail when axis range is huge (unstable signals)
    set(hlabel,'Position',pos);
    subplot(313)
    plot(y3,sstr)
    hlabel = ylabel('cascade');
    set(hlabel,'Position',pos);

    h1 = figure(2);   % compare the pole-zero plots
    set(h1,'OuterPosition',[scrsz(3)/3 1 scrsz(3)/3 scrsz(4)])
    drawnow();  % strangely, the position doesn't get set in time and the plots come out different sizes without this
    subplot(311)
    zplane(b,a)
    title(['Pole-Zero Plot, Wordlength = ' num2str(wdl_list(j))])
    ylabel('unquantized')
    subplot(312)
    zplane(qb,qa,dstr)
    ylabel('direct')
    subplot(313)
    % multiply the polynomials back out to get the zplane plot (conv = poly
    % multiplication)
    zplane(sb,sa,sstr)
    ylabel('cascade')

    h1 = figure(3);
    set(h1,'OuterPosition',[2*scrsz(3)/3 1 scrsz(3)/3 scrsz(4)])
    [hh1,ww] = freqz(b,a);  % need to compare all three
    hh2 = freqz(qb,qa);
    hh3 = freqz(sb,sa);
    subplot(311)
    plot(ww,20*log10(abs(hh1)))
    title(['H(e^{j\omega}) ,Wordlength = ' num2str(wdl_list(j))])
    ylabel('unquantized')
    subplot(312)
    plot(ww,20*log10(abs(hh2)),dstr)
    ylabel('direct')
    subplot(313)
    plot(ww,20*log10(abs(hh3)),sstr)
    ylabel('cascade')
    drawnow();

    h = msgbox('Playing Direct Form','Info','replace');
    player = audioplayer(y2/sc,fs,16);
    playblocking(player);
    close(h)

    h = msgbox('Playing Cascade','Info','replace');
    player = audioplayer(y3/sc,fs,16);
    playblocking(player);
    close(h); % get rid of old message box
    options.Default = 'Continue';
    options.Interpreter = 'tex';
    switch(wdl_list(j))
        case 16
            qtxt = {'16 bits is enough to represent the coefficients adequately in either form, so all three filters perform correctly.'  ' ' 'View another wordlength or quit?'};
            disp(textwrap(cellstr(qtxt),80))
            button = questdlg(qtxt,'Continue?','Continue','Quit',options);
            if strcmp(button, 'Quit'), break; end;

        case 11
            qtxt = {'16 down to 11 bits is enough for stability in both forms but direct form is more sensitive to quantization error, so 11 bits is NOT enough to represent the coefficients adequately in direct form. It IS enough for cascade form because the coefficient values stay in a smaller range and do not interact between sections.  You should be able to see that the pole zero plot does not match the unquantized one and that the frequency response does not match the unquantized one.  Note that both quantized filters are stable but the direct form does not produce the desired response due to coefficient quantization.'  ' ' 'You may wish to turn your volume down a little since the unstable output (coming next) is unpleasant to listen to.' ' ' 'View another wordlength or quit?'};
            disp(textwrap(cellstr(qtxt),80))
            button = questdlg(qtxt,'Continue?','Continue','Quit',options);
            if strcmp(button, 'Quit'), break; end;

        case 10
            qtxt = {'10 bits is NOT enough to preserve stability in direct form but it is enough for cascade form.  You should see (and hear) that the direct form output becomes unstable while the others work correctly.  The unstable waveform is red. You should see poles outside the unit circle in the zplane plot.  The frequency response plot is not really a frequency response since the output wouldn''t be sinusoidal, but we can still evaluate H(e^{j\omega})'  ' ' 'View another wordlength or quit?'};
            disp(textwrap(cellstr(qtxt),80))
            button = questdlg(qtxt,'Continue?','Continue','Quit',options);
            if strcmp(button, 'Quit'), break; end;

        case 7
            qtxt = {'The cascade form works down to 7 bits, and the response is reasonably close to the unquantized case.  The cascade case has become more unstable (poles farther outside the unit circle).  The frequency response plot is not really a frequency response since the output wouldn''t be sinusoidal, but we can still evaluate H(e^{j\omega})'  ' ' 'View another wordlength or quit?'};
            disp(textwrap(cellstr(qtxt),80))
            button = questdlg(qtxt,'Continue?','Continue','Quit',options);
            if strcmp(button, 'Quit'), break; end;

        case 6
            options.Default = 'Done';
            qtxt = {'The cascade form becomes marginally stable (poles on the unit circle) at 6 bits.  The result is that the filter integrates any energy at the pole location - so the output always gets bigger and never dies out, but only grows when the input is active instead of exponentially as the unstable case does.  The frequency response plots are not really frequency responses since the outputs wouldn''t be sinusoidal, but we can still evaluate H(e^{j\omega})'  ' ' 'See the command window for some of the filter coefficient values used in this demo.'};
            disp(textwrap(cellstr(qtxt),80))
            button = questdlg(qtxt,'Continue?','Done',options);  % should be the last case
            if strcmp(button, 'Quit'), break; end;

        case default
            qtxt = {'You have chosen a bit length I didn''t comment on.'  ' ' 'View another wordlength or quit?'};
            disp(textwrap(cellstr(qtxt),80))
            button = questdlg(qtxt,'Continue?','Continue','Quit',options);
            if strcmp(button, 'Quit'), break; end;

    end
wdl =

    16

ideal sos matrix

sosm =

    1.0000   -0.0761    1.0000    1.0000   -1.5758    0.6675
    1.0000   -1.3889    1.0000    1.0000   -1.5831    0.8672
    1.0000   -1.5368    1.0000    1.0000   -1.5959    0.9724


g =

    0.0160


ibits =

     2


fbits =

    14


qb =

    0.0161   -0.0482    0.0859   -0.0989    0.0859   -0.0482    0.0161


qa =

    1.0000   -4.7549   10.0434  -11.9260    8.3647   -3.2804    0.5629


qsosm =

    0.2523   -0.0192    0.2523    1.0000   -1.5758    0.6675
    0.2523   -0.3503    0.2523    1.0000   -1.5831    0.8672
    0.2523   -0.3876    0.2523    1.0000   -1.5959    0.9724

    '16 bits is enough to represent the coefficients adequately in either form, so '
    'all three filters perform correctly.'
    ''
    'View another wordlength or quit?'

wdl =

    11

ideal sos matrix

sosm =

    1.0000   -0.0761    1.0000    1.0000   -1.5758    0.6675
    1.0000   -1.3889    1.0000    1.0000   -1.5831    0.8672
    1.0000   -1.5368    1.0000    1.0000   -1.5959    0.9724


g =

    0.0160


ibits =

     2


fbits =

     9


qb =

    0.0156   -0.0488    0.0859   -0.0996    0.0859   -0.0488    0.0156


qa =

    1.0000   -4.7559   10.0430  -11.9258    8.3652   -3.2813    0.5625


qsosm =

    0.2520   -0.0195    0.2520    1.0000   -1.5762    0.6680
    0.2520   -0.3496    0.2520    1.0000   -1.5840    0.8672
    0.2520   -0.3867    0.2520    1.0000   -1.5957    0.9727

    '16 down to 11 bits is enough for stability in both forms but direct form is '
    'more sensitive to quantization error, so 11 bits is NOT enough to represent the '
    'coefficients adequately in direct form. It IS enough for cascade form because '
    'the coefficient values stay in a smaller range and do not interact between '
    'sections.  You should be able to see that the pole zero plot does not match the '
    'unquantized one and that the frequency response does not match the unquantized '
    'one.  Note that both quantized filters are stable but the direct form does not '
    'produce the desired response due to coefficient quantization.'
    ''
    'You may wish to turn your volume down a little since the unstable output '
    '(coming next) is unpleasant to listen to.'
    ''
    'View another wordlength or quit?'

wdl =

    10

ideal sos matrix

sosm =

    1.0000   -0.0761    1.0000    1.0000   -1.5758    0.6675
    1.0000   -1.3889    1.0000    1.0000   -1.5831    0.8672
    1.0000   -1.5368    1.0000    1.0000   -1.5959    0.9724


g =

    0.0160


ibits =

     2


fbits =

     8


qb =

    0.0156   -0.0469    0.0859   -0.0977    0.0859   -0.0469    0.0156


qa =

    1.0000   -4.7539   10.0430  -11.9258    8.3633   -3.2813    0.5625


qsosm =

    0.2539   -0.0195    0.2539    1.0000   -1.5742    0.6680
    0.2539   -0.3516    0.2539    1.0000   -1.5820    0.8672
    0.2539   -0.3867    0.2539    1.0000   -1.5977    0.9727

    '10 bits is NOT enough to preserve stability in direct form but it is enough for '
    'cascade form.  You should see (and hear) that the direct form output becomes '
    'unstable while the others work correctly.  The unstable waveform is red. You '
    'should see poles outside the unit circle in the zplane plot.  The frequency '
    'response plot is not really a frequency response since the output wouldn't be '
    'sinusoidal, but we can still evaluate H(e^{j\omega})'
    ''
    'View another wordlength or quit?'

wdl =

     7

ideal sos matrix

sosm =

    1.0000   -0.0761    1.0000    1.0000   -1.5758    0.6675
    1.0000   -1.3889    1.0000    1.0000   -1.5831    0.8672
    1.0000   -1.5368    1.0000    1.0000   -1.5959    0.9724


g =

    0.0160


ibits =

     2


fbits =

     5


qb =

    0.0313   -0.0625    0.0938   -0.0938    0.0938   -0.0625    0.0313


qa =

    1.0000   -4.7500   10.0313  -11.9375    8.3750   -3.2813    0.5625


qsosm =

    0.2500   -0.0313    0.2500    1.0000   -1.5625    0.6563
    0.2500   -0.3438    0.2500    1.0000   -1.5938    0.8750
    0.2500   -0.3750    0.2500    1.0000   -1.5938    0.9688

    'The cascade form works down to 7 bits, and the response is reasonably close to '
    'the unquantized case.  The cascade case has become more unstable (poles farther '
    'outside the unit circle).  The frequency response plot is not really a '
    'frequency response since the output wouldn't be sinusoidal, but we can still '
    'evaluate H(e^{j\omega})'
    ''
    'View another wordlength or quit?'

wdl =

     6

ideal sos matrix

sosm =

    1.0000   -0.0761    1.0000    1.0000   -1.5758    0.6675
    1.0000   -1.3889    1.0000    1.0000   -1.5831    0.8672
    1.0000   -1.5368    1.0000    1.0000   -1.5959    0.9724


g =

    0.0160


ibits =

     2


fbits =

     4


qb =

         0   -0.0625    0.0625   -0.1250    0.0625   -0.0625         0


qa =

    1.0000   -4.7500   10.0625  -11.9375    8.3750   -3.2500    0.5625


qsosm =

    0.2500         0    0.2500    1.0000   -1.5625    0.6875
    0.2500   -0.3750    0.2500    1.0000   -1.5625    0.8750
    0.2500   -0.3750    0.2500    1.0000   -1.6250    1.0000

    'The cascade form becomes marginally stable (poles on the unit circle) at 6 '
    'bits.  The result is that the filter integrates any energy at the pole location '
    '- so the output always gets bigger and never dies out, but only grows when the '
    'input is active instead of exponentially as the unstable case does.  The '
    'frequency response plots are not really frequency responses since the outputs '
    'wouldn't be sinusoidal, but we can still evaluate H(e^{j\omega})'
    ''
    'See the command window for some of the filter coefficient values used in this '
    'demo.'

end % for j = 1:length(wdl_list)
filter coefficients

b =

    0.0160   -0.0482    0.0860   -0.0989    0.0860   -0.0482    0.0160


a =

    1.0000   -4.7549   10.0434  -11.9260    8.3648   -3.2804    0.5629