Bruce Land
One of the most demanding applications for fast arithmetic is digital flitering. Atmel application note AVR201 shows how to use the hardware multiplier to make a multiply-and-accumulate operation (MAC). The MAC is the basis for computing the sum of terms for any digital filter, such as the "Direct Form II Transposed" implementation of a difference equation shown below (from the Matlab filter description).
a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
- a(2)*y(n-1) - ... - a(na+1)*y(n-na) The y(n) term is the current filter output and x(n) is the current filter input. The a's and b's are constant coefficients, determined by some design process, such as the Matlab butter routine. Usually, a(1) is set to one, but as we shall see below, setting a(1) to a power of two may increase filter accuracy, at the cost of slightly increased execution time. The y(n-x) and x(n-x)
are past outputs and inputs respectively. You can use the Matlab filter
design tools to determine the a's and b's.
The 2nd order IIR code (plus uart.c, uart.h)
implements an efficient MAC operation using the multfix asssembler
macro and then uses it to produce a second order IIR filter. The
time/sample is 140 cycles. It would not be reasonable to use this
fixed-point algorithm with normalized cutoff frequency below about 0.10 because of the lack of relative precision in the b coefficients. With a normalized cutoff frequency=0.25, the impulse response is accurate to within 1% of peak response. With a normalized cutoff frequency=0.10,
the impulse response is accurate within 5% of peak response. If you
need a cutoff lower than 0.10, consider changing your hardware sample
rate or using CIC filters (see below) to lower the effective sample
rate. It is also possible to scale up all of the filter coefficients by
16 to give 4 bits more percision, then divide the output of the filter
by 16 for each sample. IIRfilterScaledGCC644_macro.c is an example which executes in 190 cycles. This Matlab program allows you to simulate the effect of coefficient truncation on filter response.
Another approach used in speech synthesis is to use impulse-invariant
second order filters which are a bit quicker to execute because they
zero out two of the general second order filter coefficients (b(2) and b(3)).
These filters can also be designed on-the-fly in a C program, rather
than using Matlab. It is possible to build bandpass and lowpass filters
using this technique, but you need to be careful of aliasing effects and
the sensitivity of coefficients to 8-bit truncation. This matlab program
allows you to estimate truncation errors and design filters. Note that
the bandpass filters have a gain greater than one in the passband, so
you may need to scale output to avoid overflow of fixed point numbers.
At 8000 Hz sample rate, you can expect to get reasonable bandpass filter
accuracy down to F=250, BW=25 Hz, and lowpass accuracy down to BW=250
Hz.
Impulse invariant second order filter design process:
- Choose the bandpass center frequency,
F. If the filter is to be lowpass, setF=0. - Choose the bandpass bandwidth,
BW. If the filter is to be lowpass, the cutoff frequency isBW/2. - Compute:
C = -exp(-2*pi*BW*T)whereT=1/(sample rate) - Compute:
B = 2 * exp(-pi*BW*T) * cos(2*pi*F*T)whereT=1/(sample rate) - Compute:
A = 1 - B - C - The filter output at time i is then
output(i)=A*input(i) + B*output(i-1) + C*output(i-2) - In some circumstances generating a table of A, B and C at the
frequencies you need will be the fastest way to change between filters.
These two programs suggest how you might do this for lowpass and bandpass filters
Second order and fourth order Butterworth IIR Filters for 8:8
If you use a Butterworth approximation for lowpass (maximally
flat IIR filter) then you can simplify the general IIR filter further to
only three MAC operations. Since b(1)=b(3)=0.5*b(2), the three multiplies on the first line below can be combined into one multiply by sum of inputs,...
Est
Abhilash Patel
Question regarding filter implementation:
I want to implement IIR BPF and examine it with frequencies
I want to use direct form 1, in the conventional way a(1)y(n)=b(1)x(n)+b(2)x(n−1)-a(2)y(n−1)−a(3)y(n−2).
when feeding the filter with frequency according to the BPF center frequency, the output is well.
BUT when feeding with higher frequencies for example results in convergence section in which the signal amplitde is big
here's my code:
t=1:1000;
y=1*sin(2*pi*10/1000*t);
[b,a] = butter(1,[2.5 7.5]/(Fs/2))
b1(1)=-a(2);
b2(1)=-a(3);
% b3(1)=-a(3);
a0(1)=b(1);
a1(1)=b(2);
a2(1)=b(3);
plot(y)
y_oldestt=0;
y_oldest=0;
y_older=0;
dataOut = filter(b,a,y);
% plot(dataOut);
y_oldest=y(1);
y_older=y(2);
y_hist=[];
inn_idx=1;
for idx=3:length(y)
temp=y_older;
y_older=y_older*b1(inn_idx)+y_oldest*b2(inn_idx)+y(idx)*a0(inn_idx)+y(idx-1)*a1(inn_idx)+y(idx-2)*a2(inn_idx);
y_oldest=temp;
y_hist=[y_hist y_older];
end
figure(1);
subplot(3,1,1)
plot(y);
ylim([-1.5 1.5]);
subplot(3,1,2)
plot(dataOut);
ylim([-1.5 1.5]);
subplot(3,1,3)
plot(y_hist)