Wavelet domain wiener filter matlab implementation regarding

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

On Apr 5, 12:25&#4294967295;pm, "rsram" <raghu.coolsriram@n_o_s_p_a_m.gmail.com>
wrote:
> Dear Sir/Madam,
> &#4294967295; &#4294967295; &#4294967295; &#4294967295; 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;

On Apr 5, 1:25&#4294967295;pm, "rsram" <raghu.coolsriram@n_o_s_p_a_m.gmail.com>
wrote:
> Dear Sir/Madam,
> &#4294967295; &#4294967295; &#4294967295; &#4294967295; 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

On Apr 5, 12:25&#4294967295;pm, "rsram" <raghu.coolsriram@n_o_s_p_a_m.gmail.com>
wrote:
> Dear Sir/Madam,
> &#4294967295; &#4294967295; &#4294967295; &#4294967295; 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;