Menu:
Updates:
Spt, 2012:
Computing in cardiology conference in Krakow
Mar 7, 2012:
King's fund conference in London
ECG and Wiener filter
Last updated: 09/02/2012
1) Introduction
In this tutorial you will learn about non-causal frequency-domain Wiener filtering for its application to ECG. The Wiener filter is optimal in the minimum mean square error (MMSE) sense for stationary signals. It can be applied in a causal sense in the time domain or as a non causal filter in the frequency domain. In this tutorial we will use it in the frequency domain.
Prerequisite: Being familiar with Fourier transformation and basic digital signal processing. If you are not familiar with it, I advise you this book by Dr Steven W. Smith.
Going further than this tutorial: For some theory about the Wiener filter: Chapter 6, Advanced Digital Signal Processing and Noise reduction by Saeed V. Vaseghi.
2) Initial samples
Initial samples for this tutorial are: two 10s clean ECGs (called idealECG_1 and idealECG_2) to which I have added 50Hz noise and white Gaussian noise with $SNR = 12db$ (we will call the resulting samples observation_1 and observation_2). The two 10s samples have been extracted from the same ECG record but at different times. Figure 1 shows the idealECG_1 and observation_1 samples. Keep in mind that you are working on a case study! If you want to try your algorithms on real data, I recommend you use some of the Physionet Databases. The reason we need two samples of ideal and observation ECGs is because we need to first build a model to be used by the Wiener filter - on set-1 (idealECG_1, observation_1)- and then test the model on set-2 (idealECG_2, observation_2). You should now download set-1.
[ECG samples]
Figure 1: A few seconds of the idealECG sample and observation sample (set-1). Note the axis: seconds for the x-axis and mV for the y-axis. The 50Hz noise on the observation sample is very visible.
| Another thing that we should control is the type of error distribution that we are dealing with. Indeed, the least square error criterion is optimal for Gaussian distributed signals [2]. So in our case we should control that the error e (e=idealECG-observation) is normally distributed. Figure 2 shows the histogram we get which looks enough like Gaussian distributed (for our case study!). |
Figure 2: Error distribution obtained from set-1.
|
2) The Wiener filter
The non-causal frequency-domain Wiener filter is designed for its application in the frequency domain and is based on the relative amounts of signal and noise present at each frequency. The transfer function of the optimal (non-causal) Wiener filter is defined by: $$ H[f]={{S[f]^2}/{S[f]^2+N[f]^2}}$$where $H[f]$ is the frequency response of the filter, $N[f]$ the frequency spectra of the noise and $S[f]$ the frequency spectra of the signal.
This can also be written as: $$H[f]={{SNR(f)}/{SNR(f)+1}}$$
where $SNR(f) = {S[f]^2}$/${N[f]^2}$ is the signal to noise ratio expressed in terms of the power-spectral ratio. Also for additive noise, the Wiener filter attenuates each frequency component in proportion to an estimate of the signal to noise ratio [2].
The Wiener filter is optimal in the sense that it maximizes the ratio of the signal power to the noise power [3]. It is important to note that this filter is not an adaptive filter (it does not self-adjust its transfer function according to some optimization algorithm controlled by the error signal). Indeed, the filter assumes that the inputs are stationary. For a good mathematical derivation of the Wiener transfer function, the reader is advise to read Chapter 6, Advanced Digital Signal Processing and Noise reduction by Saeed V. Vaseghi.
3) Building the model
We first need to build a model of the noise and ECG spectrum. From this model we will compute $ H[f] $ as previously defined and use it on another 10s noisy sample (set-2) - see part 4). You should now download the Matlab functions: main_Wiener_v1.m, wienerFilter.m and pwelch_Wiener.m.[ECG Code]
Figure 3 shows the result of filtering observation_1 using the Wiener filter. As you can see, the filtering is relatively good which is expected as we are building the ECG model on observation_1 and idealECG_1 which means that we have perfect knowledge of the ECG model spectral density (ECG model = idealECG_1) and of the noise spectral density (noise = idealECG_1-observation_1). Thus this experiment shows the best type of result we can expect by using the Wiener filter on our sample. Figure 4 shows the spectral density of the ECG model (idealECG), of the observation (i.e ECG to be processed) and the result of filtering the observation using the Wiener filter (filtered ECG). The black dashed line represents the frequency content of the noise obtained by taking the Fast Fourier Transform of the noise (=idealECG_1-observation_1=50Hz noise + Gaussian white noise). As you can see the spectra of the filtered ECG (red) is close to the one of the model ECG (blue) at least until 100Hz. From the filter transfer function $H$: if the $SNR$ is very high then $H≃1$ and if the $SNR$ is low then $H≃SNR$ which can help interpreting the behaviour of the graph. In particular the 50hz noise is well filtered; obtaining these result with a Notch filter for power line noise filtering could have been possible with a very high quality factor (Q). The advantage of using the Wiener filter (for this particular case of 50Hz noise) is that there is no need to 'empirically' try/adjust the Q factor; given an ECG and noise model the filter takes care of optimizing by how much the 50Hz frequency content should be attenuated.
Figure 3: idealECG and filtered observation (set-1) using the Wiener filter.
Figure 4: Spectral density of the idealECG (model), the observation (ECG to process) and the result of the filtering (filtered ECG). The spectra of the 50Hz noise (noise model) is not represented for clarity of the plot but it would correspond to a spike at 50Hz.
4) Applying the model
The question now is: how will the filter perform on another 10s ECG chunk (called observation_2) that has not been used to define the model? For this part you will need to download set-2 (idealECG_2, observation_2) and filter observation_2 using the predefined Wiener model in step 3). Then compare the filtered result (filtered observation) with idealECG_2.[ECG Code]
As you can see on Figure 5 the result is still acceptable.
Figure 5: idealECG and filtered observation (set-2) using the Wiener filter.
5) What if there is new type of noise?
What would happen if we have to deal with a type of noise which has not been taken into account in the noise spectrum model? Take the imaginary case where you are teleported from the UK to the US and your ECG spectrum model has been initialized in UK but you continue to be monitored in the US. Let's assume that set-a corresponds to the record in UK on which the ECG model is built. Now that we are in US we are likely to encounter 60Hz power interference noise. What will happen?As you can see on Figure 6 the 60Hz noise has not been filtered at all which is expected as we did not take it into account when building the model. Although this is a case study, more realistic cases such as electrode motion are very unpredictable and might not be taken into account when first creating the model. This is why in practice it might be more convenient to build an approximate model of the ECG and noise spectrum and use this model in a more generic way.
Figure 6: idealECG and filtered observation (containing 60Hz noise) using the Wiener filter. The Wiener filter model has been built on the idealECG_1 and observation_1 (containing 50Hz noise).
6) Conclusion on the Wiener filter
It is important to note that the Wiener filter is completely determined by the characteristics of the problem. Although it is optimal in the sense that it maximizes the signal to noise ratio it does not mean that it will make the electro-cardiogram easier to interpret by the physician [3].If you are considering using the Wiener filter for an online application, then you have to look into the FIR Wiener filters formulation in the time domain.
Let me end this tutorial by encouraging you playing around with the code. I hope you found it useful and do not hesitate to send me an email if so (it always nice to know!) or if you disagree with some of the above.
References
[1] Paul Kligfield et al. Recommendations for the Standardization and Interpretation of the Electrocardiogram. Part I: The Electrocardiogram and Its Technology. Journal of the American College of Cardiology. Vol. 49, No. 10, 2007[2] Saeed V. Vaseghi. Chapter 6, Advanced Digital Signal Processing and Noise Reduction, Second Edition. Saeed V. Vaseghi Copyright © 2000 John Wiley & Sons Ltd.
[3] Steven W. Smith, Ph.D. The Scientist and Engineer's Guide to Digital Signal Processing. http://www.dspguide.com/pdfbook.htm
[4] Gari D. Clifford, Francisco Azuaje and Patrik McSharry. Advanced Methods and Tools for ECG Data Analysis. Artech House. 2006