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 & x_2 & ... & x_r \\
x_2 & x_3 & & x_{r+1} \\
x_3 & x_4 & & x_{r+2} \\
... & & & ...\\
x_n & x_{n+1} & ... & 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')