Estimate the Transfer Function of an Unknown System
Here's the way I think about a discrete Wiener Filter

Consider a sequence of observations $\mathbf{y} \in \Re^n $

Form a matrix from the input $\mathbf{x} \in \Re^{n+r-1}$ by shifting columns one sample each:
$$
X= \begin{bmatrix}
x_1 &amp; x_2 &amp; ... &amp; x_r \\ 
x_2 &amp; x_3 &amp; &amp; x_{r+1} \\ 
x_3 &amp; x_4 &amp; &amp; x_{r+2} \\
... &amp; &amp; &amp; ...\\
x_n &amp; x_{n+1} &amp; ... &amp; x_{r+n-1}\\
\end{bmatrix}
$$

Then the act of correlating the signal $\mathbf{x}$ with a filter $\mathbf{h}$ can be represented as $\mathbf{y} = \mathbf{X h}$

So $\mathbf{X^Ty}=\mathbf{X^TX~h}$ and the least squares solution for $\mathbf{h}=(\mathbf{X^TX})^{-1}\mathbf{X^Ty}$

Here's an octave/Matlab example that might help illustrate

r = 31; % filter length
n = 1e4;% signal length
SNR = 0;% AWGN measurement SNR in dB
x = randn(n,1); % random signal
h0 = randn(r,1); % unknown filter
y0 = filter(h0,1,x); % convolve x0,h0
y0 = circshift( y0, -(r-1)/2 ); % align x,y0
y = y0 + randn(size(y0)) * norm(y0)/sqrt(n) * 10^(-SNR/20); % add noise

% Wiener-Hopf solution
XX = xcorr(x,x,r-1);
Xy = xcorr(x,y,(r-1)/2 );
hrev = toeplitz(XX(r:end))  \ Xy; % least squares solution
h1 = conj(flipud(hrev));  % change from a correlating filter to a convolution filter

nmse_recovery = 20*log10( norm(h1-h0) / norm(h0));
plot([h0 h1]);
legend('truth','recovered')