function [w] = fWiener(x1,x2,p,N)
% FWIENER Calculates and returns the FIR filter's coefficients using Wiener
% method.
%
% [W] = FWIENER(X1,X2,P,N) where (W) stores Wiener's FIR filter's
% coefficients, (X1) corresponds to, in the noise cancellation case, the
% contaminated signal and (X2) to the noise proceeding from the noise
% source. We need (X1) and (X2) to be the more correlated possible. The
% filter order is (P) and (N) is the number of samples used to obtain the
% filter's coefficients.
% Copyright (C) Iván López Espejo 2009
% Version: Id: fWiener.m , v 1.0 2009
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You can obtain a copy of the GNU General Public License from
% ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to
% Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We consider real samples, so we do not take into account the conjugate of
% the autocorrelation matrix.
% Calculation of the autocorrelation of (X1) and the correlation between
% (X2) and (X1).
x1 = x1(1:N);
x2 = x2(1:N);
autocorrx1 = xcorr(x1,'biased');
autocorrx1x2 = xcorr(x2,x1,'biased');
% Calculation of the Wiener's matrix and independent coefficients vector.
vWiener = autocorrx1(ceil(length(autocorrx1)/2):ceil(length(autocorrx1)/2)+p-1);
mWiener = toeplitz(vWiener);
vWiener = autocorrx1x2(ceil(length(autocorrx1x2)/2):ceil(length(autocorrx1x2)/2)+p-1);
vWiener = vWiener';
% Wiener's FIR filter's coefficients.
w = inv(mWiener)*vWiener;