Reply by illywhacker April 6, 20112011-04-06
On Apr 5, 12:25&#2013266080;pm, "rsram" <raghu.coolsriram@n_o_s_p_a_m.gmail.com>
wrote:
> Dear Sir/Madam, > &#2013266080; &#2013266080; &#2013266080; &#2013266080; I need to perform wavelet domain wiener filtering to denoise ecg > signal as a part of my project > > First I performed wiener filtering in fourier domain > The code is as follows: > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > function output = Wienerfilt(noisy,clean) > > nois=clean-noisy; > stdnois=std(nois); > > y=fft(clean); > fftnsy=fft(noisy); > n=length(clean); > > syy=abs(y.*y); > H=syy./(syy+ n*stdnois^2);%Wiener filter transfer function > fftout=fftnsy.*H; > output=real(ifft(fftout)); > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > I wanted to extend this to wavelet domain. > My code for this is > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > function output = WienerWavelet(noisy,clean) > > [YA,YD] = dwt(clean,'coif5');%approx & detail coefficients of clean signal > [NA,ND] = dwt(noisy,'coif5');%approx and detail coeffs of noisy signal > > nois=clean-noisy; > stdnois=std(nois);%standard deviation of noise > n=length(clean); > > PAY=abs(YA.*YA);PAD=abs(YD.*YD); > > HA=PAY./(PAY+(n/2)*stdnois^2); > HD=PAD./(PAD+ (n/2)*stdnois^2); > > NA=NA.*HA; > ND=ND.*HD; > > output = idwt(NA,ND,'coif5'); > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > The problem that I am experiencing is that except R- wave(the spike like > wave of QRS complex) all other characteristics are lost in the process of > denoising. > > Is this approach of wiener filtering in the wavelet domain wrong? > Is there any algorithm or code for wiener filtering of 1-D signals in > wavelet domain? > > Please tell me the corrections to be made in this code
Rune raises a good point: do you really just want to implement the Wiener filter in the wavelet domain? This is silly, not hard in principle, but computationally messy. Or do you want to implement something analogous to the Wiener filter, using a signal model analogous to that used in the standard Wiener filter? What you are doing now in the wavelet domain is not the classical Wiener filter at all, because the signal model behind it is different. In the Fourier case, the Gaussian signal model is translation invariant, hence the covariance is diagonalized in the Fourier basis and the coefficients are all independent. In the wavelet case, modulo the points I mentioned about the mean, your prior model is again Gaussian, but now the wavelet coefficients are assumed to be independent. This is unusual, because it does not even try to impose translation invariance. It is certainly different from the classical Wiener filter. illywhacker;
Reply by Rune Allnor April 6, 20112011-04-06
On Apr 5, 1:25&#2013266080;pm, "rsram" <raghu.coolsriram@n_o_s_p_a_m.gmail.com>
wrote:
> Dear Sir/Madam, > &#2013266080; &#2013266080; &#2013266080; &#2013266080; I need to perform wavelet domain wiener filtering to denoise ecg > signal as a part of my project
What the heck is a 'wavelet domain Wiener filter'? Wavelets are nothing but a basis on which to project your signal. Nothing more nothing less. The Fourier spectrum is also a basis for the signal, but one that coincidentially has the added property that analysis of time domain operations like convolutions can easily be analyzed, interpreted, and even specified. You need to work out all the issues you otherwise find in the ubiquotous chapter 'Properties of the Fourier Transform' in DSP and system analysis texts, for the wavelet case before you can begin the task you 'need' to do. You then have to come up with a convincing argument that the excercise is worth the while. Last, you need to express the Wiener filter in this new framework of analysis. You might want to consider if there might be a reason why no one have cared about doing any of this in the past. Rune
Reply by illywhacker April 5, 20112011-04-05
On Apr 5, 12:25&#2013266080;pm, "rsram" <raghu.coolsriram@n_o_s_p_a_m.gmail.com>
wrote:
> Dear Sir/Madam, > &#2013266080; &#2013266080; &#2013266080; &#2013266080; I need to perform wavelet domain wiener filtering to denoise ecg > signal as a part of my project > > First I performed wiener filtering in fourier domain > The code is as follows: > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > function output = Wienerfilt(noisy,clean) > > nois=clean-noisy; > stdnois=std(nois); > > y=fft(clean); > fftnsy=fft(noisy); > n=length(clean); > > syy=abs(y.*y); > H=syy./(syy+ n*stdnois^2);%Wiener filter transfer function > fftout=fftnsy.*H; > output=real(ifft(fftout)); > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > I wanted to extend this to wavelet domain. > My code for this is > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > function output = WienerWavelet(noisy,clean) > > [YA,YD] = dwt(clean,'coif5');%approx & detail coefficients of clean signal > [NA,ND] = dwt(noisy,'coif5');%approx and detail coeffs of noisy signal > > nois=clean-noisy; > stdnois=std(nois);%standard deviation of noise > n=length(clean); > > PAY=abs(YA.*YA);PAD=abs(YD.*YD); > > HA=PAY./(PAY+(n/2)*stdnois^2); > HD=PAD./(PAD+ (n/2)*stdnois^2); > > NA=NA.*HA; > ND=ND.*HD; > > output = idwt(NA,ND,'coif5'); > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > > The problem that I am experiencing is that except R- wave(the spike like > wave of QRS complex) all other characteristics are lost in the process of > denoising. > > Is this approach of wiener filtering in the wavelet domain wrong? > Is there any algorithm or code for wiener filtering of 1-D signals in > wavelet domain? > > Please tell me the corrections to be made in this code > > Please help me in this regard > My email ID is raghu.coolsri...@gmail.com > > Thanks in advance > > Yours sincerely > > Sri Ram Raghu
- What are you doing about the zero-frequency component in the Fourier case? - In the wavelet case, the individual approximation coefficients may have finite means, although translation invariance means that these must all be the same. How are you estimating and allowing for that in your code? (Note that I do not mean 'average', but 'mean'. You are implicitly assuming a prior probability model for your clean signal: the mean I am talking about is a parameter of that distribution.) illywhacker;
Reply by rsram April 5, 20112011-04-05
Dear Sir/Madam,
        I need to perform wavelet domain wiener filtering to denoise ecg
signal as a part of my project

First I performed wiener filtering in fourier domain
The code is as follows:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function output = Wienerfilt(noisy,clean)

nois=clean-noisy;
stdnois=std(nois);

y=fft(clean);
fftnsy=fft(noisy);
n=length(clean);

syy=abs(y.*y);
H=syy./(syy+ n*stdnois^2);%Wiener filter transfer function
fftout=fftnsy.*H;
output=real(ifft(fftout));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I wanted to extend this to wavelet domain.
My code for this is

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function output = WienerWavelet(noisy,clean)

[YA,YD] = dwt(clean,'coif5');%approx & detail coefficients of clean signal
[NA,ND] = dwt(noisy,'coif5');%approx and detail coeffs of noisy signal

nois=clean-noisy;
stdnois=std(nois);%standard deviation of noise
n=length(clean);

PAY=abs(YA.*YA);PAD=abs(YD.*YD);

HA=PAY./(PAY+(n/2)*stdnois^2);
HD=PAD./(PAD+ (n/2)*stdnois^2);

NA=NA.*HA;
ND=ND.*HD;

output = idwt(NA,ND,'coif5');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The problem that I am experiencing is that except R- wave(the spike like
wave of QRS complex) all other characteristics are lost in the process of
denoising.

Is this approach of wiener filtering in the wavelet domain wrong?
Is there any algorithm or code for wiener filtering of 1-D signals in
wavelet domain?

Please tell me the corrections to be made in this code

Please help me in this regard 
My email ID is raghu.coolsriram@gmail.com

Thanks in advance


Yours sincerely

Sri Ram Raghu