FIR filters, IIR filters,
and the linear constant-coefficient difference equation
Causal Moving Average (FIR) Filters
We've discussed systems in which each sample of the output is a weighted sum of (certain of the) the samples of the input.
Let's take a causal "weighted sum" system, where causal means that a given output sample depends only on the current input sample and other inputs earlier in the sequence. Neither linear systems in general, nor finite impulse response systems in particular, need to be causal. However, causality is convenient for a kind of analysis that we're going to explore soon.
If we symbolize the inputs as values of a vector x, and the outputs as corresponding values of a vector y, then such a system can be written as
y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(Nb+1)*x(n-Nb)
where the b values are "weights" applied to the current and earlier input samples to get the current output sample. We can think of the expression as an equation, with the equals sign meaning "equals", or as a procedural instruction, with the equals sign meaning assignment.
Let's write the expression for each output sample as a MATLAB loop of assignment statements, where x is an N-length vector of input samples, and b is an M-length vector of weights. In order to deal with the special case at the start, we will embed x in a longer vector x_hat whose first M-1 samples are zero.
We will write the weighted summation for each y(n) as an inner product, and will do some manipulations of the the inputs (like reversing b) to this end.
This kind of system is often called a "moving average" filter, for obvious reasons.
function y = mafilt(b,x) % causal Moving Average filter % y = mafilt(b,x) % b is system impulse respose; x is input signal % Make sure x is a column vector: [r c] = size(x); if(c > r) x = x'; end N = length(x); % Make sure b is a row vector -- then reverse it [r c] = size(b); if(r > c) b = b'; end M = length(b); brev = b(M:-1:1); xhat = zeros(N+M-1,1); xhat(M:(N+M-1)) = x; y = zeros(N,1); for n=1:N, y(n) = brev * xhat(n:(n+M-1)); end return
From our earlier discussions, it should be obvious that such a system is linear and shift-invariant. Of course, it would be much faster to use the MATLAB convolution function conv() instead of our mafilt().
>> x = zeros(8,1); x(1) = 1;
>> y1 = conv([1 2 3 4 5],x); y2 = mafilt([1 2 3 4 5],x);
>> [ y1(1:8) y2 ]
ans =
1 1
2 2
3 3
4 4
5 5
0 0
0 0
0 0
Instead of considering the first M-1 samples of the input to be zero, we could consider them to be the same as the last M-1 samples. This is the same as treating the input as periodic. We'll use cmafilt() as the name of the function, a small modification of the earlier mafilt() function. In determining the impulse response of a system, there is usually no difference between these two, since all non-initial samples of the input are zero:
>> y3 = cmafilt([1 2 3 4 5],x);
>> [y1(1:8) y2 y3]
ans =
1 1 1
2 2 2
3 3 3
4 4 4
5 5 5
0 0 0
0 0 0
0 0 0
Since a system of this kind is linear and shift-invariant, we know that its effect on any sinusoid will be only to scale and shift it. Here it matters that we use the circular version!
>> s1 = makecos(16,1,1,pi/2)';
>> y1 = mafilt([.1 .2 .3 .4 .5],s1);
>> y2 = cmafilt([.1 .2 .3 .4 .5],s1);
>> [y1(1:8) y2(1:8)]
ans =
0.0000 -1.1582
0.0383 -0.8213
0.1472 -0.3594
0.3486 0.1573
0.6500 0.6500
1.0437 1.0437
1.2786 1.2786
1.3188 1.3188
>> plot(s1,':');
>> hold on
>> plot(y1,'x'); plot(y2,'o');
>> title('Input: (...), circular conv.: (o), non-circular conv: (x)')
The circularly-convolved version is shifted and scaled a bit, while the version with 'ordinary' convolution is distorted at the start...
Let's see what the exact scaling and shifting is by using an fft:
>> S1 = fft(s1); Y2 = fft(y2);
>> [ abs(fftshift(S1)) abs(fftshift(Y2)) ]
ans =
0.0000 0
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
8.0000 10.6251
0.0000 0.0000
8.0000 10.6251
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
Both input and output have amplitude only at frequencies 1 and -1, which is as it should be, given that the input was a sinusoid and the system was linear. The output values are greater by a ratio of 10.6251/8 = 1.3281. This is the gain of the system.
What about the phase? We only need to look where the amplitude is non-zero:
>> [ [angle(S1(2)) angle(S1(16))]' [angle(Y2(2)) angle(Y2(16))]' ]
ans =
-1.5708 -2.6302
1.5708 2.6302
The input has a phase of pi/2, as we requested. The output phase is shifted by an additional 1.0594 (with opposite sign for the negative frequency), or about 1/6 of a cycle to the right, as we can see on the graph.
Now let's try a sinusoid with the same frequency (1), but instead of amplitude 1 and phase pi/2, let's try amplitude 1.5 and phase 0.
>> s3 = makecos(16, 1, 1.5, 0)'; >> y3 = cmafilt([.1 .2 .3 .4 .5],s3); >> S3 = fft(s3); >> Y3 = fft(y3);
We know that only frequency 1 and -1 will have non-zero amplitude, so let's just look at them:
>> [ [abs(S3(2)) abs(S3(16))]', [abs(Y3(2)) abs(Y3(16))]' ] ans = 12.0000 15.9377 12.0000 15.9377
Again the amplitude ratio (15.9377/12.0000) is 1.3281 -- and as for the phase
>> [ [angle(S3(2)) angle(S3(16))]', [angle(Y3(2)) angle(Y3(16))]' ]
ans =
-0.0000 -1.0594
0.0000 1.0594
it is again shifted by 1.0594!
If these examples are typical, we can predict the effect of our system (impulse response [.1 .2 .3 .4 .5]) on any sinusoid with frequency 1 -- the amplitude will be increased by a factor of 1.3281 and the (positive frequency) phase will be shifted by 1.0594.
We could go on to compute the effect of this system on sinusoids of other frequencies by the same methods. But there is a much simpler way, and one which establishes the general point. Since (circular) convolution in the time domain means multiplication in the frequency domain, from
y(n) = x(n) # h(n) % # means convolution here
it follows that
Y(k) = X(k) .* H(k)
or
H(k) = Y(k) ./ X(k)
In other words, the DFT of the impulse response is the ratio of the DFT of the output to the DFT of the input.
In this relationship
H(k) = Y(k) ./ X(k)
the DFT coefficients are complex numbers. Since abs(c1/c2) = abs(c1)/abs(c2) for all complex numbers c1, c2, this equation tells us that the amplitude spectrum of the impulse response will always be the ratio of the amplitude spectrum of the output to that of the input.
In the case of the phase spectrum, angle(c1/c2) = angle(c1) - angle(c2) for all c1, c2 (with the proviso that phases differing by n*2*pi are considered equal). Therefore the phase spectrum of the impulse response will always be the difference between the phase spectra of the output and the input (with whatever corrections by 2*pi are needed to keep the result between -pi and pi).
>> h = zeros(16,1); h(1:5) = [.1 .2 .3 .4 .5];
>> H = fft(h); printfft(H);
DFT =
Real Imaginary Amplitude Phase
-8 0.30000 0 0.30000 0
-7 -0.02572 -0.26604 0.26728 -1.66716
-6 -0.25858 0.12426 0.28689 2.69361
-5 0.18088 0.31957 0.36721 1.05574
-4 0.30000 -0.20000 0.36056 -0.58800
-3 -0.40515 -0.25617 0.47934 -2.57778
-2 -0.54142 0.72426 0.90427 2.21273
-1 0.64998 1.15822 1.32814 1.05940
0 1.50000 0 1.50000 0
1 0.64998 -1.15822 1.32814 -1.05940
2 -0.54142 -0.72426 0.90427 -2.21273
3 -0.40515 0.25617 0.47934 2.57778
4 0.30000 0.20000 0.36056 0.58800
5 0.18088 -0.31957 0.36721 -1.05574
6 -0.25858 -0.12426 0.28689 -2.69361
7 -0.02572 0.26604 0.26728 1.66716
>>
>> H = fftshift(H); f = -8:7;
>> hold off; plot(f,abs(H),'x'); hold on; plot(f,abs(H),':');
>> title('Amp (x) spectrum for h = [.1 .2 .3 .4 .5]');
>> xlabel('Frequency');
>> hold off; plot(f,angle(H),'o'); hold on; plot(f,angle(H),'-');
>> title('Phase (o) spectrum for h = [.1 .2 .3 .4 .5]');
>> xlabel('Frequency');
We can see the phase effects more clearly if we unwrap the representation of phase, i.e. if we add various multiples of 2*pi as needed to minimize the "jumps" that are produced by the periodic nature of the angle() function.
>> hold off; plot(f,unwrap(angle(H)),'o'); hold on; plot(f,unwrap(angle(H)),':');
>> title('Unwrapped phase spectrum for h = [.1 .2 .3 .4 .5]');
>> xlabel('Frequency');
Although the amplitude and phase are usually used for graphical and even tabular presentation, since they are an intuitive way to think about the effects of a system on the various frequency components of its input, the complex Fourier coefficients are more useful algebraically, since they permit the simple expression of the relationship
H(k) = Y(K) ./ X(k)
The general approach we have just seen will work with arbitrary filters of the type sketched, in which each output sample is a weighted sum of some set of input samples.
As mentioned earlier, these are often called Finite Impulse Response filters, because the impulse response is of finite size, or sometimes Moving Average filters.
We can determine the frequency response characteristics of such a filter from the FFT of its impulse response, and we can also design new filters with desired characteristics by IFFT from a specification of the frequency response.
Autoregressive (IIR) Filters
There would be little point in having names for FIR filters unless there were some other kind(s) to distinguish them from, and so those who have studied pragmatics will not be surprised to learn that there is indeed another major kind of linear time-invariant filter.
These filters are sometimes called "recursive" because the value of previous outputs (as well as previous inputs) matters, although the algorithms are generally written using iterative constructs. They are also called Infinite Impulse Response (IIR) filters, because in general their response to an impulse goes on forever. They are also sometimes called "autoregressive" filters, because the coefficients can be thought of as the result of doing linear regression to express signal values as a function of earlier signal values.
The relationship of FIR and IIR filters can be seen clearly in a "linear constant-coefficient difference equation", i.e.
a(1)*y(n) + a(2)*y(n-1) + ... + a(Na+1)*y(n-Na) = b(1)*x(n) + b(2)*x(n-1) + ... + b(Nb+1)*x(n-Nb)
setting a weighted sum of outputs equal to a weighted sum of inputs. This is like the equation that we gave earlier for the causal FIR filter, except that in addition to the weighted sum of inputs, we also have a weighted sum of outputs.
If we want to think of this as a procedure for generating output samples, we need to rearrange the equation to get an expression for the current output sample y(n),
y(n) = (1/a(1)) * ( 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) )
Adopting the convention that a(1) = 1 (e.g. by scaling other a's and b's), we can get rid of the 1/a(1) term:
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)If all the a(n) other than a(1) are zero, this reduces to our old friend the causal FIR filter.
This is the general case of a (causal) LTI filter, and is implemented by the MATLAB function "filter."
>> help filter
FILTER One-dimensional digital filter.
Y = FILTER(B,A,X) filters the data in vector X with the
filter described by vectors A and B to create the filtered
data Y. The filter is a "Direct Form II Transposed"
implementation of the standard difference equation:
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)
If a(1) is not equal to 1, FILTER normalizes the filter
coefficients by a(1).
When X is a matrix, FILTER operates on the columns of X. When X
is an N-D array, FILTER operates along the first non-singleton
dimension.
[Y,Zf] = FILTER(B,A,X,Zi) gives access to initial and final
conditions, Zi and Zf, of the delays. Zi is a vector of length
MAX(LENGTH(A),LENGTH(B))-1 or an array of such vectors, one for
each column of X.
FILTER(B,A,X,[],DIM) or FILTER(B,A,X,Zi,DIM) operates along the
dimension DIM.
See also FILTER2, FILTFILT (in the Signal Processing Toolbox).
Let's look at the case where the b coefficients other than b(1) are zero (instead of the FIR case, where the a(n) are zero):
y(n) = b(1)*x(n) - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
In this case, the current output sample y(n) is calculated as a weighted combination of the current input sample x(n) and the previous output samples y(n-1), y(n-2), etc. To get an idea of what happens with such filters, let's start with the case where:
na is 1 (a "first-order" filter, looking only the immediately previous output sample), b(1) is 1, and a(2) is (say) .5.
That is, the current output sample is the sum of the current input sample and half the previous output sample.
We'll take an input impulse through a few time steps, one at a time.
Step 1
MATLAB n: ? 1 2 3 4 5 6 7 8
__________________________________________________________________
Input: 1 0 0 0 0 0 0 0
Output: 0 0 0 0 0 0 0 0 0
__________________________________________________________________
|
Result: 1
Step 2
MATLAB n: ? 1 2 3 4 5 6 7 8
_________________________________________________________________
Input: 1 0 0 0 0 0 0 0
Output: 0 1 0 0 0 0 0 0 0
_________________________________________________________________
|
Result: .5
Step 3
MATLAB n: ? 1 2 3 4 5 6 7 8
________________________________________________________________
Input: 1 0 0 0 0 0 0 0
Output: 0 1 .5 0 0 0 0 0 0
________________________________________________________________
|
Result: .25
Step 4
MATLAB n: ? 1 2 3 4 5 6 7 8
________________________________________________________________
Input: 1 0 0 0 0 0 0 0
Output: 0 1 .5 .25 0 0 0 0 0
|
Result: .125
It should be clear at this point that we can easily write an expression for the nth output sample value: it is just
.5^(n-1)
(If MATLAB counted from 0, this would be simply .5^n).
Since what we are calculating is the impulse response of the system, we have demonstrated by example that the impulse response can indeed have infinitely many non-zero samples.
To implement this trivial first-order filter in MATLAB, we could use "filter". The call will look like this:
b(n) -- input coeff. a(n) -- output coeff. input signal
filter([1 0], [1 -.5], [1 0 0 0 0 0 0 0])
and the result is:
>> filter([1 0], [1 -.5], [1 0 0 0 0 0 0 0])
ans =
Columns 1 through 7
1.0000 0.5000 0.2500 0.1250 0.0625 0.0312 0.0156
Column 8
0.0078
Is this business really still linear?
We can look at this empirically:
r = rand(1,6) s = rand(1,6) filter([1 0],[1 -.5],r+s) filter([1 0],[1 -.5],r) + filter([1 0],[1 -.5],s)
For a more general approach, consider the value of an output sample y(n).
It is x(n)+.5*y(n-1),
and y(n-1) in turn is
x(n-1)+.5*y(n-2),
and y(n-2) in turn is
x(n-2)+.5*y(n-3), and so on.
By successive substitution we could write this as
y(n)=x(n)+.5*(
x(n-1)+.5*(
x(n-2)+.5*(
x(n-3)+.5*( ...
or again as
y(n)=x(n)+.5*x(n-1)+.25*x(n-2)+.125*x(n-3)+...
This is just
y(n) = .5^0 * x(n) + .5^1 * x(n-1) + .5^2 * x(n-2) + .5^3 * x(n-3) + ...
or
This is just like our old friend the convolution-sum form of an FIR filter, with the impulse response provided by the expression .5^k, and the length of the impulse response being infinite. Thus the same arguments that we used to show that FIR filters were linear will now apply here.
So far this may seem like a lot of fuss about not much. What is this whole line of investigation good for?
We'll answer this question in stages, starting with an example.
It's not a big surprise that we can calculate a sampled exponential by recursive multiplication. Let's look at a recursive filter that does something less obvious. This time we'll make it a second-order filter, so that the call to "filter" will be of the form
input coeff. output coeff. input vector
X = filter( [1 0 0], [1 a2 a3], x )
Let's set the second output coefficient a2 to -2*cos(2*pi/40), and the third output coefficient a3 to 1, and look at the impulse response.
>> x = zeros(1,64); x(1) = 1; >> X = filter([1 0 0], [1 -2*cos(2*pi/40) 1], x); >> plot(X,'o');
Not very useful as a filter, actually, but it does generate a sampled sine wave (from an impulse) with three multiply-adds per sample! In order to understand how and why it does this, and how recursive filters can be designed and analyzed in the more general case, we need to step back and take a look at some other properties of complex numbers, on the way to understanding the z transform.