Regular ad hoc filters don’t guarantee optimal signal filtering as there is no any criteria that evaluates filter characteristics. Usually filter parameters are calculated empirically and filtering is done by best results.
In order to avoid such shortage there are optimal filters used where parameters are optimized by some criteria. The main idea of optimal filtering is to give bigger weight coefficients to signal spectra parts where signal noise has less power and true signal spectral components has bigger power.
Lets project simple Butterworth filter that will be used as comparative filter to optimal Wiener DSP filter.
Butterworth filter transfer characteristics:
Where N – indicates filter Tap number. I will skip description of butterworth filter for now as main idea is constructing optimal wiener filter. Butterworth filter characteristics are pretty plain:
Butterworth characteristics with bandpass of 70Hz
Main disadvantage of Butterworth filter is that signal is distorted on filter output. If you want minimal signal distortions it is better to use optimal Wiener filter. Filter chart looks as follows:
As you can see to make this filter functional wee need additional conditions like signal model reference. Filter coefficients are calculated using MMSE. root mean square error. Lets see how it works.
Filter input is ideal signal plus noise:
and filter output depends on filter weights of coefficients wk
Further wee need to optimize filter coefficients wk in order to meet minimal signal error e(n):
As we can see minimal root mean square error depends on ideal signal model construction.
In frequency domain optimal Wiener filter coefficients are calculated by using Wiener – Hoph equation:
Where Sη(w) – noise power spectral density;
Sd(w) – signal model power spectral density.
Where Rη(m)- noise self correlation function.
Lets see what results wee get. First of all we have made Butterworth low frequency filter. Signal sampling frequency were used – 1000Hz. We have made couple examples when filter tap numbers are 4 and 8 and cutoff frequencies are 30 and 70Hz:
It is really hard to assemble good parameters for Butterwoth filter because there is no good criteria to value filter accuracy. So following step is optimal filter construction.
First of all we need signal model reference:
According to previous formulas we calculate filter coefficients depending on signal and signal model spectral densities:
a) Signal noise spectral density; b) real signal spectral density; c) signal model spectral density ; d) Filter coefficients dependency on frequency.
When we have Wiener-Hoph equation with all required coefficients we can filter ECG signal. First of all filtering is done in frequency domain. Lets calculate signal Fourier transformation, then filtering is nothing more than multiplying signal Fourier transformation by filter coefficients. Results are in following picture:
a) real signalo Fourier transformation; b) filtered signal Fourier transformation; c) real signal; d) filtered signal
Of course you can filter signal in time domain. For this you have to make reverse Fourier transformation of filter coefficients. If wee take 30 of filter coefficients we get:
Signal filtered in time domain using 30 of filter coefficients.
Now we can compare filtered signals by using Butterworth and Wiener filters. One of comparative methods can be calculation of group delay. In following results we can see that using Butterworth filter group delay is biggest at 70Hz while Wiener filters resulting group delay oscillates around 0.
a) Butterworth filter group delay; b) Wiener filter group delay
According to results it is really hard to say which filtering method is better. We can see, that using tap 8 and 30Hz Butterworth filter it gives more clear Q segment while using Wiener filter this part of signal is hardly noticeable. Of course signal model is quite rough, this results in less precise filter coefficients.
%%%%%%%%Matlabsource%%%%%%%%%%%%%%%%%%%%%
%Optimal ECG signal filtering using Wiener filter”
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
%1. Projecting Butterworth filter
Fp=70;
Fd=1000;
eil=8;
%Cut-Off
ecg = load(‘ecg_hfn.dat’);
ecg=ecg/max(ecg(2100:2800));
figure
[B,A]=butter(eil/2,Fp/(Fd/2));
N=length(ecg);
timeall=N/Fd;
tecg=0:1/Fd:timeall;
plot(tecg(1:700),ecg(2101:2800));
grid on;
title(‘measured signal and filtered with Butterworth filter’);
xlabel(‘t, s’);
ylabel(‘Amplitude, mV’);
hold on;
fecg1 = filter(B,A,ecg);
for i=1:N
ecg1rev(i)=fecg1(N+1-i);
end
fecg2rev = filter(B,A,ecg1rev);
for i=1:N
fecg(i)=fecg2rev(N+1-i);
end
poslinkis=eil*3+1
plot(tecg(1:700),fecg(2101:2800),’red’);
xlabel(‘t, s’);
ylabel(‘Amplitude, mV’);
legend(‘Real signal and filtered signal’);
hold off;
%Butterworth characteristics
figure;
freqz(B,A,700,Fd);
figure
grpdelay(B,A,700,Fd);
%Group delay
load mas1.mat xd yd;
clear des_x; clear des_y;
des_x=[]; des_y=[];
%Model Interpolation
des_x=xd-xd(1);
tasku=round(xd(end)-xd(1));
des_y=max(ecg(2100:2800))*yd;
xi=0:tasku-1;
% Aproximation with cubic spine
yi=interp1(des_x,des_y,xi,’cubic’);
figure
plot(des_x,des_y,’o',xi,yi);
grid on;
xlabel(‘t, s’);
ylabel(‘Amplitude, mV’);
des_x = des_x(1:13);
des_y = des_y(1:13);
figure(3);
plot(des_x,des_y,’*',des_x,des_y);
hold on;
grid on;
xlabel(‘Time, ms’);
ylabel(‘Amplitude’);
axis tight;
hold off;
ans = 1;
end
%
load trm.mat
x1=trm;
L1=x1(2)-x1(1);
L2=x1(4)-x1(3);
L3=x1(6)-x1(5);
L4=x1(8)-x1(7);
L=[L1; L2; L3; L4];
Lmin=min(L);
ecgtr=[ecg(x1(2):x1(2)+Lmin), ecg(x1(3):x1(3)+Lmin), ...
ecg((x1(6)-Lmin):x1(6)),ecg(x1(7):x1(7)+Lmin)];
Tr=mean(ecgtr,2);
Trvid=mean(Tr,1);
Tr=Tr-Trvid;
%noise power spectral density
[Ptr,f]=periodogram(Tr,[],514,Fd); %512
%periodogram(Tr,[],512,Fd);
figure
psdplot(Ptr,f,’Hz’);
figure
subplot (2,2,1);psdplot(Ptr,f,’Hz’);
%Real signal spectral density
ecg=ecg-Trvid;
ecg=ecg/max(ecg(2100:2800));
[Psig,f]=periodogram(ecg(2100:2800),[],512,Fd); %512
subplot(2,2,2);psdplot(Psig,f,’Hz’);
%Model spectral density
[Pmod,f]=periodogram(yi,[],514,Fd); %512
subplot(2,2,3);psdplot(Pmod,f,’Hz’);
%Wiener coefficients
for ii=1:length(Pmod)
W(ii)=1/(1+(Ptr(ii)/Pmod(ii)));
end
ii=1:length(Pmod);
subplot(2,2,4); plot(ii,W);
grid on;
xlabel(‘f, Hz’);
ylabel(‘W’);
%Low freq filtering
ecgfft = fft(ecg(2100:2800),700);
figure
subplot(2,2,1);plot(abs(ecgfft(1:258)))
grid on;
xlabel(‘f, Hz’);
ylabel(‘Amplitude’);
for ii=1:(length(Pmod))
ecgfd(ii)=W(ii)*ecgfft(ii);
end
subplot(2,2,2);plot(abs(ecgfd));
grid on;
xlabel(‘f, Hz’);
ylabel(‘Amplitude’);
subplot(2,2,3);plot(tecg(1:700),ecg(2101:2800));
grid on;
xlabel(‘t, s’);
ylabel(‘Amplitude, mV’);
ecgiff= (ifft(ecgfd(1:length(Pmod)), 700));
ecgiff=real(ecgiff)/real(max(ecgiff))+Trvid/4;
subplot(2,2,4);plot(tecg(1:700),real(ecgiff));
grid on;
xlabel(‘t, s’);
ylabel(‘Amplitude, mV’);
figure
plot(tecg(1:700),ecg(2101:2800),’red’);
grid on;
xlabel(‘t, s’);
ylabel(‘Amplitude, mV’);
%Filtering in time domain
koef=abs(ifft(W));
ecgtd = filter(koef(1:30),1, ecg);
ecgtd=ecgtd/max(ecgtd);
figure
plot(tecg(1:700),ecgtd(1:700))
grid on;
xlabel(‘t, s’);
ylabel(‘Amplitude, mV’);
figure
subplot(1,2,1);grpdelay(B,A,700,Fd);
subplot(1,2,2);grpdelay(koef(1:30),1,700,Fd);
Files:Â mat1.mat and trm.mat files
