Contents
SECTION: INITIALIZATION
clear all;
close all;
[x,fs] = wavread('the pipe began to leak 2.wav');
disp('filter coefficients')
[b,a]=ellip(6,1,40,.2)
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);
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)
disp('ideal sos matrix')
[sosm,g] = tf2sos(b,a)
sections = size(sosm,1);
sosm(1,1:3) = sosm(1,1:3)*g^(1/sections);
sosm(2,1:3) = sosm(2,1:3)*g^(1/sections);
sosm(3,1:3) = sosm(3,1:3)*g^(1/sections);
if wdl < 2, error('This code assumes a wordlength greater than 2'); end;
ibits = 2
fbits = wdl-ibits
qb = 2^-fbits*round(b*2^fbits)
qa = 2^-fbits*round(a*2^fbits)
qsosm = 2^-fbits*round(sosm*2^fbits)
y2 = filter(qb,qa,x);
w3 = filter(qsosm(1,1:3),qsosm(1,4:6),x);
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
y2i = y2;
else
y2i = 0;
end
if max(abs(y3))<Inf
y3i = y3;
else
y3i = 0;
end
sc = max(abs([y1; y2i; y3i]));
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';
else
dstr = 'b';
end
if sunstable
sstr = 'r';
else
sstr = 'b';
end
scrsz = get(0,'ScreenSize');
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');
subplot(312)
plot(y2,dstr)
hlabel = ylabel('direct');
set(hlabel,'Position',pos);
subplot(313)
plot(y3,sstr)
hlabel = ylabel('cascade');
set(hlabel,'Position',pos);
h1 = figure(2);
set(h1,'OuterPosition',[scrsz(3)/3 1 scrsz(3)/3 scrsz(4)])
drawnow();
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)
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);
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);
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);
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
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