Contents

% prodquant.m demo
%
% The purpose of this demo is to show that products have more precision
% than their arguments and must be quantized to keep the number of bits
% from growing. The effect of quantization within the filter is the same as
% injecting noise at various points along the flowgraph, and therefore
% different filter structures will produce different noise spectra at the
% output.  The injected noise is filtered by different "downstream" filter
% elements depending on the structure.  When evaluating product
% quantization it is important to use quantized coefficients even when
% testing unquantized products, or the different frequency responses will
% contribute to the error along with the product roundoff error.
%
% This demo attempts to demonstrate the differences in roundoff noise
% between an unquantized 1st order IIR filter (no roundoff noise, so
% structure won't matter), a quantized 1st order DF2 IIR filter, and a
% quantized 1st order DF1 IIR filter.  With the assumption of a double
% precision add (available in most DSP hardware), the DF2 structure only
% requires a product quantization on the feedback product, and another at
% the output for the feedforward products.  The DF1 structure only requires
% one product quantization of the sum of all the products before the output
% is fed back into the filter.
%
% The ideal output will be used to compare and compute the error due to
% roundoff noise, and the differences in noise power and spectrum can be
% observed between the two structures.
%
% Note that outputs are not checked for overflow in this example, so signal
% levels and pole radii should be kept fairly low.

% Copyright Anderson and Padgett 2011
% Last modified 6/17/1 WTP - original

clear all
close all

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

[b,a] = butter(1,0.0333);
bits = 16;
fbits= bits-1; % assuming Q1.15, there are 15 fraction bits.
bq = q(b,fbits); % assume a max of 1, round off anything less
                  % than the smallest fractional part and restore magnitude
aq = q(a,fbits);

Calculate filtered outputs

% ideal
y_ideal = filter(bq,aq,x); % must use quantized coefficients so the filter match

% DF2
% Generally, for loops in Matlab code are to be avoided, and this is easily
% done with the fixed-point toolbox, but for this introductory demo, I want
% to be as explicit as possible and so need to avoid encapsulating
% operations in powerful functions.  So, a for loop is necessary here, with
% apologies.
N = length(x);
w = zeros(size(x)); % initialize
y_df2 = zeros(size(x));
w(1) = x(1); % assume w[-1]=0
for i = 2:N
    w(i) = x(i) + q(-aq(2)*w(i-1),fbits); % flip the sign on a for flow graph
    y_df2(i) = q(bq(1)*w(i) + bq(2)*w(i-1),fbits);
end

% DF1
y_df1 = zeros(size(x));
y(1) = q(bq(1)*x(1),fbits);
for i=2:N
    y_df1(i) = q(bq(1)*x(i)+bq(2)*x(i-1)-aq(2)*y_df1(i-1),fbits);
end

Plot the filtered outputs for comparison

figure(1)
nn = 1:N;
plot(nn,y_ideal,nn,y_df2,nn,y_df1)
title('A comparison of the ideal, DF2, and DF1 outputs - they are very similar')
legend('ideal','DF2','DF1')

zrange = 9956:9965;
figure(2)
plot(nn(zrange),y_ideal(zrange),nn(zrange),y_df2(zrange),nn(zrange),y_df1(zrange))
title('A comparison of the ideal, DF2, and DF1 outputs - they are very similar')
legend('ideal','DF2','DF1')
xlabel('Zoomed in x axis shows the tiny differences')

Compare the error spectra for DF2 and DF1

errdf2 = y_df2-y_ideal;
errdf1 = y_df1-y_ideal;
pdf2 = 10*log10(abs(fft(errdf2)).^2);
pdf1 = 10*log10(abs(fft(errdf1)).^2);

figure(3)
plot(nn,pdf2,nn,pdf1)
title('A comparison of the error spectrum of the DF2 and DF1 implementations')
legend('DF2','DF1')
xlabel('Direct calculation with the FFT produces a noisy estimate')

[Pdf2,w] = pwelch(errdf2,hamming(256));
[Pdf1,w] = pwelch(errdf1,hamming(256));

q15pow = (2/2^16)^2/12; % total power from 1 q15 quantizer
pred_df1 = q15pow*abs(1./(1+aq(2)*exp(-j*w))).^2; % power output with quantizer filtered only by poles
pred_df2 = q15pow*abs((bq(1)+bq(2)*exp(-j*w))./(1+aq(2)*exp(-j*w))).^2 + q15pow; % power output with one quantizer at output (unfiltered) and one at the input (filtered by whole filter)
figure(4)
% the predicted power must be scaled down by a factor of pi to match the
% spectral density function produced by pwelch (which should integrate to
% the total power over the radian frequency range 0 to pi).
plot(w,10*log10(Pdf2),w,10*log10(Pdf1),w,10*log10(pred_df2/pi),w,10*log10(pred_df1/pi))
legend('DF2 actual','DF1 actual','DF2 predicted','DF1 predicted')
title('A comparison of the error spectrum of the DF2 and DF1 implementations')
xlabel('Pwelch() gives a better estimate, consistent with the theory')
ylabel('Noise power spectral density in dB')

Listen to the difference in error spectra

h = msgbox('Can you hear the difference?  DF2 noise x 10000','Info');
player = audioplayer(errdf2*10^4,fs,16);
playblocking(player);
close(h)

h = msgbox('Can you hear the difference?  DF1 noise x 10000','Info');
player = audioplayer(errdf1*10^4,fs,16);
playblocking(player);
close(h)