Transfer function determination from input and output data

Question

I have some input and output data that I believe adequately includes excitation of the important dynamics of a system. I know it is at most a 4th-order transfer function.

How can I identify the transfer function? I have Python available. (I have MATLAB but do not have access to the System ID Toolbox and need to limit my dependencies on MATLAB)

I know how to do least-squares fits, but I have no idea how to apply the idea to transfer functions, where using the absolute error in the transfer function doesn't make sense (relative error seems to be the appropriate thing to fit).

I'm vaguely familiar with ARX/ARMA/ARMAX approaches but they always seem to overemphasize the high-frequency response, because they solve for the impulse response.

Answer 1

if your input excitation signal (call it $x[n]$) is arbitrarily long and so also your output $y[n]$, then you should be able to take as long of block (let's call that length $N$) of each and, using the DFT and division compute $$ H[k] = \frac{Y[k]}{X[k]} $$

in audio, we call that method of system identification the "dual-channel FFT". the DFT of $h[n]$ is equal to the the DTFT of $h[n]$ (wrapped-around and overlapped and added every $N$ samples) which is the Z-transform of the same wrapped (and aliased) $h[n]$ sampled at $N$ equally-spaced points around the unit circle. $$ H(e^{j 2 \pi k/N}) = H[k] $$

so if your 4th-order IIR transfer function is $$ G(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + b_4 z^{-4}}{1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + a_4 z^{-4}} $$

what you want to do is choose $ \{ b_0, b_1, b_2, b_3, b_4, a_1, a_2, a_3, a_4 \} $ in such a way to minimize $$ \sum_{k=0}^{N-1} \left| G(e^{j 2 \pi k/N}) - H(e^{j 2 \pi k/N}) \right|^2 \ = \ \sum_{k=0}^{N-1} \left| G(e^{j 2 \pi k/N}) - H[k] \right|^2 $$ or, perhaps a weighted version of that, $$ \sum_{k=0}^{N-1} W[k] \left| G(e^{j 2 \pi k/N}) - H[k] \right|^2 $$ where $W[k]$ is known.