This assignment is about FIR filter design, implementation, and applications.

On xSiTe, there are 10 sets of ECG signals. Use the data file according to the last digit in your student number, for example, if your student number ends with 0, you should use “ecg_0.dat” . In MATLAB, use load ‘ecg_0.dat’ to load the data. You will get a matrix called

ecg_0 having 2 columns, and the second column is the ECG data, you may use y = ecg_0(:, 2) to get a data vector.

Step 1: Display the original ECG in mV against milliseconds. Show the whole recording and one heartbeat to identify the so called PQRST waves (see Figure 1). The amplifier before the ADC has a gain of 200; the A/D converter has a 24-bit resolution and has an input range from −1.325V to +1.325V, with the sampling rate of 1KHz (see Figure 2). Develop your own formula to convert the data from integers into ECG in mV. Do NOT use any pre-defined MATLAB functions to convert it. Plot also the frequency spectrum (no mirror) of the ECG. Identify what is signal and what is noise, in both the time and frequency domains. In particular, identify the unwanted 50Hz pollution and the DC drift

of the  ECG.  Label your  figures  properly  (time  in  milliseconds,  frequency  in  Hz,  and amplitude in mV).

Figure 1. A normal ECG with P,Q,R,S,T waves

Figure 2. ECG data acquisition system

(15%)

Step 2: Write a MATLAB function which implements an FIR filter (refer to our lecture notes of Wk9, the 4 steps, see Figure  3) with the impulse response h(n) in the form value myFIR(value, hn), where both the value argument and return value are scalars and NOT vectors so that it could be used in a real-time system. You may define this function as, for example

function output = myFIR(input, hn)

where both output and input must be scalars, and hn is the impulse response h(n) of the  FIR  filter  (where  hn  is  a  vector).  No  high  level  MATLAB  commands  (eg  “filter”, “conv”, . . . ) are allowed. You can test it with the Moving Average FIR filter and make the

input with a pulse, your output should be the impulse response as shown in the lecture of Week 9.

Figure 3. FIR filter implementation

(10%)

Step  3: Design an  FIR High Pass filter which eliminates the  DC drift of the  ECG but leaves  the  ECG  intact.  Define  the  ideal  high  pass  frequency  response  and  perform analytically   an   inverse   Fourier   transform   to   obtain   the   FIR   coefficients.   Use

appropriate window functions. What is the fundamental frequency of an ECG? What is the minimum number of taps?

Use the myFIR function that you developed in Step 2 to filter the raw ECG signal. No high level MATLAB commands (for example, “fir1”, “filter”, “conv”) are allowed. If you are not sure whether a particular function is allowed or not, please speak to me before using it. Your main program should contain a loop which feeds the samples one-by-one into the filter as if the signal is acquired in real-time.  After you have filtered the ECG, check that all individual peaks are still intact (see Figure 1) and that it is not distorted.

(25%)

Step 4: Design a single FIR filter which eliminates both the DC drift and the 50Hz hum at the same time. Design an ideal filter response which has zero transmission for DC and 50Hz but leaves all other frequencies intact. Then, do the inverse Fourier transform numerically in MATLAB by defining your ideal frequency response as a vector and then

transforming it into the time domain to obtain the impulse response.

Use the myFIR function that you developed above to filter the input signal. Again, no high level MATLAB commands (eg “fir1”, “filter”, “conv”, . . . ) are allowed. Your main program should contain a loop which feeds the samples one-by-one into the filter as if the signal is acquired in real-time. What is the minimum number of taps? After you have filtered the ECG check that all individual peaks are still intact (see Figure 1) and that it is not distorted, but the DC drift and 50Hz noise has been removed.  Show zoomed-in ECG traces of one heartbeat so that it is possible to identify the different parts of a single heartbeat and that it’s possible to check if it’s still intact.

Note: In Step 3 and Step 4, when you design your filters, plot the actual amplitude response of your filters, and check the performance criteria (passband/stopband ripple,   transition   width,   and   stopband   rejection).   You   should   tweak   the parameters  like  the  cut-off  frequencies,  the  number  of  taps,  and  the  window functions to make a trade-off between the performance criteria, and maximize the filtering effect. In your report, show your final choice of parameters, the resulted actual filter amplitude response, and the output ECG signal. You should also show the zoomed-in ECG traces of one heartbeat so that it is possible to check if it’s still intact. Explain and comment your choices.

(25%)

Step 5: ECG Heartrate Detection

The task is to detect the momentary heart rate r(t) over a period of time. It is NOT the average  heart  rate  but  the  frequency  derived  from  the  times  between  adjacent heartbeats.

Design a detector which detects the individual heartbeats in the ECG, and implement this filter as a matched filter. The matched filter should detect the R-peaks of the ECG, and  the  Template  of  the Matched Filter  should  be  calculated  analytically from  a mathematical formula. Think which analytical functions look like an R-peak and then calculate the template by  sampling  such  a  function that it  matches the  R-peaks.  For example: Mexican hats, differences of Gaussians, derivatives of Gaussians, etc. Test at least two different Templates based on two functions, decide which one is best suitable and provide a reason based on the outputs of the matched filters.

Implement a program in MATLAB which uses the output of the matched filter on the filtered  ECG  (from Step 4) to  calculate the momentary heart  rate  r(t)  over time by measuring the intervals between the detected heartbeats over the whole period of the ECG. Detect the heartbeats by employing a threshold. Square the output of the detector to improve the signal to noise ratio. Use any heuristics to weed out wrong detection, for example a heartbeat is usually below 200bpm and above 30bpm. Show the momentary heart rate r(t).

(25%)