estimation - Estimate the Transfer Function of an Unknown System

Suppose you have a system, H, that you want to estimate its transfer function.

You have a finite number of complex input samples, x, and noisy complex (magnitude and phase) output samples, y:

1. 
   

   In practice, how do you estimate it's transfer function given a mean-squared error criteria?

   
   



2. 
   

   What algorithms are used to estimate transfer functions when you have limited and noisy data?

   
   


3. 
   

   How would you estimate the transfer function using a Wiener filter?

   
   

Answer

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')