Find coefficients of optimal Wiener filter of length 2

I need to find the coefficients (impulse response) of a FIR Wiener filter with length equal to 2. I have a gaussian white noise signal that is generated using the Standard Normal Distribution (mean = 0, var = 1). I also have a seperate system that "colors" the white noise signal. It alters the PSD of the signal so that is it not uniform.

Then I need to input the "colored" noise signal into the FIR Wiener filter, which needs to approximate the original white noise signal,w(n).

Let's call this white noise approximation w_hat(n). So I want to: w_hat(n) ~ w(n), where w_hat(n) is the output of the Wiener filter.

The colored noise signal is the following : v(n) = 0.6v(n-1) + w(n) , where w(n) is the gaussian white noise signal.

I know that, theoretically, the optimal Wiener filter, can be found by solving the Wiener-Hopf equation, i.e minimizing the expected value of the MSE between w(n) and w_hat(n), over the filter's impulse response. However I am having a hard time implementing this in MATLAB.

We also know that this is the Frequency Response of the optimal Wiener filter (sorry can't post an image due to low reputation),

where in the numerator there is the PSD of w(n) and v(n), (i.e the DTFT of their crosscorrelation) and in the denominator is the PSD of v(n).

What I have tried so far is the following :

1. Compute the croscorrelation between the (ideal) white noise w(n) and the filter's input v(n), using MATLAB's xcorr, store the result into r_wv
2. Finding its DTFT using abs(fftshift(fft(r_wv))), store into phi_wv
3. Compute the autocorrelation of the input signal v(n), using xcorr(), store the result into r_v
4. Finding its DTFT like in step 2, store into phi_vv
5. Get the filter's Frequency Response, as mentioned above, by doing phi_vv ./ phi_wv, store the result into freq_resp.
6. Finally, apply IDTFT, in order to compute the filter's impulse response. I use abs(ifftshift(ifft(freq_resp)))

However the impulse response I get this way is a complex one, and not a real one, and more importantly it has the same length as the signals used (1001 samples), and not the desired length of 2.

EDIT: I know that I have to solve this linear system. However, I am confused as to how I would compute, say Rvv(0) .

Any help, is greatly appreciated as I have been stuck on this for a few days.

The code :

n = 0 : 1000; % number of samples
w = randn(1, length(n)); %white noise signal
v = filter(1, [1, -0.6], w); % the "colored" white noise signal

[r_v, lags] = xcorr(v); %compute autocorrelation of v(n)
phi_vv = abs(fftshift(fft(r_v))); %compute its DTFT, i.e PSD
[r_wv,lags2] = xcorr(w,v); %compute crosscorrelation of w(n) and v(n)
phi_wv = abs(fftshift(fft(r_wv))); % compute its DTFT
freq_resp = phi_wv ./ phi_vv;       % compute the filter's Frequency Response
hW = abs(ifftshift(ifft(freq_resp))); % use IDTFT in order to get its impulse response
w_hat = filter(hW, 1, v); % filter the "colored" noise signal, to get an approximation of the white noise


You have a input signal $$x(n)$$ (the colored noise in your case) and a desired output signal $$d(n)$$. Assuming that the FIR coefficients are real. Define the filter weight vector $$\mathbf{w}$$ and the input vector $$\mathbf{x}(n)$$ $$\mathbf{w} = [w_0, w_1, \ldots, w_{N-1}]^T$$ $$\mathbf{x}(n) = [x(n), x(n-1), \ldots, x(n-N+1)]^T$$ The output of the filter is $$y(n) = \sum_{i=0}^{N-1}w_i x(n-i) = \mathbf{w}^T\mathbf{x}(n) = \mathbf{x}^T(n)\mathbf{w}$$ The error is $$e(n) = d(n) - y(n)$$ and the MSE is \begin{aligned} \xi &= E[e^2(n)] = E\{[d(n) - \mathbf{w}^T\mathbf{x}(n)] [d(n) - \mathbf{x}^T(n)\mathbf{w}]\}\\ &=E[d^2(n)] - \mathbf{w}^T E[\mathbf{x}(n)d(n)] - E[d(n)\mathbf{x}^T(n)]\mathbf{w} + \mathbf{w}^T E[\mathbf{x}(n)\mathbf{x}^T(n)]\mathrm{w} \end{aligned} We then define the cross-correlation vector $$\mathbf{p} = E[\mathbf{x}(n)d(n)]$$ and the autocorrelation matrix $$\mathbf{R} = E[\mathbf{x}(n)\mathbf{x}^T(n)]$$ we have $$\xi = E[d^2(n)] - 2\mathbf{w}^T\mathbf{p} + \mathbf{w}^T \mathbf{Rw}$$ To minimize the MSE, let $$\frac{\partial \xi}{\partial \mathbf{w}} = 2\mathbf{Rw} - 2\mathbf{p} = \mathbf{0}$$ The Wiener filter is finally derived as $$\mathbf{w} = \mathbf{R}^{-1}\mathbf{p}$$

Here's a MATLAB implementation

L = 1000; % number of samples
d = randn(1, L); % desired signal
x = filter(1, [1, -0.6], d); % input signal

N = 2;  % filter length
p = zeros(N, 1); % cross-correlation vector
R = zeros(N, N); % autocorrelation matrix

for n = N:L
xx = x(n:-1:n-N+1).'; % input vector
p = p + xx * d(n);
R = R + xx * xx.';
end
p = p / L;
R = R / L;
w = R \ p;


which gives the result w = [1.0000; -0.6000]

• Thank you so much !! I was having trouble implementing the math, into MATLAB code. Jul 6 at 13:11
• One question though.. If I want to do the same, but for say a filter of length 3 or 4, if I just change N to say 3, there seem to be problems with the matrix dimensions I think that the following line must be changed into : xx = [x(n), x(n-1), x(n-2)].'  Jul 6 at 13:14
• Then change first line in the for loop to xx = x(n:-1:n-N+1).’; Keep the definition in mind, this vector has a length of N. Jul 6 at 13:18