BME BME 4 4 52 52 Bio Bio medi medi cal Signal cal Signal Processing
Processing Lecture Lecture 5 5
Digital filtering Digital filtering
Lecture 5 Outline
In this lecture, we’ll study digital filtering methods
General considerations and filter specifications
Filtering in frequency domain
Filtering in time domain
Sum and difference (SD) filter Finite Impulse Response (FIR) filter
Infinite Impulse Response (IIR) filter using MATLAB
General considerations
As mentioned in earlier lectures, filtering is the process of keeping components of the signal with certain desired
frequencies and removing components of the signal with certain undesired frequencies
Very often, we keep the gain of the required frequency components to 1 or close to 1
And the gain of the undesired frequency components will be 0 or close to 0
In general, there are 4 types of filtering: LPF, HPF, BPF,
BSF
Filter specifications
Passband
the range of frequency components that are allowed to pass
Stopband
the range of frequency components that are suppressed
Passband ripple
ripples in the passband
the maximum amount by which attenuation in the passband may deviate from gain (which is normally 1)
Stopband ripple
Ripples in the stopband
The maximum amount by which attenuation in the stopband may deviate from gain (which is normally 0)
Stopband attenuation
the minimum amount by which frequency components in the stopband are attenuated
Transition band
The band between the passband and the stopband
Ideal filter and edge frequencies
Frequency response of ideal filters
The edge frequencies are the end frequencies of passband or stopband
fc= cut-off frequency LPF: passband: 0 f fc Stopband: fc f fs/2
HPF: passband: fc f fs/2 Stopband: 0 f fc
BPF: passband: fc1 f fc2 Stopband: 0 f fc1 and fc2
f fs/2
BSF: passband: 0 f fc1 and fc2 f fs/2
Stopband: fc1 f fc2
Actual LPF
LPF
Passes all low- frequency
components below fp and blocks all higher frequency components above fs
In reality, you can’t design ‘square’ type of filters
So, there needs to be transition betweens the bands
LPF: passband: 0f fp Stopband: fsf fs/2
Transition band: fp<f <fs
LPF example
Low-pass filter (LPF)
Eg.: Consider a combination of 3 sinusoidal signals - 2 Hz, 5 Hz and 11 Hz.
The final output signals after 2 LPF are shown
LPF at fp=3 Hz and fs=4 Hz
LPF at fp=8 Hz and fs=9 Hz
0 200 400 600 800 1000
-3 -2 -1 0 1 2 3
Combined signal
LPF, fp=8 Hz, fs=9 Hz
0 200 400 600 800 1000
-3 -2 -1 0 1 2 3
0 200 400 600 800 1000
-1.5 -1 -0.5 0 0.5 1 1.5
Only 2 Hz signal remains
Only 2 Hz and 5 Hz signals remain
0 200 400 600 800 1000
-1.5 -1 -0.5 0 0.5 1 1.5
0 200 400 600 800 1000
-1.5 -1 -0.5 0 0.5 1 1.5
0 200 400 600 800 1000
-1.5 -1 -0.5 0 0.5 1 1.5
LPF, fp=3 Hz, fs=4 Hz
2 Hz signal 5 Hz signal 11 Hz signal
HPF
Passes all high-frequency
components above fp and blocks all higher frequency components below fs
Eg.: Consider the same combination of 3 sinusoidal signals, 2 Hz, 5 Hz and 11 Hz.
The final output signals after 2 HPF are shown
HPF at fs=3 Hz and fp=4 Hz
HPF at fs=8 Hz and fp=9 Hz
0 200 400 600 800 1000
-3 -2 -1 0 1 2 3
0 200 400 600 800 1000
-3 -2 -1 0 1 2 3
Combined signal
0 200 400 600 800 1000
-1.5 -1 -0.5 0 0.5 1 1.5
Only 5 Hz and 11 Hz signals remain
Only 11 Hz signal remains HPF, fs=3 Hz, fp=4 Hz
HPF, fs=8 Hz, fp= 9 Hz
From this point
onwards, we will use Fs for sampling frequency
BPF
Passes all frequency components between edge passband frequencies,
fp1<freq(allow)<fp2 and blocks all frequencies below and above edge stopband
frequencies, freq(block)<fs1;freq(block)>fs2
Eg.: Consider the same combination of 3 sinusoidal signals, 2 Hz, 5 Hz and 11 Hz.
The final output signal after BPF at fp1=4 Hz, fp2=6 Hz, fs1=3 Hz, fs2=7 Hz is shown
0 200 400 600 800 1000
-3 -2 -1 0 1 2 3
Combined signal
BPF
0 200 400 600 800 1000
-1.5 -1 -0.5 0 0.5 1 1.5
Only 5 Hz signal remains
BSF
Band-stop filter (BSF)
Passes all frequency components
lower and higher than edge passband frequencies, freq(allow)<fp1;freq(allow)>fp2 and blocks all frequencies between fs1<freq(block)<fs2
Eg.: Consider the same combination of 3 sinusoidal signals, 2 Hz, 5 Hz and 11 Hz.
The final output signal after BSF at fp1=4 Hz, fp2=6 Hz, fs1=3 Hz, fs2=7 Hz is shown
0 200 400 600 800 1000
-3 -2 -1 0 1 2 3
Combined signal
BSF
0 200 400 600 800 1000
-1.5 -1 -0.5 0 0.5 1 1.5
5 Hz signal is filtered out, only 2 Hz and 11 Hz signals remain
Direct filtering (in frequency domain)
A simple method of doing this
Obtain the DFT of the signal (from 0 to fs)
Set to zero the values that are not in the required frequency range i.e. apply a RECTANGULAR window
Compute the Inverse Discrete Fourier Transform (IDFT)
Example: f1=8 Hz and f2=25 Hz with N=100, fs=200 Hz
Say, we wish to design a LPF with fp=10 Hz and fs=12 Hz
Compute y=fft(x) in MATLAB
Set the values y(7:95)=0 WHY? This is important!
Compute yf=ifft(y,’symmetric’)
And you get the low pass filtered signal!
What is the main disadvantage of this technique?
High computation and time
Another disadvantage??
0 20 40 60 80 100
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
0 20 40 60 80 100
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
0 20 40 60 80 100
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
+ =
0 20 40 60 80 100
-1.5 -1 -0.5 0 0.5 1 1.5
In MATLAB, you have to force conjugate symmetry, else you will get complex values due round-off errors in doing FFT and IFFT Low pass filtered signal
Direct filtering (in frequency domain) –cont.
How it works:
We can also use a different window say Hanning window, instead of rectangular window to obtain a smoother filtered output
H=hanning(12);
Hf(1:6) =H(7:12);
H(7:95)=0;
Hf(96:100)=H(2:6);
y=fft(x);
yy=y.*Hf;
yf=ifft(yy,'symmetric');
plot(yf);
0 20 40 60 80 100
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Hanning window
Filtering in time domain
Two disadvantages of direct filtering
computational time and complexity is high
there is distortion
To solve, we should filter in time domain
There are many types of time domain filtering methods
We will look at
Simple FIR filters
IIR filters using MATLAB codes
Filtering in time domain (in equation form)
The output from a IIR digital filter in made up of previous inputs and previous outputs as well
where B and A are the filter coefficients and the operation * is convolution
Convolution in time domain is equivalent to multiplying in frequency domain –do you remember that we did some window multiplication for direct filtering
Convolution operation will not be discussed in this course
The output from a FIR digital filter in made up of previous inputs only, so no feedback
N
j M
k
j n y j A k
n x k B n
y
1 1
] [
* ] [ ]
[
* ] [ ]
[
M
k
k n x k B n
y
1
] [
* ] [ ]
[
Simple low pass FIR filter –sum filter
Consider, y[n]=x[n]+x[n-1] for every data in the signal
This filter is also known as the sum filter
This simple addition acts as a LPF!
This can be proved by using z=transform but is not needed for the purpose of this course
For hardware design, the block diagram would look like
Advantages:
You just need one adder and one delay circuit – simple and cost effective (i.e.
cheap)
The filter coefficients are integer values (in this case they are 1), so no round-off errors
It is an FIR filter so it is stable – why are FIR filters stable? No feedback
Disadvantage:
Not very good LPF (see next slide)z-1 +
x(n) x(n-1) y(n)
Sum filter (cont.)
The frequency response of the sum filter: not very good as it is far from the ideal ‘square’ filter
The gain at freq=0 is 1
The stopband frequency (when gain=0) is at rad/sample or at Fs/2 Hz
So, there is no stopband
So how to define the passband frequency and transition bands?
For these cases, we use the 3 dB cut-off as the passband frequency
Figure from S.K.Mitra, DSP 3e
Magnitude is also known as gain
3 dB cut-off frequency
The 3 dB cut-off frequency is defined as the frequency when the gain drops 3 dB from maximum gain of 1, which is assumed to be 0 dB
Gain=1, 20log
10(1)=0 dB
When energy is half, i.e. gain=(1/2)
0.5=0.7071, we have 20log
10(0.7071)=- 3 dB
From the figure
We can see that the 3 dB cut-off frequency (when gain=0.7071) is approximately 0.5 rad/sample or Fs/4 Hz
This is the passband frequency
So, the passband is from 0 to Fs/4 Hz
And transition band is from Fs/4 to Fs/2
Increasing the order of sum filter
For order 1, we had y[n]=x[n]+x[n-1]
Now assume, this is fed to another same filter in cascade connection
Now, we have z[n]=y[n]+y[n-1]
Solving, we have z[n]=y[n]+y[n-1]
=x[n]+x[n-1]+x[n-1]+x[n-2]
=x[n]+2x[n-1]+x[n-2]
z-1 +
y(n) y(n-1) z(n)
z-1 +
x(n) x(n-1) y(n)
z-1 +
x(n) x(n-1) z(n)
z-1 x(n-2)
2
Verify this using x[1]=3, x[2]=2, x[3]=5
Hint:
Compute y[2] and y[3]
Compute z[3]
using single cascaded filter and two filters and compare
You will notice that z[n] will be defined only for n=3
onwards if x[1] is the starting point, i.e. for every order M, you lose M initial data points in filtering
Likewise y[n] is defined only from y[2] onwards.
Increasing the order of sum filter (cont.)
So for order M, we have
The frequency response for M=3 is given =>
The passband is about 0.302 rad or Fs/6
We can see that with increasing order, the passband is becoming smaller without any change in stopband
Also, the curve is becoming closer to the ideal ‘square’
So, we can increase/decrease M depending on our needs
The 3-dB cut-off frequency is given by fp=
cos
-1(2
-1/(2*M)) *2/pi
]0 [ )
( M x n r
r M Cr n
y
where QCr r!(QQ! r)!
Figure from S.K.Mitra, DSP 3e
What would be y[n] for N=3?
Simple high pass FIR filter
Similarly, a HPF can be designed using y[n]=x[n]-x[n-1]
For hardware design, the block diagram would look like
The stopband frequency (when gain=0) is at 0 Hz, i.e. there is no stopband
The gain at freq=Fs/2 Hz or rad is 1
The passband frequency is at 0.5 rad or Fs/4 Hz (using 3 dB cut-off)
Passband width is from Fs/4 to Fs/2 Hz
The orders, N can be increased to obtain a smaller passband width and to obtain a frequency response closer to the ideal ‘square’ filter
The 3-dB cut-off frequency is given by fp= sin-1(2-1/(2*N)) *2/pi
For order N, we have
z-1 +
x(n)
x(n-1)
- y(n)
N
r r N C r x n r n
y ( ) 0 ( 1 ) [ ]
As homework, try y[n] for N=3?
Figure from S.K.Mitra, DSP 3e
This is also known as a difference filter
Simple HPF FIR filter example
HPF can be used to remove mean values and low frequency polynomial trends, i.e. to detrend the data
As an example, we saw the passenger data plot in Lecture 4
The detrending can be simply done by using, y[n]=x[n]-x[n-1], where x[n] is the data
Figures from R.Shiavi, Introduction to applied statistical signal analysis
Simple band pass FIR filter
Similarly, a BPF can be designed using a combination of LPF and HPF
This is known as sum and difference (SD) filter
Different orders, M and N can be chosen to obtain the required frequency response
Where Gaincf is the gain at centre frequency given by
Gaincf fT N
fT M N f
GM, ( ) (2cos ) 2sin /
M N
N s M
fcentre f cos 1 2
Example (As homework, verify these later on your own)
For filter orders of M=28 and N=8 gives
The centre frequency is 40 Hz when fs=256 Hz
Approximate 3 dB bandwidth from 32 to 48 Hz (rounded to the nearest integer)
The gain amplification at 40 Hz is approximately 47.211
Simple band pass FIR filter (cont.)
Let us compute a band pass FIR filter equation:
For orders, LPF, M=4 and HPF, N=1, obtain the band pass FIR equation that expresses z[n] in terms of x[n] and
delays of x[n]
Solution (can you get the answer?)
z[n]=x[n]+3x[n-1]+2x[n-2]-2x[n-3]-3x[n-4]-x[n-5]
LPF HPF
x(n) y(n) z(n)
N
r r NCry n r n
z( ) 0( 1) [ ]
] 0 [
)
( M x n r
r MCr n
y
IIR filter
The problem with FIR filter is that you do not get ‘square’
passband like in ideal filters unless you use a very high order
To solve this problem, we can use IIR filter
The disadvantages of IIR filters are
The are not stable (due to feedback)
The filter coefficients are not normally integer, so can have round-off errors
Hardware design is more complicated
Most importantly, their phase response is not linear
Phase response?
In you refresh your memory, you’ll know that frequency responses have two parts:
magnitude and phase responses
In the diagrams, the magnitude and phase responses of Butterworth and Elliptic IIR filters are shown
You can see that the gains are relatively stable at 1 for the bandpass range of 0.2 to 0.6
But the phase response is not linear, i.e.
not a straight line
In MATLAB, this problem is solved by filtering twice, once forward and once reverse
By doing so, the magnitude response is squared while the phase response
becomes zero
0 0.2 0.4 0.6 0.8 1
-2000 -1000 0 1000
Normalized Frequency ( rad/sample)
Phase (degrees)
0 0.2 0.4 0.6 0.8 1
-400 -200 0 200
Normalized Frequency ( rad/sample)
Magnitude (dB)
Butterworth
Elliptic
0 0.2 0.4 0.6 0.8 1
-500 0 500
Normalized Frequency ( rad/sample)
Phase (degrees)
0 0.2 0.4 0.6 0.8 1
-80 -60 -40 -20 0
Normalized Frequency ( rad/sample)
Magnitude (dB)
IIR filter design
It is not possible to design IIR filters digitally
So, the approach is to design analogue IIR filters, then use a bilinear method to transform the filter to digital
This ‘bilinear’ method will not be discussed in this course as it requires z –transofrm knowledge
So, we will use MATLAB functions directly
There many types of IIR filters
Butterworth
Elliptic (Cauer)
Chebyshev I
Chebyshev II
Bessel
But we will look at two only: Butterworth and Elliptic filters
Because Butterworth filter gives flat magnitude responses in the passband and stopband (i.e. no ripples or very little ripple) while Elliptic filter requires the lowest order among all the IIR filters
Butterworth IIR filter
Remember the IIR filter equation
We have to determine the orders M and N and then determine the coefficients A and B
Very often, we use M=N
N
j M
k
j n y j A k
n x k B n
y
1 1
] [
* ] [ ]
[
* ] [ ]
[
In MATLAB, we can find the required minimum order for our specification using
buttord (Wp, Ws, Rp, Rs)
where Wp is the passband edge frequency, Ws is the stopband edge frequency, Rp is maximum ripple in passband and Rs is the minimum attenuation in
stopband
Rp and Rs will be in dB, while Ws and Wp will be in normalised radian/sample
Low pass Butterworth filter
with different orders Figure from S.K.Mitra, DSP 3e
Butterworth IIR filter (cont.)
After finding the order using buttord function, we have to find out the filter coefficients B and A using
[B,A]=butter (N,Wp)
where N is the filter order chosen earlier
Finally to actually do the filtering, we use filtfilt (B,A,x) where x is the signal to be filtered
Filtfilt function filter the signal twice and uses convolution operations
Example: Design a lowpass filter with Fs=200 Hz
Passband = 0 to 40 Hz
Passband ripple = less than 3 dB
Stopband = 50 Hz to the Fs/2
Stopband attenuation = at least 30 dB
Plot the filter's frequency response – use freqz (B,A) function for this purpose.
Butterworth IIR filter example - solution
Fs=200 Hz
Rp=3 and Rs=30
Wp=40/100
Ws=50/100
Use buttord (Wp, Ws, Rp, Rs)
N=buttord(40/100, 50/100, 3, 30) which gives N=11
Next obtain the coefficients, B and A using [B,A]=butter (11, 40/100) which gives
B = 0.0002 0.0026 0.0129 0.0386 0.0772 0.1081 0.1081 0.0772 0.0386 0.0129 0.0026 0.0002
A =1.0000 -2.1931 3.5467 -3.6414 2.9012 -1.7020 0.7749 -0.2625 0.0658 -0.0114 0.0012 -0.0001
Using freqz(B,A), we get
Butterworth IIR filter example applied to reduce 50 Hz noise from ECG
Assume we the use the designed Butterworth filter to reduce 50 Hz noise from the ECG below
plot(ecg);
ecgf=filtfilt(B,A,ecg);
plot(ecgf);
