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

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"

• The only way you can restrict certain frequencies, is by designing $H$ as matrix that performs filtering operation. In other words, if $H$ is a FIR low-pass Butterworth filter, then $H$ would be a Toeplitz matrix and the entries of $H$ would be the filter coefficients of the low-pass filter. Jul 21 '19 at 1:21
• I think you misunderstand me. $\mathbf{X}$ is a Toeplitz matrix. Least Squares can not solve this problem because as soon as use a non-white signal $x[n]$, $\mathbf{X}$ becomes very badly conditioned. My question is, which method (non Least Squares) can be used to restrict the solution to a particular frequency range?
– divB
Jul 22 '19 at 5:02
• PS: I also corrected a minor error in my posting
– divB
Jul 22 '19 at 5:05
• I don't fully understand what you mean by "frequency range in $h$" (seems $h$ are time-domain coefficients?) but in general when $X$ is not ideally conditioned, you may resort to some kind of $\min_h \|y-Xh\| + \lambda \cdot r(h)$, where $r(h)$ is a regularizer that regularizes towards your contraint set... Jul 22 '19 at 7:24
• $h$ implicitely encode the frequency behavior. But anyway, I edited my posting again and hope to have clarified what I mean.
– divB
Jul 22 '19 at 8:25

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

• I think some regularization is needed. I haven't had enough time to think about it. So I don't have the proper regularization. It might be better to analyze it in the Fourier Domain (So regularization will be done with multiplication with some rows of the DFT Matrix).
– Royi
Mar 7 '20 at 18:06
• @Dan Boschen: I think you just repeated my question in different words (at least, in my opinion) ;-) It is precicely my question: HOW to do the regularization (if any)? I mentioned that various regularization attempts I tried did not help me (so far).
– divB
Mar 9 '20 at 15:50
• @divB my point is if your x[n] is sufficiently white in your band of interest then no further regularization is necessary and the results within that band should be accurate to the least squared solution just as much as if x[n] was white using the logic I provided. Mar 9 '20 at 15:55
• Sorry if that wasn’t clear Mar 9 '20 at 15:56
• @DanBoschen: Hmm, I think this is the part where I'd appreciate if you could elaborate a bit, maybe write explicitely what you mean or give an example. In my post I show that LS estimates the vector using the matrix X'X. This matrix is NOT well conditioned if x[n] is not sufficiently white; in fact, LS does not know anything about the "bandwidth of interest" or the system to be identified! Whether it's an allpass or has -100dB at fs/4 ... (X'X) is independent of that. I just don't see how (X'X) becomes better conditioned with non-white-ish x[n].
– divB
Mar 10 '20 at 4:22

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 +")

• No need to use the pseudo inverse on ${\left( {X}^{T} X \right)}$. You can use ${X}{^\dagger}$ directly.
– Royi
Apr 23 '20 at 19:40
• @Royi D'oh! Added that too.
– Peter K.
Apr 23 '20 at 21:51