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?

• You talk about a transfer function, so you assume that the system is LTI, right? Why would you have complex output samples (for real-valued input samples and a real-valued system)? Mar 28 '14 at 8:03
• Search for "System Identification". Matlab has a corresponding toolbox. Mar 28 '14 at 8:10
• Is any real system actually LTI? Mar 28 '14 at 18:03
• An LTI system can be a very good approximation to a real system, and they are so much easier to deal with than other more complicated models. Anyway, 'transfer function' implies LTI. Mar 28 '14 at 18:05
• I think you guys might be using the term "real" in two different ways: first, to refer to real-valued signals (as opposed to complex-valued signals), and second, to refer to "realizable", or "real-world". Complex-valued systems are common in some areas of "real-world" signal processing, so I don't see that as a problem with the OP. Mar 28 '14 at 18:54

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')
• Looks promising. Any ideas for how to implement it? Apr 2 '14 at 16:12
• What's happening on the line where you use circshift? How would you put a minimum mse criteria on the weights? Apr 3 '14 at 0:59
• The filter command has an (n-1)/2 sample group delay. The circshift undoes that. Apr 3 '14 at 1:57
• Here's code to efficiently compute the MISO Wiener filter: mathworks.com/matlabcentral/fileexchange/… mathworks.com/matlabcentral/fileexchange/… Apr 3 '14 at 7:03

One approach you can take is to try to fit your data to an ARMA model.

There are several implementations (as that link suggests).

Also a good reference (if you're mathematically minded) is Lennart Ljung's book, System Identification: Theory for the User

Most algorithms mentioned in the links use a mean square error criterion. Many algorithms work well (enough) with small amounts of data.

• Looks promising. Any idea how to implement it? Apr 3 '14 at 1:00