Question

In an FIR digital filter, each sample in the output signal is found by multiplying M samples from the input signal by M predetermined coefficients, and adding the products. The characteristics of these filters (high-pass, low-pass, etc.) are determined by the coefficients used. For this problem, assume M = 5000, and that single precision floating point math is used. a. How many math operations (the number of multiplications plus the number of additions) need to be conducted to calculate each point in the output signal? b. If the output signal has an average amplitude of about one-hundred, what is the expected error on an individual output sample? Assume that the round-off errors combine by addition. Give your answer both as an absolute number, and as a fraction of the number being represented. c. Repeat (b) for the case that the round-off errors combine randomly. Explain in detail?

Answer 1

This section provides an example of how you can create a filter with integer coefficients. In this example, a raised-cosine filter with floating-point coefficients is created, and the filter coefficients are then converted to integers.
Define the Filter Coefficients
To illustrate the concepts of using integers with fixed-point filters, this example will use a raised-cosine filter:
b = rcosdesign(.25, 12.5, 8, 'sqrt');
The coefficients of b are normalized so that the passband gain is equal to 1, and are all smaller than 1. In order to make them integers, they will need to be scaled. If you wanted to scale them to use 18 bits for each coefficient, the range of possible values for the coefficients becomes:
[?2
?17
,2
17
?1]==[?131072,131071]
Because the largest coefficient of b is positive, it will need to be scaled as close as possible to 131071 (without overflowing) in order to minimize quantization error. You can determine the exponent of the scale factor by executing:
B = 18; % Number of bits
L = floor(log2((2^(B-1)-1)/max(b))); % Round towards zero to avoid overflow
bsc = b*2^L;
Alternatively, you can use the fixed-point numbers autoscaling tool as follows:
bq = fi(b, true, B); % signed = true, B = 18 bits
L = bq.FractionLength;
It is a coincidence that B and L are both 18 in this case, because of the value of the largest coefficient of b. If, for example, the maximum value of b were 0.124, L would be 20 while B (the number of bits) would remain 18.
Build the FIR Filter
First create the filter using the direct form, tapped delay line structure:
h = dfilt.dffir(bsc);
In order to set the required parameters, the arithmetic must be set to fixed-point:
h.Arithmetic = 'fixed';
h.CoeffWordLength = 18;
You can check that the coefficients of h are all integers:
all(h.Numerator == round(h.Numerator))
ans =
1
Now you can examine the magnitude response of the filter using fvtool:
fvtool(h, 'Color', 'white')

This shows a large gain of 117 dB in the passband, which is due to the large values of the coefficients— this will cause the output of the filter to be much larger than the input. A method of addressing this will be discussed in the following sections.
Set the Filter Parameters to Work with Integers
You will need to set the input parameters of your filter to appropriate values for working with integers. For example, if the input to the filter is from a A/D converter with 12 bit resolution, you should set the input as follows:
h.InputWordLength = 12;
h.InputFracLength = 0;
The info method returns a summary of the filter settings.
info(h)
Discrete-Time FIR Filter (real)   
-------------------------------   
Filter Structure : Direct-Form FIR   
Filter Length : 101
Stable : Yes
Linear Phase : Yes (Type 1)
Arithmetic : fixed
Numerator : s18,0 -> [-131072 131072)
Input : s12,0 -> [-2048 2048)   
Filter Internals : Full Precision
Output : s31,0 -> [-1073741824 1073741824) (auto determined)
Product : s29,0 -> [-268435456 268435456) (auto determined)  
Accumulator : s31,0 -> [-1073741824 1073741824) (auto determined)
Round Mode : No rounding   
Overflow Mode : No overflow
In this case, all the fractional lengths are now set to zero, meaning that the filter h is set up to handle integers.
Create a Test Signal for the Filter
You can generate an input signal for the filter by quantizing to 12 bits using the autoscaling feature, or you can follow the same procedure that was used for the coefficients, discussed previously. In this example, create a signal with two sinusoids:
n = 0:999;
f1 = 0.1*pi; % Normalized frequency of first sinusoid
f2 = 0.8*pi; % Normalized frequency of second sinusoid
x = 0.9*sin(0.1*pi*n) + 0.9*sin(0.8*pi*n);
xq = fi(x, true, 12); % signed = true, B = 12
xsc = fi(xq.int, true, 12, 0);
Filter the Test Signal
To filter the input s
you will notice that the number of bits in the output is considerably greater than in the input. Because such growth in the number of bits representing the data may not be desirable, you may need to truncate the wordlength of the output. The best way to do this is to discard the least significant bits, in order to minimize error. However, if you know there are unused high order bits, you should discard those bits as well.
To determine if there are unused most significant bits (MSBs), you can look at where the growth in WordLength arises in the computation. In this case, the bit growth occurs to accommodate the results of adding products of the input (12 bits) and the coefficients (18 bits). Each of these products is 29 bits long (you can verify this using info(h)). The bit growth due to the accumulation of the product depends on the filter length and the coefficient values- however, this is a worst-case determination in the sense that no assumption on the input signal is made besides, and as a result there may be unused MSBs. You will have to be careful though, as MSBs that are deemed unused incorrectly will cause overflows.
Suppose you want to keep 16 bits for the output. In this case, there is no bit-growth due to the additions, so the output bit setting will be 16 for the wordlength and –14 for the fraction length.
Since the filtering has already been done, you can discard some bits from ysc:
yout = fi(ysc, true, 16, -14);
Alternatively, you can set the filter output bit lengths directly (this is useful if you plan on filtering many signals):
specifyall(h);
h.OutputWordLength = 16;
h.OutputFracLength = -14;
yout2 = filter(h, xsc);
You can verify that the results are the same either way:
norm(double(yout) - double(yout2))
ans =
0
However, if you compare this to the full precision output, you will notice that there is rounding error due to the discarded bits:
norm(double(yout)-double(ysc))
ans =
1.446323386867543e+005
In this case the differences are hard to spot when plotting the data, as seen below:
stem(n(yidx), [double(yout(yidx)'), double(ysc(yidx)')]);
axis([850 950 -2.5e8 2.5e8]);
legend('Scaled Input', 'Output');
set(gcf, 'color', 'white');

Scale the Output
Because the filter in this example has such a large gain, the output is at a different scale than the input. This scaling is purely theoretical however, and you can scale the data however you like. In this case, you have 16 bits for the output, but you can attach whatever scaling you choose. It would be natural to reinterpret the output to have a weight of 2^0 (or L = 0) for the LSB. This is equivalent to scaling the output signal down by a factor of 2^(-14). However, there is no computation or rounding error involved. You can do this by executing the following:
yri = fi(yout.int, true, 16, 0);
stem(n(idx)', [xscext, double(yri(yidx)')]);
axis([800 950 -1.5e4 1.5e4]);
legend('input', 'rescaled output');

This plot shows that the output is still larger than the input. If you had done the filtering in double-precision floating-point, this would not be the case— because here more bits are being used for the output than for the input, so the MSBs are weighted differently. You can see this another way by looking at the magnitude response of the scaled filter:
[H,w] = freqz(h);
plot(w/pi, 20*log10(2^(-14)*abs(H)));

This plot shows that the passband gain is still above 0 dB.
To put the input and output on the same scale, the MSBs must be weighted equally. The input MSB has a weight of 2^11, whereas the scaled output MSB has a weight of 2^(29–14) = 2^15. You need to give the output MSB a weight of 2^11 as follows:
yf = fi(zeros(size(yri)), true, 16, 4);
yf.bin = yri.bin;
stem(n(idx)', [xscext, double(yf(yidx)')]);
legend('input', 'rescaled output');

This operation is equivalent to scaling the filter gain down by 2^(-18).
[H,w] = freqz(h);
plot(w/pi, 20*log10(2^(-18)*abs(H)));

The above plot shows a 0 dB gain in the passband, as desired.