DSP Project 3 Report
Digital Filters BY
Hani Saleh
1
Table of Contents Table of Contents ............................................................................................................... 2 Table of Figures ................................................................................................................. 4 List of Tables ...................................................................................................................... 5 1
Abstract....................................................................................................................... 6
2
Introduction................................................................................................................ 6
3
Technical Background............................................................................................... 6 3.1
FIR Filters..................................................................................................................7
3.2
Designing FIR Filters using the Fourier Series Method..........................................7
3.2.1 Use of a window function..................................................................................................... 10 Rectangular window ........................................................................................................................... 10 Hann window...................................................................................................................................... 11 Hamming window............................................................................................................................... 11 Blackman window .............................................................................................................................. 11 3.2.2 Kaiser window...................................................................................................................... 11
3.3
Designing FIR Filters using the Frequency Sampling Method .............................12
3.4
Transformation of FIR filters transfer function....................................................14
3.5
IIR Filters.................................................................................................................15
3.5.1 3.5.1.1 3.5.2 3.5.3 3.5.3.1 3.5.3.2 3.5.4 3.5.5
4
Bilinear Transformation........................................................................................................ 15 Frequency Warping.......................................................................................................... 16 Design of Digital Filters using the bilinear transform .......................................................... 18 Analog IIR Filters Types ...................................................................................................... 18 Butterworth Filters........................................................................................................... 18 Chebyshev Filters ............................................................................................................ 20 Transformation of analog filters transfer function................................................................ 22 Other families of Analog filters [Ref-3] ............................................................................... 23
3.6
Filtering using Fast Fourier Transform.................................................................23
3.7
Comparison of IIR and FIR filters ........................................................................24
Implementation methodology:................................................................................. 24 4.1.1 FIR Filtering ......................................................................................................................... 24 4.1.1.1 FIR Using the Fourier series ............................................................................................ 24 4.1.1.2 FIR Using the Frequency Sampling Method ................................................................... 25 4.1.2 IIR Filtering using Butterworth and Chebyshev Type I ....................................................... 25 4.1.3 Filtering using Fast Fourier Transform................................................................................. 25
5
Computer Simulation............................................................................................... 25 5.1.1 The Input biomedical signals................................................................................................ 25 5.1.2 The Generated noise signals ................................................................................................. 28 5.1.3 Adding the noise and then filtering it out ............................................................................. 29 5.1.3.1 Filters specs...................................................................................................................... 30 5.1.3.2 Filtering using Fourier series (Dolph-Chebyshev Window) ............................................ 30 5.1.3.3 Filtering using Fourier series (Kaiser window )............................................................... 32 5.1.3.4 Filtering using the Frequency sampling method (Using the Kaiser window order) ........ 35 5.1.3.5 Filtering using the IIR Butterworth filter ......................................................................... 36 5.1.3.6 Filtering using the IIR Chebyshev Type I filter ............................................................... 38 5.1.3.7 Filtering using the FFT method ....................................................................................... 39 5.1.3.8 Results comparison Table ................................................................................................ 41
2
5.1.3.9
Conclusions:..................................................................................................................... 41
6
References ................................................................................................................ 41
7
Matlab Code ............................................................................................................. 42 7.1
Zero Order Bessel Function ...................................................................................42
7.2
Generate the Kaiser Window ....................................................................................42
7.3
Generate the Dolph-Chebyshev Window ..................................................................43
7.4
FIR Filtering using Fourier Series and Kaiser or Chebyshev windows (LP & HP).....45
7.5
FIR Filtering using Fourier Series and Kaiser or Chebyshev windows (BP & BS).....50
7.6
Create White noise ...................................................................................................54
7.7
FIR Filtering using the Frequency Sampling Approach (LP, HP, BP & BS) ..............55
7.8
Butterworth Normalized IIR Filter (Low pass)..........................................................58
7.9
Chebyshev Normalized IIR Filter (Low pass) ...........................................................61
7.10
Read the Biomedical signals .....................................................................................65
7.11
Filtering using FFT filter kernel function..............................................................67
7.12
Radix-2 FFT used by the FFT Filter......................................................................69
7.13 Add 60 Hz noise to the ECG Signal and then use all of the developed filters to remove it..............................................................................................................................72 7.14
Add white noise to the ECG signal and use all of the developed filters to remove it 82
7.15
Radix-2 IFFT used by the FFT Filter ....................................................................90
7.16
Test the FFT Filter (LP, HP, BP & BS) .................................................................90
7.17
Test the developed Filters.......................................................................................98
3
Table of Figures Figure 1: Ideal FIR Filter Implementation [Ref-3] ............................................................. 7 Figure 2: Ideal low pass filter response .............................................................................. 8 Figure 3: Sinc Function....................................................................................................... 8 Figure 4: - Frequency response of length 25 filter using rectangular window [Ref-3] ...... 9 Figure 5: Fig 4 re-plotted using dB gain [Ref-3] ................................................................ 9 Figure 6: Frequency response of length 25 filter using Hanning window [Ref-3]........... 10 Figure 7: Basic windows including a rectangular, Hamming, Hann, Blackman, and Kaiser (β=2.0, 5.0) along with their magnitude frequency responses. [Ref-5].................... 12 Figure 8: FIR design objective. [Ref-5]........................................................................... 13 Figure 9: Direct Synthesis Methodology [Ref-5] ............................................................. 13 Figure 10: Typical mirror and complement FIR filter’s relationship to a parent FIR [Ref5] ............................................................................................................................... 14 Figure 11: Relationship between the continuous- and discrete-frequency axes under the bi-linear z-transform, [Ref-5].................................................................................... 17 Figure 12: Design of a discrete-time IIR from an analog model using a bilinear ztransform. .................................................................................................................. 18 Figure 13: Ideal Lowpass Filter Frequency response [Ref-2] .......................................... 19 Figure 14: Butterworth Lowpass Filter Frequency response for n=1 -6 [Ref-2].............. 20 Figure 15: Chebyshev type I Lowpass Filter Frequency response for n=1 -6 [Ref-2] ..... 21 Figure 16: Chebychev type II Lowpass Filter Frequency response [Ref-2]..................... 22 Figure 17: MCL signal...................................................................................................... 26 Figure 18: Respiratory System signal ............................................................................... 26 Figure 19: Original Blood Pressure signal........................................................................ 27 Figure 20: Original ECG signal ........................................................................................ 27 Figure 21: White Noise signal .......................................................................................... 28 Figure 22: 60 Hz Noise signal .......................................................................................... 29 Figure 23: Noisy ECG signal (60 Hz noise added)........................................................... 30 Figure 24: FIR by Fourier Series (Dolph-Chebyshev window) Filter coefficients (N=375)..................................................................................................................... 31 Figure 25: FIR by Fourier Series (Dolph-Chebyshev window) Freq Response (N=375)31 Figure 26: FIR by Fourier Series (Dolph-Chebyshev window) filtering results (N=375) ................................................................................................................................... 32 Figure 27: FIR by Fourier Series (Kaiser Window) Filters Coefficients (N=229)........... 33 Figure 28: FIR by Fourier Series (Kaiser Window) Freq response (N=229) ................... 33 Figure 29: FIR by Fourier Series (Kaiser Window) Freq response (N=229) ................... 34 Figure 30: FIR by Fourier Series (Kaiser Window) filtering results (N=229) ................. 34 Figure 31: FIR by Freq Sampling Filter Freq response (N=229) ..................................... 35 Figure 32: FIR by Freq Sampling Filtering results (N=229) ............................................ 36 Figure 33: Butterworth Filter Freq response (N=25)........................................................ 37 Figure 34: Butterworth Filtering results (N=25)............................................................... 37 Figure 35: Chebyshev Filter Freq response (N=10) ......................................................... 38 Figure 36: Chebyshev Filtering result (N=10).................................................................. 39 Figure 36: FFT Filtering Results ...................................................................................... 40 Figure 36: FFT Filter Kernel............................................................................................. 40
4
List of Tables Table 1 –Windows effect of FIR filters designed by the Fourier series method [Ref-4] . 11 Table 2 - Transformations of a unit cutoff frequency lowpass filter to other specifications [Ref-3]....................................................................................................................... 23 Table 3 – Analog Filters response comparison................................................................. 23 Table 4 – All filters results comparison table ................................................................... 41
5
1 Abstract The main objective of this project is to study and design FIR and IIR filters; to achieve this goal FIR and IIR filters are designed and used to remove various types of noise from few biomedical signals obtained from the MIT database for biomedical signals.
2 Introduction The Digital Filter Design problem involves the determination of a set of filter coefficients to meet a set of design specifications. These specifications typically consist of the width of the passband and the corresponding gain, the width of the stopband(s) and the attenuation therein; the band edge frequencies (which give an indication of the transition band) and the peak ripple tolerable in the passband and stopband(s). Two types of digital filters exist according to their impulse response: • •
The FIR (Finite Impulse Response). The IIR (Infinite Impulse Response)
FIR filters possess certain properties, which make them the preferred design choices in numerous situations over IIR filters. Most notably, FIR filters (all zero system function) are always stable, with a realization existing for each FIR filter. Another feature exclusive to FIR filters is that of a linear phase response. The design of IIR filters is closely related to the design of analog filters, which is a widely studied topic. An analog filter is usually designed and a transformation is carried out into the digital domain. Two transformations exist – the impulse invariant transformation and the bilinear transformation. This project addresses the following topics: 1. Noise generation 2. FIR filters 3. FIR Filtering Using FFT 4. Butterworth and Chebyshev filters Matlab programs were developed to perform filtering on a real life biomedical signal from the MIT database. The filtration results were analyzed and then the analysis results are summarized in this report.
3 Technical Background A filter is a device which passes some signals 'more' than others. In this chapter a short background is given about designing FIR and IIR Digital filters.
6
3.1 FIR Filters The difference equation corresponding to the following block diagram x(n)
= unit delay
b(0)
b(M) y(n)
Figure 1: Ideal FIR Filter Implementation [Ref-3] is M
y (n ) = ∑ b(k )x(n − k ) k =0
This is called a feedforward or non-recursive, or transversal, or FIR filter, because it has a Finite Impulse Response (h (n) =0 for n>M). Such an Mth order FIR filter, has M+1 coefficients b(k), k=0,...,M. One of the M+1 degrees of freedom provided by M+ 1 coefficient is the overall filter gain. The other M determines the variation of gain with frequency. 3.2 Designing FIR Filters using the Fourier Series Method Given the desired frequency response H(Ω) of a filter, we can compute an appropriate inverse Fourier transform to obtain its ideal impulse response. Since the coefficients of an FIR filter equate to its impulse response, this would produce an “ideal” FIR filter. However, this “ideal” impulse is not actually constrained to be of finite length, and it may be non-causal (i.e. have non-zero response at negative time). Somehow we must generate an impulse response which is of limited duration, and causal. In more detail, H(Ω) (which is periodic in Ω with period 2π) can be represented using a Fourier series representation
H (Ω ) =
∞
∑ h(k )exp(− jkΩ)
k = −∞
in which the h(k) are computed from H(Ω) using the standard formula for Fourier series coefficients: π 1 h(k ) = H (Ω ) exp( jkΩ )dΩ 2π −∫π If the "ideal" filter coefficients h(k) are to be real-valued, then H(Ω) must be conjugate symmetric, i.e. H(-Ω) = H*(Ω) . We will consider the simplest case, a frequency response which is pure real (zero phase), and therefore symmetric about zero frequency.
7
For example, consider an ideal lowpass response, H(Ω)=1, |Ω|
Figure 2: Ideal low pass filter response The ideal filter coefficients (Fourier series coefficients) can in this case be calculated: π 1 A sin (kA) h(k ) = H (Ω ) exp( jkΩ )dΩ = ∫ 2π −π π kA
This 'Sinc' response is symmetric about sample k=0, and infinite in extent: 1
0.8
0.6
0.4
0.2
0
-0.2
-0.4 -10
-8
-6
-4
-2
0
2
4
6
8
10
Figure 3: Sinc Function
To implement an order-M FIR filter, assume we select only a finite length section of d(k). For the Sinc response shown above, the best section to select (that is, the one which gives minimum total squared error) is symmetric about 0, i.e. h(k), -M/2 ≤ k ≤ M/2. The resulting filter is non-causal, but it can be made causal simply by shifting the response (adding delay). This selection operation is equivalent to multiplying the filter coefficients by a rectangular window extending from -M/2 to M/2. 8
We can compute the resulting filter frequency response, which is a “truncated Fourier series approximation” of H(Ω), given by B(Ω ) =
A
π
M 2
∑
k =− M
sin (kA) exp(− jkΩ ) kA 2
This truncated response is illustrated in Figure 6.1 for the case M=24 (length 25) and A= π /2 (cutoff frequency = 0.25 x sample frequency). Note the well known Gibb's phenomenon (an oscillatory error, increasing in magnitude close to any discontinuities in H(Ω) ). 1.5
D
1 0.5 0 -0.5
0
0.5
1
1.5
2
2.5
3
Normalised radian frequency
Figure 4: - Frequency response of length 25 filter using rectangular window [Ref-3] Because the filter is still symmetric about zero frequency, the response in Figure 4 is still zero phase. The actual filter would require an added delay of M/2 samples, which does not affect the amplitude response, but introduces a linear phase term to the frequency response. Let us re-plot Figure 5 using a dB amplitude scale. The sidelobes due to the rectangular window can be clearly seen: 20
Mainlobe
First sidelobe
0
Sidelobes
|D| (dB) -20 -40 -60 -80
0
0.5
1
1.5
2
2.5
3
Normalised radian frequency
Figure 5: Fig 4 re-plotted using dB gain [Ref-3]
9
The high sidelobe level close to the passband, and the slow decay of sidelobe level away from the passband, make this an unsatisfactory response for most purposes.
3.2.1
Use of a window function
To reduce the ripples in both of the passband and the stop band we create the required finite number of filter coefficients by multiplying the infinite-length coefficient vector h(k) by a finite-length window with non-rectangular shape. Consider a raised cosine (Hann or Hanning) window function 0.5 1+ cos(kπ (1 + M 2 )) , giving the filter frequency response
(
B(Ω ) =
)
A
π
kπ sin (kA) 1 1 + cos exp(− jkΩ ) 22 1 + M 2 kA
M 2
∑
k =−M
(6.5)
Which is illustrated in Figure 6, again for the case of length 25 and A= π /2. 1.5 1 D 0.5 0 -0.5
0
0.5
1
1.5
2
2.5
3
Normalised radian frequency
Figure 6: Frequency response of length 25 filter using Hanning window [Ref-3] Several window functions have been devised, which offer different tradeoffs of transition width, sidelobe level, and so on. As well as the rectangular window and the Hann or Hanning window, common window functions are the Hamming window, the Blackman window, and the Kaiser window. The latter includes a 'ripple control' parameter ß which allows the designer to tradeoff passband ripple against transition width. Their discrete-time window specifications are summarized below. Rectangular window
A rectangular window is given by
(N − 1) 1if k ≤ w [k ] = 2 0otherwise 10
Hann window
A Hann window is given by
2πk (N − 1) if k ≤ 0.5 + 0.5 cos w [k ] = 2 (N − 1) 0otherwise
Hamming window
A Hamming window is given by
2πk (N − 1) if k ≤ 0.54 + 0.46 cos w [k ] = 2 (N − 1) 0otherwise
Blackman window
A Blackman window is given by
2πk 4πk (N − 1) + 0.08 cos if k ≤ 0.42 + 0.5 cos w [k ] = 2 (N − 1) (N − 1) 0otherwise
3.2.2 Kaiser window A Kaiser window is given by 2 I 0 β 1 − 2k N − 1 w[k ] = if k ≤ ( N − 1) I0β 2 0 otherwise where I0(β) is a 0th order Bessel function. Choosing a suitable window function can be done with the aid of published data such as this that shown in Table 1. Name of Transition Passband ripple Main lobe Maximum window width/ sample (dB) relative to side stopband function frequency lobe (dB) attenuation (dB) Rectangular 0.9 / N 0.75 13 21 Hann(ing) 3.1/N 0.055 31 44 Hamming 3.3/N 0.019 41 53 Blackman 5.5/N 0.0017 57 74 0.0274 50 Kaiser (β=4.54) 2.93/N 0.000275 90 Kaiser (β=8.96) 5.71/N Table 1 –Windows effect of FIR filters designed by the Fourier series method [Ref-4]
11
Figure 7: Basic windows including a rectangular, Hamming, Hann, Blackman, and Kaiser (β=2.0, 5.0) along with their magnitude frequency responses. [Ref-5]
3.3 Designing FIR Filters using the Frequency Sampling Method Normally, FIRs are designed to achieve a pre-specified frequency domain response. If |Hd(ejω)| denotes a desired magnitude frequency response, then the design objective becomes one of deriving an FIR with a frequency response |H(ejω)| which closely approximates |Hd(ejω)| in some acceptable manner. The mathematical metrics used to compare the desired to the realized filter response is generally defined in terms of a frequency dependent error
( )
( )
e(ϖ ) = H e jϖ − Hd e jϖ
The criterion of optimization normally takes on one of the following forms: •
Minimum squared error (MSE) criterion: minimize φ (ϖ ) =
∑ω e(ϖ )
2
∀
•
Minimax error criterion: minimize (maximum( ε (ϖ )));∀ω
Both design criteria are motivated in Figure 8. The advantage of the MSE method is that it is generally easy to compute a solution when applied to the design of an LTI filter. The problem with MSE method is that local errors can be very large even though the overall MSE error is small. The minimax error criterion is usually more difficult to satisfy but will guarantee that the worst case error has been reduced to a quantifiable value. This, if it can be achieved, is generally desired.
12
Figure 8: FIR design objective. [Ref-5] The simplest FIR synthesis technique is called the Frequency sampling direct method. It is based on the use of the DFT and IDFT. It involves computing the M-sample inverse DFT (IDFT) of a specified desired frequency response Hd(ejω). An Nth-order Direct FIR, where N≤M, is defined by an N-sample time-series extracted from the M-sample IDFT, symmetrically oriented about the center tap. The Direct FIR filter design procedure is shown in Figure 8.
Figure 9: Direct Synthesis Methodology [Ref-5] The z-transform representation of the realized FIR filter is given by N -1
H (z ) =
∑ h[k ]z
-k
k =0
13
3.4 Transformation of FIR filters transfer function The transfer function of an Nth order FIR has the form: M
H (z ) =
∑h
mz
−m
6
m = −M
where N=2M+1. It is often desired to develop variations of these filters which maintain a close mathematical relationship to the original but exhibit different frequency-domain behavior. Examples of this would be the creation of a filter which mirrors the frequency response of another. If the original, or parent filter is lowpass, having a normalized passband width ∆ relative to the sampling frequency, then a mirror filter would be a highpass filter also having a bandwidth ∆. This can be achieved simply by modulating the impulse response h[k] by a sinusoid running at the Nyquist frequency, namely c[k]=cos(2π(N/2)k/N)= cos(πk)=(-1)k. This action will translate the frequency response of the original FIR up to the Nyquist frequency. Therefore, the mirror version of an FIR having an impulse response h[k], is simply given by hmirror [k ] = ( −1)k h[k ]
7
Physically, the only difference between the mirror and the original is the alternating sign on the tap-weight coefficients. The spectral relationship between a mirror and its parent is graphically interpreted in Figure 7.
Figure 10: Typical mirror and complement FIR filter’s relationship to a parent FIR [Ref5] Another useful filter type is called the complement filter. A complement filter is mathematically related to the original odd order H(z) through H comp (e jω ) + H (e jω ) = 1
That is, the original and complement filter, when additively combined, form an all-pass filter. This can be particularly useful in designing banks of filter which must cover a
14
frequency range but must do so with a set of band limited filters. The implementation of a complement filter is very simple. Consider again Equation 8 rewritten as Hcomp (e jω ) = 1 − H (e jω )
which simply states the center tap coefficient h[0] of the original FIR is subtracted from unit to form the center tap for the complement FIR. The other complement FIR coefficients are simply the negation of the original. If the initial sample of the complement FIR’s impulse response to be indexed at k=0, rather than h[-M] as in Equation 9, then Equation 9 can be expressed as Hcomp_M ( z ) = z M Hcomp ( z ) = z −M (1 − H ( z )) = z −M − H0 ( z )
which states that the compensated filter can be obtained simply with the addition of a Mdelay shift register which is a multiplier less operation.
3.5 IIR Filters Digital Filters are designed by using the values of both the past outputs and the present input, an operation brought about by convolution. If such a filter is subjected to an impulse then its output need not necessarily become zero. The impulse response of such a filter can be infinite in duration. Such a filter is called an Infinite Impulse Response filter or IIR filter. The infinite impulse response of such a filter implies the ability of the filter to have an infinite impulse response. This indicates that the system is prone to feedback and instability. IIR filters are designed essentially by many methods. The most useful method in practice is the bilinear transform
3.5.1 Bilinear Transformation The Bilinear Transformation method overcomes the effect of aliasing that is caused to due the analog frequency response containing components at or beyond the Nyquist Frequency. The bilinear transform is a method of compressing the infinite, straight analogue frequency axis to a finite one long enough to wrap around the unit circle once only. This is also sometimes called frequency warping. This introduces a distortion in the frequency. This is undone by pre-warping the critical frequencies of the analog filter (cutoff frequency, center frequency) such that when the analog filter is transformed into the digital filter, the designed digital filter will meet the desired specifications. Consider an analog filter: b H (S ) = S+a
15
Then this filter could be mapped into the digital domain using:
H ( z ) = H ( s)
s=
2 × ( z − 1) T × ( z + 1)
This transformation is known as the Bilinear or Tustin Transformation. The Laplace Transforms in the filter expressions are replaced by the corresponding z-transforms.
Replacing s = σ + j Ω and performing algebraic manipulations, substituting z = ejw we get
ω = 2 tan-1(ΩT/2) It can be seen that analog dc (s = 0) maps to digital dc (z = 1) and the highest analog frequency (s = ∞) maps to the highest digital frequency (z = -1). It is easy to show that the entire jw axis in the s plane is mapped exactly once around the unit circle in the z plane. Therefore, it does not alias. With (2/T) real and positive, the left-half s plane maps to the interior of the unit circle, and the right-half s plane maps outside the unit circle. The constant provides one remaining degree of freedom that can be used to map any particular finite frequency the jw axis in the s plane to a particular desired location on the unit circle ejw in the z plane. All other frequencies will be warped. In particular, approaching half the sampling rate, the frequency axis compresses more and more. Filters having a single transition frequency, such as lowpass or highpass filters, map beautifully under the bilinear transform; you simply map the cut-off frequency where it belongs, and the response looks great. In particular, ``equal ripple'' is preserved for optimal filters of the elliptic and Chebyshev types because the values taken on by the frequency response are identical in both cases; only the frequency axis is warped. 3.5.1.1 Frequency Warping
Assume that the frequency response of an analog filter is defined by Ha(jΩ). The frequency response of a digital filter is given by H(ejω),where ϖ∈[-π,π]. The discretetime frequency axis can be calibrated with respect to the actual sampling frequency using ϖ←ω/Ts. The problem at hand is to construct a mapping of the jΩ-axis in the s-plane to the unit circle in the z-plane, in a manner which eliminates aliasing. The bilinear ztransform is capable of achieving this objective. To demonstrate this claim, consider evaluating the Bilinear Transform Equation for a given s=jΩ and z=ejϖ, namely
16
2 e jϖ − 1 2 j sin(ϖ / 2) 2 jΩ = jϖ = = j tan(ϖ / 2) Ts e + 1 Ts cos(ϖ / 2) Ts
This, upon simplification becomes: 2 Ω = tan(ϖ / 2) or ϖ = tan−1(ΩTs / 2) Ts
This is interpreted in Figure 11.
Figure 11: Relationship between the continuous- and discrete-frequency axes under the bi-linear z-transform, [Ref-5]. The relationships found in the above equations are called the warping equation (Ω→ϖ) and pre-warping equation (ϖ→Ω). It can also be seen that the relationship between the analog- and discrete-frequency axes is non-linear. As a result, the critical frequencies of an analog filter model will not in general align themselves exactly with the critical frequencies of a bilinear z-transformed discrete-time filter. What is of significant importance, however, is noting that the bilinear z-transform overcomes the leakage or aliasing problem associated with the impulse invariant design method. It can also be seen that as Ω→∞, ϖ→π found at the end of the unit circle. As a result, no aliasing is possible in a bilinear-z transformed spectrum.
17
Because of this, the bilinear z-transform is well suited to the task of converting classic analog filter into a digital IIR model which preserves the shape of the magnitude frequency response of the parent analog system.
3.5.2
Design of Digital Filters using the bilinear transform
The design paradigm is, unfortunately, not a straightforward process due to the fact that the analog prototype is defined in terms of a set of analog frequencies that differ from those possessed by the target digital filter. The difference is governed by the warping equations. The possible distortion between the frequency axes can, however, be accounted for by using the procedure outlined below and summarized graphically in Figure 12: 1. Specify the desired digital filter specifications. 2. Pre-warp the digital critical frequencies into corresponding analog frequencies. From the derived set of analog filter specifications, estimate analog filter order. 3. Design the analog prototype filter H(s). 4. Convert the analog prototype H(s) into a target analog filter H(s) using frequencyfrequency transforms. 5. Design a digital filter H (z) using a bilinear z-transform of H(s). The pre-warping introduced in Step 2 is corrected by the warping introduced by the bilinear z-transform. While this method may initially seem to be complicated, it is a sequential procedure that can be reduced to a digital computer program.
Figure 12: Design of a discrete-time IIR from an analog model using a bilinear ztransform.
3.5.3
Analog IIR Filters Types
3.5.3.1 Butterworth Filters The ideal low pass filter frequency response is as shown in Figure 1.
18
Figure 13: Ideal Lowpass Filter Frequency response [Ref-2] The Butterworth response is:
H ( jw)
2
=
Ho w 1 + wc
2n
This is known as the nth order Butterworth low pass response. A plot of the response for many values of n is shown in Figure 2.
19
Figure 14: Butterworth Lowpass Filter Frequency response for n=1 -6 [Ref-2]
Note that the more we increase n the closer we get to the ideal low pass filter. Butterworth filters are causal in nature and of various orders, the lowest order being the best (shortest) in the time domain, and the higher orders being better in the frequency domain. Butterworth or maximally flat filters have a monotonic amplitude frequency response which is maximally flat at zero frequency response (Figure 2) and the amplitude frequency response decreases logarithmically with increasing frequency. The Butterworth filter has minimal phase shift over the filter's band pass when compared to other conventional filters. The general expression for the transfer function of an n-th order Butterworth lowpass filter is given by: 1 1 H ( s) = n = ∏ i (1 + s i ) (1 + s1 )(1 + s 2 )(1 + s 3 )....(1 + s n )
Where
2i + n − 1 2i + n − 1 si = cos π + j sin π 2n 2n
3.5.3.2 Chebyshev Filters Chebyshev filters are of two types: Chebyshev I filters are all pole filters which are equiripple in the passband and are monotonic in the stopband, the frequency response is shown Figure 3.
20
Figure 15: Chebyshev type I Lowpass Filter Frequency response for n=1 -6 [Ref-2] The transfer function for it is given by: 2 1 H ( jw) = 2 1 + ε × Tn2 ( w)
ε 2 = 10 r / 10 − 1 Tn ( w) = Chebychev polynomial ororder n r = pasband ripple in dB The general expression for the transfer function of an n-th order Chebychev I lowpass filter is given by: H Ho H (s) = n o = ∏ i (1 + si ) (1 + s1 )(1 + s2 )(1 + s3 )....(1 + sn ) Where
si = σ i + jωi
1 − γ γ σi = 2
(2i − 1) × sin π 2n
21
1 − γ γ ωi = 2
(2i − 1) × cos π 2n
1+ 1+ ε 2 γ = ε
r
1/ n
ε = 1010 − 1
Chebyshev II filters contain both poles and zeros exhibiting a monotonic behavior in the passband and equiripple in the stopband, a sample frequency response is shown in Figure 4.
Figure 16: Chebychev type II Lowpass Filter Frequency response [Ref-2]
3.5.4
Transformation of analog filters transfer function
An analogue lowpass filter with cutoff frequency ωc = 1 and Laplace variable s' may be transformed to other filters, with Laplace variable s, as shown in Table 2.
22
Lowpass transformed to: Lowpass Highpass Bandpass Bandstop
Use substitution s' → s' →
New cutoff (s)
s
ωc
ωc ωc
ωc
s s 2 + ω lω u s' → s (ω u − ω l ) s (ω − ω l ) s' → 2 u s + ω lω u
(Use for changing cutoff frequency)
ωl lower, ωu upper
ωl lower, ωu upper
Table 2 - Transformations of a unit cutoff frequency lowpass filter to other specifications [Ref-3]
3.5.5
Other families of Analog filters [Ref-3]
There is a wide range of closed form analogue filters. Some are all-pole; others have zeros. Some have monotonic responses; some equiripple. For a given band edge frequency, ripple specification, and filter order, narrower transition bandwidth can be traded off against worse phase linearity (i.e. more group delay variation), as summarized in Table 3 below: Filter family Bessel Optimal (L) Butterworth Chebyshev type I Chebyshev type II Elliptic
Passband response type monotonic monotonic monotonic equiripple monotonic equiripple
Stopband response type monotonic monotonic monotonic monotonic equiripple equiripple
Table 3 – Analog Filters response comparison 3.6 Filtering using Fast Fourier Transform The DFT generate the spectrum of the input signal; the frequency content of the input signal will be present in the signal.
Digital filtering can be easily done with the FFT spectrum, by: 1) Calculate the FFT of the input signal.
23
2) Chang the FFT results values which correspond to the frequencies you wish to change. 3) Perform IFFT to get the filtered time domain signal. 3.7
Comparison of IIR and FIR filters
If the desired filter is highly selective (that is, its frequency response has small transition bandwidths or "steep sides"), then the impulse response will be long in the time domain. Examples include narrowband filters and lowpass /highpass /bandpass filters with steep cutoffs. For an FIR filter, a long impulse response means the filter is long (high order), so it requires many multiplications, additions and delays per sample. An IIR filter has active poles as well as zeros. Poles, acting as high-Q resonators, can provide highly selective frequency responses (hence long impulse responses) using much lower filter order than the equivalent FIR filter, hence much lower computational cost. Although it is still true that a more selective response requires a higher order filter. On the other hand, the closer to linear the phase is required to be, the higher the order of IIR filters that is needed. Also the internal word lengths in IIR filters need generally to be higher than those in FIR filters; this may increase the implementation cost (e.g in VLSI). An FIR filter is inherently stable, unlike an IIR filter. Hence an FIR implementation involving inaccurate (finite precision, or 'quantized') coefficients will be stable, whereas an IIR one might not. (However, it is desirable in either case to compute the actual frequency response of the filter, using the actual quantized values of the coefficients, to check the design.)
4 Implementation methodology: This section explains the algorithm used to design the filters. 4.1.1
FIR Filtering
4.1.1.1 FIR Using the Fourier series
Algorithm: 1. Compute the required filter Order 2. Use the equations provided in section 3.2 compute the filter coefficients 3. Compute the selected window values 4. Multiply the window by the coefficients 5. Do the filtering
24
4.1.1.2 FIR Using the Frequency Sampling Method
Algorithm: 1. Decide on the filter order (Use the Kaiser window computed order) 2. Build your kernel in the Frequency domain 3. Compute the IFFT of the filter kernel with number of points equal to the filter order 4. Do the filtering 4.1.2 IIR Filtering using Butterworth and Chebyshev Type I Algorithm: 1. Specify the desired digital filter specifications. 2. Pre-warp the digital critical frequencies into corresponding analog frequencies. From the derived set of analog filter specifications, estimate analog filter order. 3. Design the analog prototype filter H(s). 4. Convert the analog prototype H(s) into a target analog filter H(s) using frequencyfrequency transforms. 5. Design a digital filter H (z) using a bilinear z-transform of H(s). The pre-warping introduced in Step 2 is corrected by the warping introduced by the bilinear ztransform. 4.1.3 Filtering using Fast Fourier Transform Since the FFT generate the spectrum of the input signal; it was used to perform digital filtering by getting ride off the unneeded frequency components.
Digital filtering using FFT was performed as follows: • Calculate the FFT of the input signal. • Generate a frequency response in the frequency domain that corresponds to the Filter requirements • Multiply the filter vector by the FFT result of the input signal • Perform IFFT to get the filtered time domain signal.
5 Computer Simulation In this section I summarize the computer simulation results. 5.1.1 The Input biomedical signals The following figures shows the time domain and the spectrum of the four biomedical signals used for this project simulation:
25
The MCL signal Time domain plot 0.4
Amplitude (mV)
0.2
0
-0.2
-0.4
-0.6
10
10.5
11
11.5
12 Time (sec)
12.5
13
13.5
14
The spectrum of the MCL Signal 400 350 300
MCL Signal
Amplitude
250 200 150 100 50 0 -80
-60
-40
-20
0 Freq. (Hz)
20
40
60
80
Figure 17: MCL signal
The RESP signal Time domain plot 0.5 0.4 Respiratory System Signal 0.3
Amplitude (mV)
0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 7
7.5
8
8.5 Time (sec)
9
9.5
10
The spectrum of the RESP Signal 250
Amplitude
200
150
100
50
0 -80
-60
-40
-20
0 Freq. (Hz)
20
40
60
80
Figure 18: Respiratory System signal
26
The ABP signal Time domain plot 40 Blood Pressure Original SignalSignal
Amplitude (mmhg)
30
20
10
0
-10
7.5
8
8.5
9 Time (sec)
9.5
10
10.5
The spectrum of the ABP Signal 18000 16000 14000
10000 8000 6000 4000 2000 0 -80
-60
-40
-20
0 Freq. (Hz)
20
40
60
80
Figure 19: Original Blood Pressure signal The ECG signal Time domain plot 1.5
Original ECG Signal
Amplitude (mV)
1
0.5
0
-0.5 13
14
15
16
17
18
19
Time (sec)
The spectrum of the ECG Signal 70 60 50 Amplitude
Amplitude
12000
40 30
20 10
0 -80
-60
-40
-20
0 Freq. (Hz)
20
40
60
80
Figure 20: Original ECG signal
27
5.1.2 The Generated noise signals The following figures shows the time domain and the spectrum for the generated white noise and 60Hz noise signals. The White Noise signal Time domain plot 1.2 1
White Noise
0.8
Amplitude (mV)
0.6 0.4 0.2 0 -0.2 -0.4 -0.6 0
5
10
15
20
25
30
Time (sec)
The spectrum of the Noise Signal 40 35 30
Amplitude
25 20 15 10 5 0 -80
-60
-40
-20
0 Frequency (Hz)
20
40
60
80
Figure 21: White Noise signal
28
The 60 Hz Noise signal Time domain plot 1 0.8 0.6
Amplitude (mV)
0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1
0
0.02
0.04
0.06
0.08
0.1 Time (sec)
0.12
0.14
0.16
0.18
0.2
The spectrum of the Noise Signal 2500
2000
Amplitude
60 Hz Noise Signal 1500
1000
500
0 -80
-60
-40
-20
0 Frequency (Hz)
20
40
60
80
Figure 22: 60 Hz Noise signal 5.1.3 Adding the noise and then filtering it out The following figure shows the time domain and the spectrum of the ECG signal after adding the noise.
29
The ECG signal with 60 Hz Noise Time domain plot 1.5
Amplitude (mV)
1
0.5
0
-0.5
7
7.5
8
8.5
9
9.5
10
10.5
11
11.5
Time (sec)
The spectrum of the ECG signal + Noise 450 400 The Noisy ECG Signal (60 Hz noise added)
350
Amplitude
300 250 200 150 100 50 0 -80
-60
-40
-20
0 Frequency (Hz)
20
40
60
80
Figure 23: Noisy ECG signal (60 Hz noise added)
5.1.3.1 Filters specs
The requirements for all of the designed filters were: F_pass = 40 Hz F_stop = 45 Hz Dp = 3 dB Ds = 70 dB Fs = 125 Hz 5.1.3.2 Filtering using Fourier series (Dolph-Chebyshev Window)
The following figures show the frequency response, pole-zero plot, phase response and the filtered signal compared to the original signal.
30
Poles & Zeros plot
Filter Coeifficints after windowing 0.7
1
0.6 0.5 0.4 Amplitude
Imag
0.5
374
0
0.3 0.2 0.1
-0.5
0 -0.1
-1 -1.5
-1
-0.5
0 Real Part
0.5
1
-0.2
1.5
The Selected winow 1 0.9
100
150
200 Time
250
300
350
400
300
350
400
Impulse response after windowing
0.5
0.7
0.4
0.6
Amplitude
Amplitude
50
0.6
0.8
0.5 0.4
0.3 0.2 0.1
0.3
0
0.2
-0.1
0.1 0
0
ERROR Measurements for FIR Filtering using Fourier Series (Dolph-Cheby window) Error-Max = 0.0683814 Error-RMS = 0.00167725 SNR = 43.4112 N = 375 0.7
0
50
100
150
200 Time
250
300
350
-0.2
400
0
50
100
150
200 Time
250
Figure 24: FIR by Fourier Series (Dolph-Chebyshev window) Filter coefficients (N=375) Freq response befor windowing
Freq response after windowing
10
50
0 0
-10
-50 Amplitude (dB)
Amplitude (dB)
-20 -30 -40 -50
ERROR Measurements for FIR Filtering using Fourier Series (Dolph-Cheby window) Error-Max = 0.0683814 Error-RMS = 0.00167725 SNR = 43.4112 N = 375
-100
-150
-60 -70
-200
-80 -90
0
0.1
0.2
0.3
0.4
0.5 Freq
0.6
0.7
0.8
0.9
1
-250
0
0.1
0.2
4
3
3
2
2
1 0 -1
-3
0.2
0.3
0.4
0.5 Freq
0.6
0.7
0.7
0.8
0.9
1
0.8
0.9
1
0
-2
0.1
0.6
-1
-3
0
0.5 Freq
1
-2
-4
0.4
Phase response after windowing
4
Phase in Radians
Phase in Radians
Phase response befor windowing
0.3
0.8
0.9
1
-4
0
0.1
0.2
0.3
0.4
0.5 Freq
0.6
0.7
Figure 25: FIR by Fourier Series (Dolph-Chebyshev window) Freq Response (N=375)
31
The Noisy ECG signal Time domain plot
The Filtered ECG signal Time domain plot
3 3
2.5
2.5
2
2 Amplitude (mV)
Amplitude (mV)
1.5 1 0.5 0
1.5 1 0.5 0
-0.5
-0.5
-1 -1.5
ERROR Measurements for FIR Filtering using Fourier Series (Dolph-Cheby window): Error-Max = 0.0683814 Error-RMS = 0.00167725 SNR = 43.4112 N = 375
-1 0
0.2
0.4
0.6
0.8
1 1.2 Time (sec)
1.4
1.6
1.8
2
1.6
1.8
The spectrum of the ECG Noisy Signal
2
2.2
2.4 2.6 Time (sec)
2.8
3
3.2
3.4
The spectrum of the Filtered Noise Signal
450
70
400
60
350 50 Amplitude
Amplitude
300 250 200
40 30
150 20 100 10
50 0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
Figure 26: FIR by Fourier Series (Dolph-Chebyshev window) filtering results (N=375) ERROR Measurements for FIR Filtering using Fourier Series (Cheby window) Error-Max = 0.0683814 Error-RMS = 0.00167725 SNR = 43.4112 N = 375 5.1.3.3 Filtering using Fourier series (Kaiser window )
The following figures show the frequency response, pole-zero plot, phase response and the filtered signal compared to the original signal.
32
Poles & Zeros plot
Filter Coeifficints after windowing 0.7
1
0.6 Kaiser Window, N = 229 0.5 0.4 Amplitude
Imag
0.5
228
0
0.3 0.2 0.1
-0.5
0 -0.1
-1 -1.5
-1
-0.5
0 Real Part
0.5
-0.2
1
0
50
The Selected winow
200
250
200
250
Impulse response after windowing
1
0.7 0.6
0.8
0.5
0.7
0.4 Amplitude
0.6
Amplitude
150 Time
0.9
0.5 0.4
0.3 0.2 0.1
0.3
0
0.2
-0.1
0.1 0
100
0
50
100
150
200
-0.2
250
0
50
100
Time
150 Time
Figure 27: FIR by Fourier Series (Kaiser Window) Filters Coefficients (N=229) Freq response befor windowing
Freq response after windowing
20
20 0
0 -20
Fourier Series, Order = 229 -40
Kaiser Window, Order = 229
-40
Amplitude (dB)
Amplitude (dB)
-20
-60 -80
-60 -100 -80 -120 -100
0
0.1
0.2
0.3
0.4
0.5 Freq
0.6
0.7
0.8
0.9
1
-140
0
0.1
0.2
4
3
3
2
2
1 0 -1
-3
0.2
0.3
0.4
0.5 Freq
0.6
0.7
0.7
0.8
0.9
1
0.8
0.9
1
0
-2
0.1
0.6
-1
-3
0
0.5 Freq
1
-2
-4
0.4
Phase response after windowing
4
Phase in Radians
Phase in Radians
Phase response befor windowing
0.3
0.8
0.9
1
-4
0
0.1
0.2
0.3
0.4
0.5 Freq
0.6
0.7
Figure 28: FIR by Fourier Series (Kaiser Window) Freq response (N=229)
33
The Noisy ECG signal Time domain plot
The Filtered ECG signal Time domain plot
3 3
2.5
2 Amplitude (mV)
Amplitude (mV)
1.5 1 0.5 0
1.5 1 0.5 0
-0.5
-0.5
-1 -1.5
ERROR Measurements for FIR Filtering using Fourier Series (Kaiser Window) Error-Max = 0.0683609 Error-RMS = 0.0016757 SNR = 43.45 N = 229
2.5
2
-1 0
0.2
0.4
0.6
0.8
1 1.2 Time (sec)
1.4
1.6
1.8
2
1
1.2
The spectrum of the ECG Noisy Signal
1.4
1.6
1.8 2 Time (sec)
2.2
2.4
2.6
2.8
The spectrum of the Filtered Noise Signal
450
70
400
60
350 50 Amplitude
Amplitude
300 250 200
40 30
150 20 100 10
50 0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
0 -80
80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
Figure 29: FIR by Fourier Series (Kaiser Window) Freq response (N=229) The Noisy ECG signal Time domain plot
The Filtered ECG signal Time domain plot
3 3
2.5
2 Amplitude (mV)
Amplitude (mV)
1.5 1 0.5 0
1.5 1 0.5 0
-0.5
-0.5
-1 -1.5
ERROR Measurements for FIR Filtering using Fourier Series (Kaiser Window) Error-Max = 0.0683609 Error-RMS = 0.0016757 SNR = 43.45 N = 229
2.5
2
-1 0
0.2
0.4
0.6
0.8
1 1.2 Time (sec)
1.4
1.6
1.8
2
1
1.2
The spectrum of the ECG Noisy Signal
1.4
1.6
1.8 2 Time (sec)
2.2
2.4
2.6
2.8
The spectrum of the Filtered Noise Signal
450
70
400
60
350 50 Amplitude
Amplitude
300 250 200
40 30
150 20 100 10
50 0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
Figure 30: FIR by Fourier Series (Kaiser Window) filtering results (N=229)
34
ERROR Measurements for FIR Filtering using Fourier Series (Kaiser Window) Error-Max = 0.0683609 Error-RMS = 0.0016757 SNR = 43.45 N = 229 5.1.3.4 Filtering using the Frequency sampling method (Using the Kaiser window order)
The following figures show the frequency response, pole-zero plot, phase response and the filtered signal compared to the original signal. FIR (Freq-Sampling) Poles & Zeros plot
FIR (Freq-Sampling) Freq response of the filter
1
0
0.8 -10
0.6 Amplitude (dB)
0.4 Imag
0.2 228
0 -0.2
-20
ERROR Measurements for FIR Filtering using Frequency Sampling, (Same order as Kaiser window order): Error-Max = 0.0661603 Error-RMS = 0.00165881 SNR = 43.5344 N = 229
-30 -40
-0.4 -50
-0.6 -0.8
-60
-1 -70 -1
-0.5
0 Real Part
0.5
1
1.5
0.1
FIR (Freq-Sampling) Impulse response of the filter 4
0.6
3
0.5
0.4 Freq x π
0.5
0.6
0.7
2 Phase in Radians
0.4 Amplitude
0.3
FIR (Freq-Sampling) Phase response of the filter
0.7
0.3 0.2 0.1
1 0 -1 -2
0
-3
-0.1 -0.2
0.2
0
50
100
150 Time
200
250
-4
0
0.1
0.2
0.3
0.4
0.5 Freq x π
0.6
0.7
0.8
0.9
1
Figure 31: FIR by Freq Sampling Filter Freq response (N=229)
35
The Noisy ECG signal Time domain plot
The Filtered ECG signal Time domain plot
3
ERROR Measurements for FIR Filtering using Frequency Sampling (Same order as Kaiser window order) Error-Max = 0.0661603 Error-RMS = 0.00165881 SNR = 43.5344 N = 229
3 2.5 2.5 2 Amplitude (mV)
Amplitude (mV)
2 1.5 1 0.5
1.5 1 0.5
0 0 -0.5 -0.5 -1 0.2
0.4
0.6
0.8
1 1.2 Time (sec)
1.4
1.6
1.8
-1
2
1
1.2
The spectrum of the ECG Noisy Signal
1.4
1.6
1.8 2 Time (sec)
2.2
2.4
2.6
2.8
The spectrum of the Filtered Noise Signal
450
70
400
60
350 50 Amplitude
Amplitude
300 250 200
40 30
150 20 100 10
50 0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
Figure 32: FIR by Freq Sampling Filtering results (N=229) ERROR Measurements for FIR Filtering using Frequency Sampling (Same as Kaiser order Error-Max = 0.0661603 Error-RMS = 0.00165881 SNR = 43.5344 N = 229 5.1.3.5 Filtering using the IIR Butterworth filter
The following figures show the frequency response, pole-zero plot, phase response and the filtered signal compared to the original signal.
36
Butterworth Poles & Zeros plot
Butterworth Freq response of the filter 0
1 0.8
-50
0.6 0.4 Amplitude (dB)
-100
Imag
0.2 0 -0.2
ERROR Measurements for FIR Filtering using Butterworth: Error-Max = 0.0451924 Error-RMS = 0.00154873 SNR = 44.2513 N = 25
-150
-200
-0.4 -0.6
-250
-0.8 -1 -1.5
-1
-0.5 0 Real Part
0.5
-300
1
0
0.05
0.1
Butterworth Impulse response of the filter
0.4
3
0.3
2 Phase in Radians
4
Amplitude
0.2 0.1 0
-3
60
80
100
120
0.3
0.35
0
-2
40
0.35
-1
-0.2
20
0.3
1
-0.1
0
0.25
Butterworth Phase response of the filter
0.5
-0.3
0.15 0.2 Freq x π
140
160
-4
180
0
0.05
0.1
0.15 0.2 Freq x π
Time
0.25
Figure 33: Butterworth Filter Freq response (N=25) The Noisy ECG signal Time domain plot
The Filtered ECG signal Time domain plot
3
3 2.5
2
2
1.5
1.5 Amplitude (mV)
Amplitude (mV)
2.5
ERROR Measurements for FIR Filtering using Butterworthr Error-Max = 0.0451924 Error-RMS = 0.00154873 SNR = 44.2513 N = 25
1 0.5
1 0.5
0
0
-0.5
-0.5
-1
-1
-1.5
-1.5
0
0.2
0.4
0.6
0.8
1 1.2 Time (sec)
1.4
1.6
1.8
2
0
0.2
The spectrum of the ECG Noisy Signal
0.4
0.6
0.8
1 1.2 Time (sec)
1.4
1.6
1.8
2
The spectrum of the Filtered Noise Signal
450
70
400
60
350 50 Amplitude
Amplitude
300 250 200
40 30
150 20 100 10
50 0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
Figure 34: Butterworth Filtering results (N=25)
37
ERROR Measurements for FIR Filtering using Butterworth Error-Max = 0.0451924 Error-RMS = 0.00154873 SNR = 44.2513 N = 25
5.1.3.6 Filtering using the IIR Chebyshev Type I filter
The following figures show the frequency response, pole-zero plot, phase response and the filtered signal compared to the original signal.
Chebyshev Poles & Zeros plot
Chebyshev Freq response of the filter 0
0.8
-20
0.6
-40
0.4
-60 Amplitude (dB)
1
Imag
0.2 0 -0.2 -0.4
-80 Chebyshev I Freq response Error-Max = 0.0344709 Error-RMS = 0.00113185 SNR = 43.9214 N = 10
-100 -120 -140
-0.6 -160 -0.8 -180
-1 -1
-0.5
0 Real Part
0.5
1
0
0.1
0.2
Chebyshev Impulse response of the filter
0.3
0.4
0.5 Freq x π
0.6
0.7
0.8
0.9
Chebyshev Phase response of the filter
0.3
4
0.25
3
0.2 2 Phase in Radians
Amplitude
0.15 0.1 0.05 0
1 0 -1
-0.05 -2 -0.1 -3
-0.15 -0.2
0
200
400
600
800 Time
1000
1200
1400
1600
-4
0
0.1
0.2
0.3
0.4
0.5 Freq x π
0.6
0.7
0.8
0.9
1
Figure 35: Chebyshev Filter Freq response (N=10)
38
The Filtered ECG signal Time domain plot 3
2.5
2.5
2
2
1.5
1.5 Amplitude (mV)
Amplitude (mV)
The Noisy ECG signal Time domain plot 3
1 0.5
1 0.5
0
0
-0.5
-0.5
-1
-1
-1.5
-1.5
0
0.2
0.4
0.6
0.8
1 1.2 Time (sec)
1.4
1.6
1.8
2
Filtering the 60 Hz Noize using Chebyshev Type 1 IIR Filter Error-Max = 0.0344709 Error-RMS = 0.00113185 SNR = 43.9214 Order = 10
0
0.2
The spectrum of the ECG Noisy Signal
0.4
0.6
0.8
1 1.2 Time (sec)
1.4
1.6
1.8
2
The spectrum of the Filtered Noise Signal
450
70
400
60
350 50 Amplitude
Amplitude
300 250 200
40 30
150 20 100 10
50 0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
0 -80
-60
-40
-20
0 20 Frequency (Hz)
40
60
80
Figure 36: Chebyshev Filtering result (N=10) ERROR Measurements for FIR Filtering using Chebyshev I Error-Max = 0.0344709 Error-RMS = 0.00113185 SNR = 43.9214 N = 10 5.1.3.7 Filtering using the FFT method
The following figures show the frequency response, of the used FFT filter kernel and the filtered signal compared to the original signal.
39
The Noisy ECG signal Time domain plot
The Filtered ECG signal Time domain plot
3
3
2.5
2.5 2 Amplitude (mV)
Amplitude (mV)
2 1.5 1 0.5 0 -0.5
1 0.5 0 -0.5
-1 -1.5
ERROR Measurements using FFT Method Error-Max = 0.0344709 Error-RMS = 0.00113185 SNR = 43.9214
1.5
-1 0
0.5
1 Time (sec)
1.5
-1.5
2
0
The spectrum of the ECG Noisy Signal
0.5
1 Time (sec)
1.5
2
The spectrum of the Filtered Noise Signal
500
70 60
400 Amplitude
Amplitude
50 300 200
40 30 20
100 10 0 -100
-50
0 Frequency (Hz)
50
100
0 -100
-50
0 Frequency (Hz)
50
100
Figure 37: FFT Filtering Results The FFT Filter Kernel Spectrum 1.5 ERROR Measurements using FFT Method Error-Max = 0.0344709 Error-RMS = 0.00113185 SNR = 43.9214
Amplitude
1
0.5
0 -60
-40
-20
0 Frequency (Hz)
20
40
60
Figure 38: FFT Filter Kernel
40
ERROR Measurements using FFT Method Error-Max = 0.0344709 Error-RMS = 0.00113185 SNR = 43.9214
5.1.3.8 Results comparison Table
The following table summarizes the filtering results from all developed filters: Method Four. Series (Chebyshev) Four. Series (Kaiser) Frequency Sampling Butterworth (IIR) Chebyshev I (IIR) FFT Filtering
Order 375
Max Error 0.0683814
RMS error 0.00167725
SNR 43.4112
229 229 25 10
0.0683609 0.0661603 0.0451924 0.0344709
0.0016757 0.00165881 0.00154873 0.00113185
43.45 43.5344 44.2513 43.9214
NA 0.0344709 0.00113185 Table 4 – All filters results comparison table
43.9214
5.1.3.9 Conclusions:
• • • •
The IIR filters supercedes the FIR filters in performance for much lower orders Chebyshev Type I ranked the best interms performance with the minimum order Kaiser window gave the best performance for the lowest order Filtering using FFT was very effective in achieving the desired output. Due to the need of calculating the FFT of the input signal and then the need to calculate the IFFT of the filtered spectrum such an approach might not be suitable for high speed real time signal processing applications.
• DSP IS FUN; I LOVE IT ☺ 6 References 1. 2. 3. 4. 5.
Prof. Sos Agaian DSP Lecture notes Digital Filters Designer’s Handbook, C. BrittonRorabaugh, 1st edition http://www-sigproc.eng.cam.ac.uk/~ad2/3f3/ Digital Signal Processing" by Ifeachor and Jervis, Addison-Wesley http://www.hsdal.ufl.edu/Projects/IntroDSP/
41
7 Matlab Code 7.1
Zero Order Bessel Function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Function to calculate the zero order % Bessel function for input x, N defines % the number of terms to use %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y = f_bessel(x,N); temp = 1; for k = 1:N temp = temp + (1/factorial(k))*(x/2)^k; end %return the result y = temp; 7.2 Generate the Kaiser Window function [order,win]=f_kaiser(f_sample,domega,delta) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% % Generate a Kaiser window for certin specs %----------------------------------------------------------------% Input: % fsample = sampling frequency % f_pass = 1st edge of the transition band % f_stop = 2nd edge of the transition band % ds = delta stop % dp = delta pass % % Output: % order = order of the FIR filter % win = the Kaiser window %---------------------------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define variables
42
fs = f_sample; %%%%%%%%%%%%%%%%%%%%%%%%%%% % Design the Kaiser window % low pass A = -20 * log(delta); disp('Attenuation (db) ='); disp(A); order = ceil(((A - 7.95 )/(2.285*domega)) + 2); %calculate the order odd if (rem(order,2)== 0) order = order + 1; end disp(order) %calculate beta if (A > 50) beta = .1102*(A - 8.7); elseif (A >= 21) beta = 0.5842*(A - 21)^.4 + .07886*(A - 21); else beta = 0 end alpha = order /2; %Compute the Kaiser window for n = 0:order-1 x= beta *(1 - (((n-alpha)/alpha)^2))^.5; win(n+1) =f_bessel(x,50)/f_bessel(beta,50); end %Return the results 7.3 Generate the Dolph-Chebyshev Window function [order,win]=f_cheby(f_sample,domega,delta) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% % Generate a Chebychev window for certin specs % Written by Hani Saleh for the DSP project 3 %----------------------------------------------------------------% Input: % fsample = sampling frequency % f_pass = 1st edge of the transition band % f_stop = 2nd edge of the transition band
43
% ds = delta stop % dp = delta pass % % Output: % order = order of the FIR filter % win = the Kaiser window %---------------------------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % define variables fs = f_sample; %%%%%%%%%%%%%%%%%%%%%%%%%%% % Design the cheby window % low pass Ts = 1/fs; % Ts = sampling interval. r = delta; A = -20 * log(r); disp('Attenuation (db) ='); disp(A); N = ceil((2.056*A - 16.4)/(2.85*domega))+1 if rem(N,2) == 0 N = N+1; end; disp('The length (order) of the filter is') order = N; disp(order) %Calculate the middle point M = (N-1)/2; % Calculate gama % gama = 10 ^(-1*A/20); gama = r; % % Calculating and plotting the window function. % if N > 1 beta = cosh((1/(N-1))*acosh(1/gama)); else beta = 0; end
44
%%% % for n = -1*M:1:M % win(n+M+1) = (1/N)*(1/gama); % for i = 1:M % x = beta*cos(i*pi/N); % if x >= -1 & x <= 1 % T = cos((N-1)*acos(x)); % else % T = cosh(acosh(x)); % end % win(n+M+1) = win(n+M+1) + (1/N)*2*T*cos(2*n*pi*i/N); % end % end %%%%% for k = 1:N win(k) = (1/N)*(1/gama); for i = 1:M x = beta*cos(i*pi/N); if x >= -1 & x <= 1 T = cos((N-1)*acos(x)); else T = cosh(acosh(x)); end win(k) = win(k) + (1/N)*2*T*cos(2*(k-1-(N-1)/2)*pi*i/N); end end scale=1/win((N-1)/2); for k = 1:N win(k) = scale*win(k); end %Use the matlab built-in window sice the above implimentation is not %working well %Return the results win = chebwin(N,A)'; 7.4
FIR Filtering using Fourier Series and Kaiser or Chebyshev windows (LP & HP) function d=fir_fseries(f_sample, f_pass,f_stop,dp,ds,type,wtype) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% % fir_fseries_lp_hp.m % FIR Filtering using the Fourier (Low Pass & High Pass) % Written by Hani Saleh for the DSP project 3
45
%----------------------------------------------------------------% Input: % fsample = sampling frequency % f_pass = Passband cutoff frequency for low pass and high pass % f_stop = Stopband cutoff frequencyfor low pass and high pass % ds = delta stop % dp = delta pass % type = type of filter % lo = low pass % ho = high pass % Window = "kais" for kaiser windwo or "cheb" % for doll-chebychev window %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % define variables fs = f_sample; delta = min(ds,dp); disp (delta); domega = 2*pi*(f_stop - f_pass)/fs; %%%%%%%%%%%%%%%%%%%%%%%%%%% % Design the Kaiser window if (wtype == 'cheb') [order,win] = f_cheby(fs,domega,delta); else [order,win] = f_kaiser(fs,domega,delta); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Design the filter disp(order) N = ((order-1)/2); L = 2*N+1; fc = (f_pass + ((f_stop - f_pass)/2)); wc = 2*pi*fc/fs; k=[];
if (type == 'lo') % low pass for n = -1*N:1: 1*N if(n ~=0 )
46
k(n+N+1) = (1/(1))*sin(n*wc)/(pi*n); else k(n+N+1)= wc/pi; end end else % high pass if (type == 'hi') for n = -1*N:1: 1*N if(n ~=0 ) k(n+N+1) = -1*sin(n*wc)/(pi*n); else k(n+N+1)= 1-wc/pi; end end end end b = k .* win; %%%%%%%%%%%%%%%%%%% % plot(win); % figure; [h1,w1]=freqz(k,1,[]); % figure; [h2,w2]=freqz(b,1,[]); %Normalize Freq to PI w1 = w1 /pi; w2 =w2/pi; %Return the results d= b; %impulse response [h3,t3] = impz(b,1); % zplane(b,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots the results figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t4,s4); % hold on; zplane(b,1); %plot(k,'LineWidth',2, 'Color','m');
47
% hold on; % plot(t1,abs(s1),'-g'); % axis([0 0.1 -1.5 1.5]); % % hold off; % % axis([0 0.1 -3 3]); grid; xlabel('Real Part') ylabel('Imag') title('Poles & Zeros plot'); %%%%%%%%%%%%%% subplot(2,2,2); %stem(t3,s3); % hold on; plot(b,'LineWidth',2, 'Color','k'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Time') ylabel('Amplitude') title('Filter Coeifficints after windowing'); %%%%%%%%%%%%%% subplot(2,2,3); plot(win,'LineWidth',2, 'Color','r'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Time') ylabel('Amplitude') title('The Selected winow'); %%%%%%%%%%%%%% subplot(2,2,4); % stem(t2,s2); % hold on; plot(t3,h3,'LineWidth',2, 'Color','c'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Time') ylabel('Amplitude') title('Impulse response after windowing');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
48
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots the results figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t2,s2); % hold on; plot(w1,20*log10(abs(h1)),'LineWidth',2, 'Color','c'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Freq') ylabel('Amplitude (dB)') title('Freq response befor windowing'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t4,s4); % hold on; plot(w2,20*log10(abs(h2)),'LineWidth',2, 'Color','m'); % hold on; % plot(t1,abs(s1),'-g'); % axis([0 0.1 -1.5 1.5]); % % hold off; % % axis([0 0.1 -3 3]); grid; xlabel('Freq') ylabel('Amplitude (dB)') title('Freq response after windowing'); %%%%%%%%%%%%%% %%%%%%%%%%%%%% subplot(2,2,3); plot(w1,angle(h1/pi),'LineWidth',2, 'Color','r'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Freq') ylabel('Phase in Radians') title('Phase response befor windowing'); %%%%%%%%%%%%%% subplot(2,2,4); % stem(t2,s2); % hold on; plot(w2,angle(h2/pi),'LineWidth',2, 'Color','c');
49
% hold off; % axis([0 0.1 -3 3]); grid; xlabel('Freq') ylabel('Phase in Radians') title('Phase response after windowing'); %%%%%%%%%%%%%% 7.5
FIR Filtering using Fourier Series and Kaiser or Chebyshev windows (BP & BS) function d=fir_fseries(f_sample, f_pass1,f_stop1,f_pass2,f_stop2,dp,ds,type,wtype) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% % filt_kern_fft.m % Filtering using FFT % Written by Hani Saleh for the DSP project 3 %----------------------------------------------------------------% Input: % data = data to be filtered % fsample = sampling frequency % f_pass1/f_stop1 = The first edge of the band pass and band reject % f_pass2/f_stop2 = The second edge of the band pass and band reject, % % ds = delta stop % dp = delta pass % wtype = "kais" for kaiser windwo or "cheb" % for doll-chebychev window %---------------------------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define variables fs = f_sample fc1 = (f_pass1 + (abs((f_stop1 - f_pass1)/2))) fc2 = (f_pass2 + (abs((f_stop2 - f_pass2)/2))) wc1 = 2*pi*fc1/fs wc2 = 2*pi*fc2/fs delta = min(dp,ds); domega1 = 2*pi*abs(f_stop1 - f_pass1)/fs; domega2 = 2*pi*abs(f_stop2 - f_pass2)/fs; domega = min(domega1,domega2); delta = min(dp,ds); %%%%%%%%%%%%%%%%%%%%%%%%%%%
50
% Design the Kaiser window if (wtype == 'cheb') [order,win] = f_cheby(fs,domega,delta); else [order,win] = f_kaiser(fs,domega,delta); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Design the filter disp('The length (order) of the filter is') disp(order) disp(win) disp(type) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Design the filter N = ((order-1)/2); L = 2*N+1; k=[];
if (type == 'bp') % bandpass for n = -1*N:1: 1*N if(n ~=0 ) k(n+N+1) = (sin(n*wc2)/(pi*n)) - (sin(n*wc1)/(pi*n)); else k(n+N+1)= wc2/pi - wc1/pi; end end elseif (type == 'bs') % bandstop for n = -1*N:1: 1*N if(n ~=0 ) k(n+N+1) = -1* sin(n*wc2)/(pi*n) + sin(n*wc1)/(pi*n); else k(n+N+1)= 1 - wc2/pi + wc1/pi; end end end disp(k) b = k .* win %%%%%%%%%%%%%%%%%%% % plot(win); % figure; [h1,w1]=freqz(k,1,[]); w1 = w1./pi; % figure;
51
freqz(k,1,[]); [h2,w2]=freqz(b,1,[]) w2 = w2./pi; %Return the results d= b; %impulse response [h3,t3] = impz(b,1); % zplane(b,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots the results figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t4,s4); % hold on; zplane(b,1); %plot(k,'LineWidth',2, 'Color','m'); % hold on; % plot(t1,abs(s1),'-g'); % axis([0 0.1 -1.5 1.5]); % % hold off; % % axis([0 0.1 -3 3]); grid; xlabel('Real Part') ylabel('Imag') title('Poles & Zeros plot'); %%%%%%%%%%%%%% subplot(2,2,2); %stem(t3,s3); % hold on; plot(b,'LineWidth',2, 'Color','k'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Time') ylabel('Amplitude') title('Filter Coeifficints after windowing'); %%%%%%%%%%%%%% subplot(2,2,3); plot(win,'LineWidth',2, 'Color','r'); % hold off; % axis([0 0.1 -3 3]); grid;
52
xlabel('Time') ylabel('Amplitude') title('The Selected winow'); %%%%%%%%%%%%%% subplot(2,2,4); % stem(t2,s2); % hold on; plot(t3,h3,'LineWidth',2, 'Color','c'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Time') ylabel('Amplitude') title('Impulse response after windowing');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots the results figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t2,s2); % hold on; plot(w1,20*log10(abs(h1)),'LineWidth',2, 'Color','c'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Freq x \pi') ylabel('Amplitude (dB)') title('Freq response befor windowing'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t4,s4); % hold on; plot(w2,20*log10(abs(h2)),'LineWidth',2, 'Color','m'); % hold on; % plot(t1,abs(s1),'-g'); % axis([0 0.1 -1.5 1.5]); % % hold off; % % axis([0 0.1 -3 3]);
53
grid; xlabel('Freq x \pi') ylabel('Amplitude (dB)') title('Freq response after windowing'); %%%%%%%%%%%%%% %%%%%%%%%%%%%% subplot(2,2,3); plot(w1,angle(h1/pi),'LineWidth',2, 'Color','r'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Freq x \pi') ylabel('Phase in Radians') title('Phase response befor windowing'); %%%%%%%%%%%%%% subplot(2,2,4); % stem(t2,s2); % hold on; plot(w2,angle(h2/pi),'LineWidth',2, 'Color','c'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Freq x \pi') ylabel('Phase in Radians') title('Phase response after windowing'); %%%%%%%%%%%%%% 7.6 Create White noise %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Function to generate a normalized white noise vector % N = length of vector %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function wn = f_wnoise(N); noise = randn(1,N); wn = noise ./ max(abs(noise)); % plot(abs(fft(wn))) % sound(wn,fs); % figure % hist(wn,t);
54
7.7 FIR Filtering using the Frequency Sampling Approach (LP, HP, BP & BS) function d=f_freq_sampling(f_sample, order,f_pass1,f_stop1,f_pass2,f_stop2,type) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% % filt_kern_fft.m % Filtering using Frequency Sampling method % Written by Hani Saleh for the DSP project 3 %----------------------------------------------------------------% Input: % fsample = sampling frequency % f_pass1/f_stop1 = The first edge of the band pass and band reject % f_pass2/f_stop2 = The second edge of the band pass and band reject, % % type = type of filter % lo = low pass % ho = high pass % bp = bandpass % br = band reject %---------------------------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define variables N = order; N2 = floor(order/2); fs = f_sample; fc1 = (f_pass1 + (abs((f_stop1 - f_pass1)/2))); fc2 = (f_pass2 + (abs((f_stop2 - f_pass2)/2))); wc1 = 2*pi*fc1/fs; wc2 = 2*pi*fc2/fs; domega1 = 2*pi*abs(f_stop1 - f_pass1)/fs; domega2 = 2*pi*abs(f_stop2 - f_pass2)/fs; domega = min(domega1,domega2); % delta = min(dp,ds);
fc1_num = floor(N * (fc1/fs)); fc2_num = floor(N * (fc2/fs)); f_pass1_num = floor(N * (f_pass1/fs)); f_stop1_num = floor(N * (f_stop1/fs)); f_pass2_num = floor(N * (f_pass2/fs)); f_stop2_num = floor(N * (f_stop2/fs));
55
kern(1:N) = 1; % low pass if (type == 'lo') kern(fc1_num:N2)=0; kern(N2+1:N-fc1_num)=0; end % high pass if (type == 'hi') kern(1:fc1_num)=0; kern(N-fc1_num:N)=0; end % bandpass if (type == 'bp') kern(1:fc1_num-1)=0; kern(fc2_num+1:N2-1)=0; kern(N2:N-fc2_num-1)=0; kern(N-fc1_num+1:N)=0; end % bandpass if (type == 'bs') kern(fc1_num:fc2_num)=0; kern(N-fc2_num:N-fc1_num)=0; end % Add samples into the transition band % kern(f_pass1_num) = 0.9; % kern(fc1_num) = .5; % kern(f_stop1_num) = 0.1; % % kern(f_pass2_num) = 0.9; % kern(fc2_num) = .5; % kern(f_stop2_num) = 0.1; % kern(N)=kern(2) if(rem(N,2)==1) M = (N-1)/2+1; for k = 1:N coef(k)=(1/N)*kern(1); for n = 2:M
56
coef(k)=coef(k)+(1/N)*2*kern(n)*cos((2*pi/N)*(k-M)*(n-1)); end end else M = N/2; for k = 1:N coef(k)=(1/N)*kern(1); for n = 2:M coef(k)=coef(k)+(1/N)*2*kern(n)*cos((2*pi/N)*(k-1-(N-1)/2)*(n-1)); end end end win = kern; %Return the results b= coef; %%%%%%%%%%%%%%%%%%% % plot(win); % figure; [h1,w1]=freqz(b,1,[]); w1 = w1./pi; %Return the results d= b; %impulse response [h3,t3] = impz(b,1); % zplane(b,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots the results figure; %%%%%%%%%%%%%% subplot(2,2,1); zplane(b,1); grid; xlabel('Real Part') ylabel('Imag') title('FIR (Freq-Sampling) Poles & Zeros plot'); %%%%%%%%%%%%%% subplot(2,2,2); plot(w1,20*log10(abs(h1)),'LineWidth',2, 'Color','c'); grid; xlabel('Freq x \pi') ylabel('Amplitude (dB)') title('FIR (Freq-Sampling) Freq response of the filter'); %%%%%%%%%%%%%%
57
subplot(2,2,4); plot(w1,angle(h1/pi),'LineWidth',2, 'Color','r'); grid; xlabel('Freq x \pi') ylabel('Phase in Radians') title('FIR (Freq-Sampling) Phase response of the filter'); %%%%%%%%%%%%% subplot(2,2,3); plot(t3,h3,'LineWidth',2, 'Color','c'); grid; xlabel('Time') ylabel('Amplitude') title('FIR (Freq-Sampling) Impulse response of the filter'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
7.8 Butterworth Normalized IIR Filter (Low pass) function [b,a] = f_butter(Wp, Ws, dp, ds) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%% % Design a normalized Butterworth low pass filter %Wp is the %pass band edge and Ws is the stop band edge. Frequencies are %normalized to [0,1], corresponding to the range [0,Fs/2]. %dp = pass band ripple in dbs %ds = stop band ripple in dbs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Normalized Sample period T = 2; % Wc = Wp + (Ws-Wp)/2; % warp the target frequencies according to the bilinear transform Ws = (2/T)*tan(pi*Ws./T); Wp = (2/T)*tan(pi*Wp./T); % compute minimum n which satisfies all band edge conditions % the factor 1/length(Wp) is an artificial correction for the
58
% band pass/stop case, which otherwise significantly overdesigns. qs = log(10^(ds/10) - 1); qp = log(10^(dp/10) - 1); N = ceil(max(0.5*(qs - qp)./log(Ws./Wp))/length(Wp)); disp('Order = '); disp(N); % compute -3dB cutoff given Wp, Rp and n Wc = exp(log(Wp) - qp/2/N); % unwarp the returned frequency Wc = atan(T/2*Wc)*T/pi; %%%%%%%%%%%%% % Calculate the coefficients % Prewarp to the band edges to s plane Wc = 2/T*tan(pi*Wc/T); % Generate splane poles for the prototype butterworth filter % source: Kuc C = 1; % default cutoff frequency pole = C*exp(1i*pi*(2*[1:N] + N - 1)/(2*N)); if mod(N,2) == 1, pole((N+1)/2) = -1; end % pure real value at exp(i*pi) zero = []; gain = C^N; % compensate for amplitude at s=0 gain = prod(-pole); % if n is even, the ripple starts low, but if n is odd the ripple % starts high. We must adjust the s=0 amplitude to compensate. if (rem(N,2)==0) gain = gain/10^(dp/20); end stop = 0; % S plane Low-pass to Low-passs frequency transform p = length(pole); z = length(zero); s_gain = gain * (C/Wc)^(z-p); s_pole = Wc * pole / C; s_zero = Wc * zero / C; % Use bilinear transform to convert S poles to the Z plane p = length(s_pole); z = length(s_zero);
59
Z_gain = real(s_gain * prod((2-s_zero*T)/T) / prod((2-s_pole*T)/T)); Z_pole = (2+s_pole*T)./(2-s_pole*T); if isempty(s_zero) Z_zero = -ones(size(Z_pole)); else Z_zero = [(2+s_zero*T)./(2-s_zero*T)]; Z_zero = postpad(Z_zero, p, -1); end % convert to the correct output form b = real(Z_gain*poly(Z_zero)); a = real(poly(Z_pole)); %%%%%%%%%%%%%%%%%%% % figure; [h1,w1]=freqz(b,a,[]); w1 = w1./pi; %Normalize Freq to PI w1 = w1 /pi; %impulse response [h3,t3] = impz(b,a); % zplane(b,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots the Freq Response figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t4,s4); % hold on; zplane(b,a); %plot(k,'LineWidth',2, 'Color','m'); % hold on; % plot(t1,abs(s1),'-g'); % axis([0 0.1 -1.5 1.5]); % % hold off; % % axis([0 0.1 -3 3]); grid; xlabel('Real Part') ylabel('Imag') title('Butterworth Poles & Zeros plot');
60
%%%%%%%%%%%%%% subplot(2,2,2); % stem(t2,s2); % hold on; plot(w1,20*log10(abs(h1)),'LineWidth',2, 'Color','c'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Freq x \pi') ylabel('Amplitude (dB)') title('Butterworth Freq response of the filter'); %%%%%%%%%%%%%%
%%%%%%%%%%%%%% subplot(2,2,4); plot(w1,angle(h1/pi),'LineWidth',2, 'Color','r'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Freq x \pi') ylabel('Phase in Radians') title('Butterworth Phase response of the filter'); %%%%%%%%%%%%%% subplot(2,2,3); % stem(t2,s2); % hold on; plot(t3,h3,'LineWidth',2, 'Color','c'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Time') ylabel('Amplitude') title('Butterworth Impulse response of the filter');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 7.9 Chebyshev Normalized IIR Filter (Low pass) function [b,a] = f_cheby1(Wp, Ws, dp, ds) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%% % Design a normalized Chebychev Type 1 low pass filter
61
%Wp is the %pass band edge and Ws is the stop band edge. Frequencies are %normalized to [0,1], corresponding to the range [0,Fs/2]. %dp = pass band ripple in dbs %ds = stop band ripple in dbs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Normalized Sample period T = 2; Wc = Wp + (Ws-Wp)/2; % warp the target frequencies according to the bilinear transform Ws = (2/T)*tan(pi*Ws./T); Wp = (2/T)*tan(pi*Wp./T); Wa = Ws/Wp; % compute minimum order which satisfies all band edge conditions stop_atten = 10^(abs(ds)/10); pass_atten = 10^(abs(dp)/10); N = ceil(acosh(sqrt((stop_atten-1)/(pass_atten-1)))/acosh(Wa)); disp(N) % Calculate the coefficients % Prewarp to the band edges to s plane Wc = 2/T*tan(pi*Wc/T); % Generate splane poles and zeros for the chebyshev type 1 filter C = 1; % default cutoff frequency epsilon = sqrt(10^(dp/10) - 1); v0 = asinh(1/epsilon)/N; pole = exp(1i*pi*[-(N-1):2:(N-1)]/(2*N)); pole = -sinh(v0)*real(pole) + 1i*cosh(v0)*imag(pole); zero = []; % compensate for amplitude at s=0 gain = prod(-pole); % if n is even, the ripple starts low, but if n is odd the ripple % starts high. We must adjust the s=0 amplitude to compensate. if (rem(N,2)==0) gain = gain/10^(dp/20); end
62
stop = 0; % S plane low pass to low pass frequency transform p = length(pole); z = length(zero); s_gain = gain * (C/Wc)^(z-p); s_pole = Wc * pole / C; s_zero = Wc * zero / C; % Use bilinear transform to convert the S poles to the z plane p = length(s_pole); z = length(s_zero); Z_gain = real(s_gain * prod((2-s_zero*T)/T) / prod((2-s_pole*T)/T)); Z_pole = (2+s_pole*T)./(2-s_pole*T); if isempty(s_zero) Z_zero = -ones(size(Z_pole)); else Z_zero = [(2+s_zero*T)./(2-s_zero*T)]; Z_zero = postpad(Z_zero, p, -1); end % convert to the correct output form b = real(Z_gain*poly(Z_zero)); a = real(poly(Z_pole)); % Plot the normalized freq responce % freqz(b,a,[]); %%%%%%%%%%%%%%%%%%% % figure; [h1,w1]=freqz(b,a,[]); w1 = w1./pi; %impulse response [h3,t3] = impz(b,a); % zplane(b,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots the results figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t4,s4); % hold on;
63
zplane(b,a); %plot(k,'LineWidth',2, 'Color','m'); % hold on; % plot(t1,abs(s1),'-g'); % axis([0 0.1 -1.5 1.5]); % % hold off; % % axis([0 0.1 -3 3]); grid; xlabel('Real Part') ylabel('Imag') title('Chebyshev Poles & Zeros plot'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t2,s2); % hold on; plot(w1,20*log10(abs(h1)),'LineWidth',2, 'Color','c'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Freq x \pi') ylabel('Amplitude (dB)') title('Chebyshev Freq response of the filter'); %%%%%%%%%%%%%%
%%%%%%%%%%%%%% subplot(2,2,4); plot(w1,angle(h1/pi),'LineWidth',2, 'Color','r'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Freq x \pi') ylabel('Phase in Radians') title('Chebyshev Phase response of the filter'); %%%%%%%%%%%%%% subplot(2,2,3); % stem(t2,s2); % hold on; plot(t3,h3,'LineWidth',2, 'Color','c'); % hold off; % axis([0 0.1 -3 3]); grid; xlabel('Time') ylabel('Amplitude')
64
title('Chebyshev Impulse response of the filter');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 7.10 Read the Biomedical signals %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This scripts reads the biomedical signals from the data files %%%%%%%%%%%%%%%%%%%%%% % Developed By Hani Saleh for Proj 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Read time vector in seconds (fsample = 125 Hz, Ts = 0.008 Seconds) fid = fopen('time.dat','r'); t = fscanf(fid,'%g'); status = fclose(fid);
% Read ECG signal_1 vector in mV fid = fopen('ecg.dat','r'); ecg_mv = fscanf(fid,'%g'); status = fclose(fid); % Read ECG signal_2 vector in mV fid = fopen('abp_mmhg.dat','r'); abp_mmhg = fscanf(fid,'%g'); status = fclose(fid); % Read ECG signal_2 vector in mV fid = fopen('mcl1_mv.dat','r'); mcl1_mv = fscanf(fid,'%g'); status = fclose(fid); % Read ECG signal_2 vector in mV fid = fopen('resp_mv.dat','r'); resp_mv = fscanf(fid,'%g'); status = fclose(fid); fs1 = 125; num_samples = 4096; % Remove the DC component mcl1_mv = mcl1_mv - mean(mcl1_mv); resp_mv = resp_mv - mean(resp_mv); abp_mmhg = abp_mmhg - mean(abp_mmhg); f_axis1 = (-1/2)*fs1:(fs1/num_samples):(1/2)*fs1-(fs1/num_samples);
65
%%%%%%%%%%%%%%%%%%%%%%%%%%%Plotting the signals figure subplot(2,1,1); plot(t,mcl1_mv,'LineWidth',2, 'Color','k'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The MCL signal Time domain plot'); %%%%%% subplot(2,1,2); plot(f_axis1,fftshift(abs(fft(mcl1_mv)))) xlabel('Freq. (Hz)') ylabel('Amplitude') title('The spectrum of the MCL Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,1,1); plot(t,resp_mv); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The RESP signal Time domain plot'); subplot(2,1,2); plot(f_axis1,fftshift(abs(fft(resp_mv))),'LineWidth',2, 'Color','r') xlabel('Freq. (Hz)') ylabel('Amplitude') title('The spectrum of the RESP Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,1,1); plot(t,abp_mmhg); xlabel('Time (sec)') ylabel('Amplitude (mmhg)') title('The ABP signal Time domain plot'); subplot(2,1,2); plot(f_axis1,fftshift(abs(fft(abp_mmhg))),'LineWidth',2, 'Color','g') xlabel('Freq. (Hz)') ylabel('Amplitude') title('The spectrum of the ABP Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,1,1); plot(t,ecg_mv); xlabel('Time (sec)') ylabel('Amplitude (mV)')
66
title('The ECG signal Time domain plot'); subplot(2,1,2); plot(f_axis1,fftshift(abs(fft(ecg_mv))),'k') xlabel('Freq. (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7.11 Filtering using FFT filter kernel function function d=filter_kern_fft(data,f_sample, f_cutoff1,f_cutoff2,type) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% % filt_kern_fft.m % Filtering using FFT % Written by Hani Saleh for the DSP project 2 %----------------------------------------------------------------% Input: % data = data to be filtered % fsample = sampling frequency % f_cutoff1 = cutoff frequency for low pass and high pass and it is % the first edge of the band pass and band reject % f_cutoff2 = The second edge of the band pass and band reject, % set to 0 for low pass and high pass % % type = type of filter % lo = low pass % ho = high pass % bp = bandpass % br = band reject %---------------------------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define variables N = length(data); fs = f_sample; fc1 = f_cutoff1; fc2 = f_cutoff2; N2 = int16(N/2); fc1_num = int16(N * (fc1/fs)); fc2_num = int16(N * (fc2/fs));
67
% low pass if (type == 'lo') data(fc1_num:N2)=0; data(N2+1:N-fc1_num)=0; end % high pass if (type == 'hi') data(1:fc1_num)=0; data(N-fc1_num:N)=0; end % bandpass if (type == 'bp') data(1:fc1_num-1)=0; data(fc2_num+1:N2-1)=0; data(N2:N-fc2_num-1)=0; data(N-fc1_num+1:N)=0; end % bandpass if (type == 'bs') data(fc1_num:fc2_num)=0; data(N-fc2_num:N-fc1_num)=0; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% %%%Plot the filter kernel kern(1:N) = 1; % low pass if (type == 'lo') kern(fc1_num:N2)=0; kern(N2+1:N-fc1_num)=0; end % high pass if (type == 'hi') kern(1:fc1_num)=0; kern(N-fc1_num:N)=0; end
68
% bandpass if (type == 'bp') kern(1:fc1_num-1)=0; kern(fc2_num+1:N2-1)=0; kern(N2:N-fc2_num-1)=0; kern(N-fc1_num+1:N)=0; end % bandpass if (type == 'bs') kern(fc1_num:fc2_num)=0; kern(N-fc2_num:N-fc1_num)=0; end fs1 = fs; num_samples = length(data); f_axis1 = (-1/2)*fs1:(fs1/num_samples):(1/2)*fs1-(fs1/num_samples); %%%%%%%%%%%%%%%%%%%%% figure; plot(f_axis1,fftshift(kern),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The FFT Filter Kernel Spectrum'); axis([(-fs1/2) (fs1/2) -.1 1.5]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%%%%% %Return the results d= data; 7.12 Radix-2 FFT used by the FFT Filter function d=fft_radix2_hsaleh(data_io,number_samples) % fft_radix2.m % Radix-2 FFT, Decimation in Frequency % % Custom Discrete Fourier transform.' % FFT(data_io) is the discrete Fourier transform (DFT) of vector data_io. % % For length N input vector data_io, the DFT is a length N vector data_io,' % with elements' % N' % data_io(k) = sum data_io(n)*edata_iop(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.' % n=1' % The inverse DFT (computed by IFFT) is given by' % N'
69
% data_io(n) = (1/N) sum data_io(k)*edata_iop( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.' % k=1' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 1/2/3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Detrmine if your input is real data vector if (mean(imag(data_io))== 0) real_data = 1; else real_data = 0; end no_add = 0; no_mult = 0; N = number_samples; M = log2(number_samples); % Compute the FFT using decimation in Freq for L = 1: M L; LE = 2^(M+1-L); no_add = no_add+2; no_mult = no_mult + (M-L); %no_mult = no_mult + 1;
LE2 = LE/2; no_mult = no_mult + 1; % Wiggle factors for 1st level U = 1.0; %Compute the Wigle factors S = cos(double(pi/LE2)) - j* sin(double(pi/LE2)); no_add = no_add+2; for J = 1:LE2 J; for I = J:LE:N I;
70
IP = I + LE2; T = data_io(I) + data_io(IP); data_io(IP) = (data_io(I) - data_io(IP)) * U; data_io(I) = T; no_add = no_add+3; no_mult = no_mult + 1; end U = U * S; no_mult = no_mult + 1; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Perform the bit reversal algorithm (re-order X(K) to get the correct results ND2 = N/2; NM1 = N -1; no_add = no_add+1; no_mult = no_mult + 1; %start with an index of 1 J = 1; for I=1:NM1 if I < J I; % Swap data_io(J) with data_io(I) T = data_io(J); data_io(J) = data_io(I); data_io(I) = T; end K = ND2; while (K < J) J = J -K; K = K /2; no_add = no_add+1; no_mult = no_mult + 1; end %Increment the index J = J + K;
71
no_add = no_add+3; end %For real data input make sure than data_io(N/2+m) = data_io(m) for m=2 ... N/2 % Due to the limitation of floating point accuracy there might be very % small diffrences, this part of the function remove these if (real_data==1) for INDEdata_io = 2:N/2 data_io(N+2-INDEdata_io)=data_io(INDEdata_io)'; end end
%Print adds 7 Mult no report cmplx=N*log(N); no_add; no_mult; % fprintf('-----------------\n'); % fprintf('N^2, NLog(N) , No of Adds , No of Mult\n'); % fprintf('%g , %g , %g , %g\n',N^2,cmplx,no_add,no_mult); % fprintf('-----------------\n'); %Return the results d= data_io; 7.13 Add 60 Hz noise to the ECG Signal and then use all of the developed filters to remove it. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% %%% This script test the filters on the biomedical signals % It injects 60 Hz nise and then removes it %%%%%% % Developed by Hani Saleh for DSP project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%% % Define variables fs1 = 125; num_samples = 4096; f_axis1 = (-1/2)*fs1:(fs1/num_samples):(1/2)*fs1-(fs1/num_samples);
%%%%%%%%%%%%%%%%%%% %%%% Read the read_biomed script; it defines the signals and create a time %%%% axis (t) disp('Call the read Biomed script');
72
read_biomed %%%%% disp('Create a normalized White noise vector'); wnoise = f_wnoise(num_samples); disp ('Plot the White noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,1,1); plot(t,wnoise,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The White Noise signal Time domain plot'); subplot(2,1,2); plot(f_axis1,fftshift(abs(fft(wnoise))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Create a 60 Hz Noize sigal freq1= 60; n60hz = cos(2*pi*freq1*t); %High sampling rate time vector for the time domain plot t_plot = 0:0.001:0.2; n60hz_plot = cos(2*pi*freq1*t_plot); disp ('Plot the 60 Hz noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,1,1); plot(t_plot,n60hz_plot,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The 60 Hz Noise signal Time domain plot'); subplot(2,1,2); plot(f_axis1,fftshift(abs(fft(n60hz))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% % Add noize to the ECG Signal ecg_n60 = ecg_mv + .2*n60hz; disp ('Plot the noisy ECG noise');
73
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,1,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The ECG signal with 60 Hz Noise Time domain plot'); % axis([0 .2 -1.5 1.5]); subplot(2,1,2); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG signal + Noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%% FIR by Fourier Series (Dolph-Cheby Window)%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% disp ('Filter the 60 Hz noise by FIR Filter (Fourier Series, Cheby)'); coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','cheb'); % Filter the signal ecg_n60_filt = filter(coef,1,ecg_n60); ecg_n60_orig = filter(coef,1,ecg_mv); disp ('Plot the Filtered noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot');
74
axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%% s4 error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error); ecg_snr = 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2); out_s1 = sprintf(' ERROR Measurements for FIR Filtering using Fourier Series (Cheby window)\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% %%%%% FIR by Fourier Series (Kaiser Window)%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% disp ('Filter the 60 Hz noise by FIR Filter (Fourier Series, kaiser)'); coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','kais'); % Filter the signal ecg_n60_filt = filter(coef,1,ecg_n60); ecg_n60_orig = filter(coef,1,ecg_mv); disp ('Plot the Filtered noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot');
75
axis([0 2 -1.5 3]); subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%% ECG error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error); ecg_snr = 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2); out_s1 = sprintf(' ERROR Measurements for FIR Filtering using Fourier Series (Kaiser Window)\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% %%%%% FIR by Frequency Sampling (Using Order of a Kaiser Windwo)%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% disp ('First: Filter the 60 Hz noise by FIR Frequency Sampling (Same as Kaiser order)'); %coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','kais');
76
order = 229; coef=f_freq_sampling(fs1,order,40,45,0,0,'lo'); % Filter the signal ecg_n60_filt = filter(coef,1,ecg_n60); ecg_n60_orig = filter(coef,1,ecg_mv); disp ('Plot the Filtered noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%% ECG error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error); ecg_snr = 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2); out_s1 = sprintf(' ERROR Measurements for FIR Filtering using Frequency Sampling (Same as Kaiser order\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr);
77
Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% %%%%% IIR Filtering - Butterworth %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% disp ('Filter the 60 Hz noise by IIR Butterworth'); %coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','kais'); [b,a] = f_butter((80/125),(90/125),7,69) % Filter the signal ecg_n60_filt = filter(b,a,ecg_n60); ecg_n60_orig = filter(b,a,ecg_mv); disp ('Plot the Filtered noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%
78
%%%%%%%%%%%%%%%%% ECG error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error); ecg_snr = 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2); out_s1 = sprintf(' ERROR Measurements for FIR Filtering using Butterworthr\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% %%%%% IIR Filtering - Chebyshev I %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% disp ('Filter the 60 Hz noise by IIR Chebyshev I'); %coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','kais'); [b,a] = f_cheby1((80/125),(90/125),7,69) % Filter the signal ecg_n60_filt = filter(b,a,ecg_n60); ecg_n60_orig = filter(b,a,ecg_mv); disp ('Plot the Filtered noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r');
79
xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%% ECG error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error); ecg_snr = 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2); out_s1 = sprintf(' ERROR Measurements for FIR Filtering using Chebyshev I\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% %%%%% Filtering Using FFT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% disp ('Filter the 60 Hz noise by The FFT method'); %coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','kais'); [b,a] = f_cheby1((80/125),(90/125),7,69) ecg_n60_fft = fft_radix2_hsaleh(ecg_n60,num_samples); ecg_n60_filt = filter_kern_fft(ecg_n60_fft,fs1,40,0,'lo'); % Filter the signal ecg_n60_filt = filter(b,a,ecg_n60); ecg_n60_orig = filter(b,a,ecg_mv); disp ('Plot the Filtered noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1);
80
plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%% ECG error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error); ecg_snr = 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2); out_s1 = sprintf(' ERROR Measurements using FFT Method\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
81
7.14 Add white noise to the ECG signal and use all of the developed filters to remove it %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% %%% This script test the filters on the biomedical signals % Using White noise % %%%%%% % Developed by Hani Saleh for DSP project 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%% % Define variables fs1 = 125; num_samples = 4096; f_axis1 = (-1/2)*fs1:(fs1/num_samples):(1/2)*fs1-(fs1/num_samples);
%%%%%%%%%%%%%%%%%%% %%%% Read the read_biomed script; it defines the signals and create a time %%%% axis (t) disp('Call the read Biomed script'); read_biomed %%%%% disp('Create a normalized White noise vector'); wnoise = f_wnoise(length(t))'; disp ('Plot the White noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,1,1); plot(t,wnoise,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The White Noise signal Time domain plot'); subplot(2,1,2); plot(f_axis1,fftshift(abs(fft(wnoise))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Create a 60 Hz Noize sigal freq1= 60; n60hz = cos(2*pi*freq1*t); %High sampling rate time vector for the time domain plot
82
t_plot = 0:0.001:0.2; n60hz_plot = cos(2*pi*freq1*t_plot); disp ('Plot the 60 Hz noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,1,1); plot(t_plot,n60hz_plot,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The 60 Hz Noise signal Time domain plot'); subplot(2,1,2); plot(f_axis1,fftshift(abs(fft(n60hz))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% % Add noize to the ECG Signal ecg_n60 = ecg_mv + .2*wnoise; disp ('Plot the noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,1,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The ECG signal with 60 Hz Noise Time domain plot'); % axis([0 .2 -1.5 1.5]); subplot(2,1,2); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG signal + Noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%% FIR by Fourier Series (Dolph-Cheby Window)%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% disp ('Filter the 60 Hz noise by FIR Filter (Fourier Series, Cheby)'); coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','cheb'); % Filter the signal ecg_n60_filt = filter(coef,1,ecg_n60); ecg_n60_orig = filter(coef,1,ecg_mv); disp ('Plot the Filtered noisy ECG noise');
83
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%% s4 error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error); ecg_snr = 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2); out_s1 = sprintf(' ERROR Measurements for FIR Filtering using Fourier Series (Cheby window)\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
84
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% %%%%% FIR by Fourier Series (Kaiser Window)%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% disp ('Filter the 60 Hz noise by FIR Filter (Fourier Series, kaiser)'); coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','kais'); % Filter the signal ecg_n60_filt = filter(coef,1,ecg_n60); ecg_n60_orig = filter(coef,1,ecg_mv); disp ('Plot the Filtered noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%% ECG error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error);
85
ecg_snr
= 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2);
out_s1 = sprintf(' ERROR Measurements for FIR Filtering using Fourier Series (Kaiser Window)\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% %%%%% FIR by Frequency Sampling (Using Order of a Kaiser Windwo)%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% disp ('First: Filter the 60 Hz noise by FIR Frequency Sampling (Same as Kaiser order)'); %coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','kais'); order = 229; coef=f_freq_sampling(fs1,order,40,45,0,0,'lo'); % Filter the signal ecg_n60_filt = filter(coef,1,ecg_n60); ecg_n60_orig = filter(coef,1,ecg_mv); disp ('Plot the Filtered noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot');
86
axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%% ECG error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error); ecg_snr = 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2); out_s1 = sprintf(' ERROR Measurements for FIR Filtering using Frequency Sampling (Same as Kaiser order\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% %%%%% IIR Filtering - Butterworth %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% disp ('Filter the 60 Hz noise by IIR Butterworth'); %coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','kais'); [b,a] = f_butter((80/125),(90/125),7,69) % Filter the signal ecg_n60_filt = filter(b,a,ecg_n60); ecg_n60_orig = filter(b,a,ecg_mv); disp ('Plot the Filtered noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot'); axis([0 2 -1.5 3]);
87
subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%% ECG error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error); ecg_snr = 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2); out_s1 = sprintf(' ERROR Measurements for FIR Filtering using Butterworthr\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% %%%%% IIR Filtering - Chebyshev I %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% disp ('Filter the 60 Hz noise by IIR Chebyshev I'); %coef=f_fir_fseries_lp_hp(fs1,40,45,.5,.001,'lo','kais'); [b,a] = f_cheby1((80/125),(90/125),7,69) % Filter the signal
88
ecg_n60_filt = filter(b,a,ecg_n60); ecg_n60_orig = filter(b,a,ecg_mv); disp ('Plot the Filtered noisy ECG noise'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure subplot(2,2,1); plot(t,ecg_n60,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Noisy ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,3); plot(f_axis1,fftshift(abs(fft(ecg_n60))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the ECG Noisy Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% subplot(2,2,2); plot(t,ecg_n60_filt,'LineWidth',2, 'Color','r'); xlabel('Time (sec)') ylabel('Amplitude (mV)') title('The Filtered ECG signal Time domain plot'); axis([0 2 -1.5 3]); subplot(2,2,4); plot(f_axis1,fftshift(abs(fft(ecg_n60_filt))),'LineWidth',2, 'Color','k') xlabel('Frequency (Hz)') ylabel('Amplitude') title('The spectrum of the Filtered Noise Signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%% %%%%%%%%%%%%%%%%% ECG error measurment ecg_error = (ecg_n60_orig) - (ecg_n60_filt); ecg_error_max = max (abs(ecg_error)); ecg_error_rms = std (ecg_error); ecg_snr = 10*log10 ( (std(ecg_n60_filt))^2/ (std(ecg_error))^2); out_s1 = sprintf(' ERROR Measurements for FIR Filtering using Chebyshev I\n'); out_s2 = sprintf('Error-Max = %g \n', ecg_error_max); out_s3 = sprintf('Error-RMS = %g \n', ecg_error_rms); out_s4 = sprintf('SNR = %g \n', ecg_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s2, out_s3, out_s4); display (Errror_measures);
89
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 7.15 Radix-2 IFFT used by the FFT Filter function d=ifft_radix2_hsaleh(Y,N) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ifft_radix2.m % % Custom Discrete Inverse Fourier transform.' % FFT(X) is the discrete Inverse Fourier transform (DFT) of vector X. % % The inverse DFT (computed by IFFT) is given by' % N' % x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.' % k=1' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 1/2/3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the complex conjugate of Input vector for INDEX = 1:N Y(INDEX) = Y(INDEX)'; end % Perform fft on the complex conjugate Y_FFT = fft_radix2_hsaleh(Y,N); % Scale the results of FFT to get the correct time domain values for INDEX = 1:N Y_FFT(INDEX) = (Y_FFT(INDEX)')/N; end % Retur the results d= Y_FFT; 7.16 Test the FFT Filter (LP, HP, BP & BS) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
90
% This m-file is used to test the Filters % % % Written by Hani Saleh for DSP project 2/3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Creat 4 time axisbases, at four sampling rates. fs1 = 2000; % Define the number of samples num_samples = 1024; %Create the time sampling vectors t1 = 0:(1/fs1):num_samples*(1/fs1)-(1/fs1); %High sampling rate time vector for the time domain plot t_plot = 0:0.001:0.1; % Four samples of a signal that is composed of two cosines functions freq1 = 25; freq2 = 200; freq3= 400; signal1 = cos(2*pi*freq1*t1); signal2 = cos(2*pi*freq2*t1); signal3 = cos(2*pi*freq3*t1); s1 = signal1 + signal2 + signal3; s_plot = cos(2*pi*freq1*t_plot) + 2*cos(2*pi*freq2*t_plot);
S1 = fft_radix2_hsaleh(s1,num_samples); f_axis1 = (-1/2)*fs1:(fs1/num_samples):(1/2)*fs1-(fs1/num_samples); %filter the data S1_filt = filter_kern_fft(S1,fs1,100,0,'lo'); S2_filt = filter_kern_fft(S1,fs1,100,0,'hi'); S3_filt = filter_kern_fft(S1,fs1,100,300,'bp'); S4_filt = filter_kern_fft(S1,fs1,100,300,'bs'); % reconstruct the time domain data s1_reconst = ifft_radix2_hsaleh(S1_filt,num_samples); s2_reconst = ifft_radix2_hsaleh(S2_filt,num_samples); s3_reconst = ifft_radix2_hsaleh(S3_filt,num_samples); s4_reconst = ifft_radix2_hsaleh(S4_filt,num_samples); s1_orig = signal1;
91
s2_orig = signal2+signal3; s3_orig = signal2; s4_orig = signal1+signal3; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% s1 error measurment s1_error = s1_orig - s1_reconst; s1_error_max = max (abs(s1_error)); s1_error_rms = std (s1_error); s1_snr = 10*log10 ( (std(s1_orig))^2/ (std(s1_error))^2); %%%%%%%%%%%%%%%%% s2 error measurment s2_error = s2_orig - s2_reconst; s2_error_max = max (abs(s2_error)); s2_error_rms = std (s2_error); s2_snr = 10*log10 ( (std(s2_orig))^2/ (std(s2_error))^2); %%%%%%%%%%%%%%%%% s3 error measurment s3_error = s3_orig - s3_reconst; s3_error_max = max (abs(s3_error)); s3_error_rms = std (s3_error); s3_snr = 10*log10 ( (std(s3_orig))^2/ (std(s3_error))^2); %%%%%%%%%%%%%%%%% s4 error measurment s4_error = s4_orig - s4_reconst; s4_error_max = max (abs(s4_error)); s4_error_rms = std (s4_error); s4_snr = 10*log10 ( (std(s4_orig))^2/ (std(s4_error))^2); out_s1 = sprintf(' ERROR Measurements, Error = Original-Signal - IFFT-of-Radix2Form-FFT \n'); out_s1a = sprintf(' Low \t\t High \t\t BP \t\t BS\n'); out_s2 = sprintf('Error-Max = %g \t %g \t %g \t %g \t \n', s1_error_max,s2_error_max,s3_error_max,s4_error_max); out_s3 = sprintf('Error-RMS = %g \t %g \t %g \t %g \t \n', s1_error_rms,s2_error_rms,s3_error_rms,s4_error_rms); out_s4 = sprintf('SNR = %g \t %g \t %g \t %g \t \n', s1_snr,s2_snr,s3_snr,s4_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s1a, out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots for the low pass figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t1,s1,'-*');
92
% hold on; plot(t1,s1,'-g'); % hold off; axis([0 0.1 -3 3]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,2,3); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input Signal Spectrum FFT'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t4,s4); % hold on; plot(t1,abs(s1_reconst),'-r'); hold on; plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('Low pass Reconstructed signal (green- Original, Red - Reconst'); %%%%%%%%%%%%%% subplot(2,2,4); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S1_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Low pass Filtered signal fft');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots for the high pass figure; %%%%%%%%%%%%%% subplot(2,2,1);
93
% stem(t1,s1,'-*'); % hold on; plot(t1,s1,'-g'); % hold off; axis([0 0.1 -3 3]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,2,3); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t4,s4); % hold on; plot(t1,abs(s2_reconst),'-r'); hold on; plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('High pass Reconstructed signal'); %%%%%%%%%%%%%% subplot(2,2,4); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S2_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('High pass Filtered signal fft'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots for the BandPass pass figure; %%%%%%%%%%%%%% subplot(2,2,1);
94
% stem(t1,s1,'-*'); % hold on; plot(t1,s1,'-g'); % hold off; axis([0 0.1 -3 3]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,2,3); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t4,s4); % hold on; plot(t1,abs(s3_reconst),'-r'); hold on; plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('Band pass Reconstructed signal'); %%%%%%%%%%%%%% subplot(2,2,4); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S3_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Band pass Filtered signal fft'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots for the Band Stop pass figure; %%%%%%%%%%%%%% subplot(2,2,1);
95
% stem(t1,s1,'-*'); % hold on; plot(t1,s1,'-g'); % hold off; axis([0 0.1 -3 3]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,2,3); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t4,s4); % hold on; plot(t1,abs(s4_reconst),'-r'); hold on; plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('Band Stop Reconstructed signal'); %%%%%%%%%%%%%% subplot(2,2,4); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S4_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Band Stop Filtered signal fft'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots of all the samples combined figure; %%%%%%%%%%%%%% subplot(2,4,1);
96
% stem(t1,s1,'-*'); % hold on; plot(t1,signal1,'-g'); % hold off; axis([0 0.1 -1.5 1.5]); grid; title('Signal1, Freq = 25 Hz'); %%%%%%%%%%%%%% subplot(2,4,4); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% subplot(2,4,2); % stem(t4,s4); % hold on; plot(t1,signal2,'-r'); % hold on; % plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('Signal2, Freq = 200 Hz'); %%%%%%%%%%%%%% subplot(2,4,3); % stem(t4,s4); % hold on; plot(t1,signal3,'-r'); % hold on; % plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('Signal3, Freq = 400 Hz'); %%%%%%%%%%%%%% subplot(2,4,5); %stem(t3,s3); % hold on;
97
plot(f_axis1,fftshift(abs(S1_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Low pass Filtered signal fft'); %%%%%%%%%%%%%% subplot(2,4,6); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S2_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('High pass Filtered signal fft'); %%%%%%%%%%%%%% subplot(2,4,7); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S3_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Band pass Filtered signal fft'); %%%%%%%%%%%%%% subplot(2,4,8); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S4_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Band reject Filtered signal fft'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 7.17 Test the developed Filters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % This Script is used to test the Filters while they were under development % % % Written by Hani Saleh for DSP project 3
98
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% % Creat time axis bases, at one sampling rates. fs1 = 1000; % Define the number of samples num_samples = 1024; %Create the time sampling vectors t1 = 0:(1/fs1):num_samples*(1/fs1)-(1/fs1); %High sampling rate time vector for the time domain plot t_plot = 0:0.001:0.1; % Four samples of a signal that is composed of two cosines functions freq1 = 50; freq2 = 200; freq3= 300; signal = cos(2*pi*freq1*t1); noize = cos(2*pi*freq2*t1) + cos(2*pi*freq3*t1); s1 = signal + .25*noize; s1 = s1+ f_wnoise(length(s1)); s_plot = cos(2*pi*freq1*t_plot) + 2*cos(2*pi*freq2*t_plot);
S1 = fft_radix2_hsaleh(s1,num_samples); f_axis1 = (-1/2)*fs1:(fs1/num_samples):(1/2)*fs1-(fs1/num_samples); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% %filter the data ; Enable the filter you like to test below: %%%%%%%%%%%%%%%%%%%%%%%%%% % coef=f_fir_fseries_lp_hp(fs1,100,150,.1,.01,'lo','cheb'); % coef=f_fir_fseries_lp_hp(fs1,100,150,.1,.01,'hi','kais'); coef=f_fir_fseries_bp_bs(fs1,100,150,250,300,.1,.01,'bs','kais'); % order = 121; % coef=f_freq_sampling(fs1,order,100,150,250,300,'lo'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%% %%%%%Use this to perform freq transformation for the coefficients % for k = 1:length(coef) % shift(k) =sin(2*pi*k/4); % end
99
% coef = shift(k) * coef; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%% coef_padded = coef; if(length(coef) < num_samples) coef_padded((length(coef)+1):num_samples) = 0; end coef_fft = fft_radix2_hsaleh(coef_padded,num_samples); s1_filt = filter(coef,1,s1); s1_filt_conv = coef_fft .* S1; S_filt = fft_radix2_hsaleh(s1_filt,num_samples); figure plot(f_axis1,abs(fftshift(S_filt))) figure plot(f_axis1,abs(fftshift(S1))) figure plot(f_axis1,abs(fftshift(s1_filt_conv)))
100
