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)