help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Wiener-Kolmogorov filter


From: Fredrik Lingvall
Subject: Re: Wiener-Kolmogorov filter
Date: Thu, 12 Nov 2009 11:03:56 +0100
User-agent: Thunderbird 2.0.0.23 (X11/20091019)

address@hidden wrote:
> Hi friends,
> I would like to ask you if there is an Octave program for wiener-Kolmogorov 
> filter for signal extraction.
> Thank you
>
>   

% Model:
%
%  y = conv(h,u) + e
%
% where it is assued that you know the impulse response h, the noise
variance (sigma_e^2)
% and  input signal variance (sigma_u^2).

sigma_e = 1.0;
sigma_u = 1.0;
mu =  (sigma_e/sigma_u)^2;

% Length of the estimated input signal.
u_len = 1000;

1) Frequency-domain deconvolution (Wiener filtering):

h_c = fft(h,length(h)+length(y)-1); % Zero-pad to avoid aliasing.
yc = fft(y,length(h)+length(y)-1); % Zero-pad to avoid aliasing.
u_hat = real(ifft( (conj(h_c) .* yc) ./ (h_c.*conj(h_c) + mu)) );
u_hat = u_hat(1:u_len);

2) Time-domain version (u_hat = inv(H'*H + mu*I)*H'*y):

% The matched filter part (U'*y):

% Alt. 1:
%h_t = [h_hat(:)' zeros(1,u_len-1)];
%H = toeplitz(h_t,[h_t(1) zeros(1,u_len-1)]);
%Hty = H'*y;
 
% Alt. 2:
%u_t = [h(:)' zeros(1,u_len-1)];
%U = toeplitz(h_t,[h_t(1) zeros(1,u_len-1)]);
%Hty = (y'*H)'; % More mem efficient.
 
% Alt. 3:
hty = fftconv(h_hat(end:-1:1),y);  % Even faster fft-based method.
Hty = hty(length(h_hat):length(h_hat)+u_len-1); % Pick the relevant part.

% The  H'*H + mu*I part:

I = eye(u_len,u_len);
HtH = H'*H + mu*I;

% The Wiener (Linear MMSE) solution:

u_hat = cholinv(HtH) * Hty;

/Fredrik





reply via email to

[Prev in Thread] Current Thread [Next in Thread]