4
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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)? $\endgroup$ – Matt L. Mar 28 '14 at 8:03
  • 1
    $\begingroup$ Search for "System Identification". Matlab has a corresponding toolbox. $\endgroup$ – Matt L. Mar 28 '14 at 8:10
  • $\begingroup$ Is any real system actually LTI? $\endgroup$ – Seth Mar 28 '14 at 18:03
  • $\begingroup$ 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. $\endgroup$ – Matt L. Mar 28 '14 at 18:05
  • $\begingroup$ 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. $\endgroup$ – Jason R Mar 28 '14 at 18:54
4
$\begingroup$

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

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.

$\endgroup$
  • $\begingroup$ Looks promising. Any idea how to implement it? $\endgroup$ – Seth Apr 3 '14 at 1:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.