System Identification with a Limited Bandwidth Input Signal and Region of Interest

Question

Given a FIR filter $h[n]$. Its action can described as:

$$ \mathbf{y} = \mathbf{H} \mathbf{x} \\ \mathbf{y} = \mathbf{X} \mathbf{h} $$

where $\mathbf{H}$ and $\mathbf{X}$ is a Toeplitz matrix. If $h$ is unknown, Least Squares with a white Gaussian input signal $x[n]$ can be used to find the unknown coefficients:

$$ \hat{\mathbf{h}} = (\mathbf{X}^{T}\mathbf{X})^{-1} \mathbf{X}^T \mathbf{y} $$

Caveat: $x[n]$ must be white; otherwise the regression matrix $\mathbf{X}^T \mathbf{X}$ is badly conditioned.

The frequency domain information is encoded in the coefficients $h$. However, as can be seen above, the LS algorithm estimates the coefficients with zero prior knowledge; the estimation ONLY depends on the input signal $x$. It does not matter if the system to be identified is an allpass filter or has a notch of 200dB attenuation at $\pi/2$.

Now my question: What do I do if I only care about a small frequency range in $h$ and hence my input signal $x[n]$ does not need to be white?

Example: My Nyquist rate is 10kHz. My unknown system is a lowpass with -3dB at 300 Hz. It has some "weird" frequency behavior around 300 Hz which I want to estimate. I do NOT care about anything beyond, say, 500 Hz. Additionally, my measurement setup prevents me from using a white input signal. I have a bandwidth limitation of 500 Hz. I cannot change the Nyquist rate.

With Least Squares I cannot identify the system because $x$ is not white (persistently exciting). Regularization/SVD does not help me: It provides a biased solution and still gives me $h$ values that try to estimate the entire frequency range. But I really want to say "Give me the $h$ that describes the unknown system best up to 500 Hz with a 500 Hz input signal"

Answer 1

The reason x[n] must be white is because the solution will effectively spectrally weight the channel response based on the amount of energy present in each spectral frequency location. A white noise source provides equal weight to all frequencies. If energy is not present in any particular frequency bin, a proper solution cannot be found for that frequency.

If you only care about a small band of signals, then I would argue that you can still use the least squared approach. The reason for this is to consider a system that provides band pass filtering: for such a system I could excite the input with a white noise source and use the least squares approach (Wiener-Hopf equations) of comparing the input and output signals which would accurately provide a least squares estimate of the channel. So if you have a band limited signal, as long as that signal is white over your band of interest, it will still provide the accurate solution over that band (and you ignore everything else).

Answer 2

Can't you just use the pseudo-inverse? That'll mean instead of:

$$ \hat{\mathbf{h}} = (\mathbf{X}^{T}\mathbf{X})^{-1} \mathbf{X}^T \mathbf{y} $$

you use

$$ \hat{\mathbf{h}}_{\tt pseudo} = (\mathbf{X}^{T}\mathbf{X})^{\dagger} \mathbf{X}^T \mathbf{y} $$

or

$$ \hat{\mathbf{h}}_{\tt pseudo 2} = \mathbf{X}^{\dagger} \mathbf{y} $$

The plot below shows what happens when I do a short example in python (code below). The $x$ in this case is just a sinusoid. Though the length may be too short, it still gives a decent answer whereas inv (the usual inverse) has problems because of the singularity of $\mathbf{X}^{T}\mathbf{X}$.

---

Code below

from numpy import random, zeros, arange, cos
from scipy import pi
from scipy.linalg import toeplitz, inv, pinv
from pylab import figure, clf, plot, xlabel, ylabel, xlim, ylim, title, grid, axes, show, subplot

N = 5
h = [0.2,1,-1,0.6,1]

# x = random.normal(0, 0.01, N)
x = cos(2*pi*0.01234*arange(N) + 2*pi*random.uniform(-1,1))
X = toeplitz(x, zeros(N)) # Need to in fill with zeros.

H = toeplitz(h, zeros(N)) # Need to in fill with zeros.
y = H @ x
y2 = X @ h

h_hat = pinv(X.transpose() @ X) @ X.transpose() @ y
h_hat2 = pinv(X.transpose() @ X) @ X.transpose() @ y2
h_hat3 = pinv(X) @ y

figure(1,  figsize=(20, 6))
subplot(1, 3, 1)
plot(h)
title("True FIR filter")

subplot(1, 3, 2)
plot(y)
plot(y2,'r.')
title("$\mathbf{Xh}$ (red) and $\mathbf{Hx}$ (blue) of filter")

subplot(1, 3, 3)

plot(h)
plot(h_hat,'ro')
plot(h_hat2,'g.')
plot(h_hat3,'k+',markersize=10)
title("True (blue) and estimated (red) filter just pseudo +")