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 w_k

Further wee need to optimize filter coefficients w_k 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;

S_d(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:

butterworth_tap4_cut30Hz.jpg

butterworth_tap4_cut70Hz.jpg

butterworth_tap8_cut30Hz.jpg

butterworth_tap8_cut70Hz.jpg

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

Blogsphere: TechnoratiFeedsterBloglines
Bookmark: Del.icio.usSpurlFurlSimpyBlinkDigg
RSS feed for comments on this post  |  TrackBack URI for this post

This entry was posted on Tuesday, January 2nd, 2007 at 2:49 pm and is filed under DSP Filtering. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

New on WinAVR Tutorial Running TX433 and RX433 RF modules with AVR microcontrollers, Sometimes in embedded design you may want to go wireless. Might be you will want to log various readi … Weiterlesen » Programming AVR ADC module with WinAVR, Most of AVR microcontrollers have Analog to Digital Converter (ADC) integrated in to chip. Such solut … Weiterlesen »

New on WinARM Tutorial What are differences between WinARM and WinAVR, Everyone who is working with AVR microcontrollers knows this powerful tool – WinAVR (http://win … Weiterlesen » LPC2000 watchdog timer, As in all microcontrollers watchdog timers purpose isto reset microcontroller after reasonable amount … Weiterlesen »