[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