Fitting high order rational function to frequency response measurement data

Question

I need to model measurement data of a frequency response with physically meaningful band-pass filters. All Measurements happening in die range of 20Hz to 20000Hz

I have to work with python.

The Problem:

Imagine a rational transfer function as an undefined number of second order band-pass filters in parallel:

$$H(s) = \sum_{n=1}^{N}(\frac{b_{n1} s}{a_{n0} + a_{n1} s + a_{n2} s^{2}})\tag{1}$$

I'm looking for an algorithm that takes a frequency response measurement $f(s)$ to calculate all coefficients in (1) while finding the correct order of $f(s)$.

What I did so far:

Least-squares

$$H(s) = f(s)$$ or $$ \frac{\mbox{num}(H(s))}{\mbox{den}(H(s))} = f(s) $$ Solving $Ax = b$ with least squares.

As an example, for just one band-pass filters: $$A_k = [s_{k}b_1, -f(s_k)a_0 , -f(s_k) s_k a_1, -f(s_k) s_{k}^{2} a_2]$$

and

$$ b_k = f(s_k)$$

and

$$ x = [b1,a_0,a_1,a_2]$$

This is very useful for $ N<2 $ or order 4 of $ H(s) $ respectively. For higher orders, the fitting gets bad because of the high orders of $s$ being multiplied with $f(s)$

Same problem does exsit with the invfreqs() function know from matlab.

Vector Fitting

As shown in this paper by Bjørn Gustavsen, a rational function can be fitted to high accuracy with a model like:

$$H(s) = \sum_{n=1}^{N}(\frac{c_{n}}{(s-p_n)})\tag{2}$$

For resonant peaks in the measurement data, the poles $p_n$ will form complex conjugate pairs. The same is true for the corresponding residues $c_n$

As an exsample for oder 2, this leads to a second order system of the following form:

$$H(s) = \frac{b_0 + b_{1} s}{a_{0} + a_{1} s + a_{2} s^{2}}\tag{3}$$

As you can see, the extra coefficient $b_0$ allows the system to perform a low-pass behavior. In Fact, you can describe this filter as a low-pass and a band-pass in parallel, both share the same poles.

I cant ensure that the Vector Fitting converges on a set of band-pass filters. Especially on band-pass filters with positive gain at the resonance frequency (which is a must for physically meaningful filters)

As you can see in this plot of $|H(s)|$:

The blue curve has been generated by 4 band-pass filters in parallel and is used as input data. I split the results in it's band-pass (BP) and low-pass (LP) parts to show how the response is combined by all of these components. Note the negative gain of the purple band-pass filter and the cancellation of all low-pass filters in low frequencies.

$H(s)$ reassembles $f(s)$ perfectly, but not in a physically meaningful way

My Question:

Can anyone help me find a good working algorithm for my problem? Or do you spot possibilities to modify the Vector Fitting algorithm to work for my problem?

Thanks a LOT!

Answer 1

I think it should be possible to modify the "Vector Fitting" algorithm in a way such that the individual first order filters add up to proper (second-order) bandpass filters, without the lowpass component.

Note that in any case, for each first order section (leaving out the section index for convenience)

$$\frac{c}{s-p}\tag{1}$$

you must have another section of the form

$$\frac{c^*}{s-p^*}\tag{2}$$

We assume that the pole $p$ is complex-valued, with non-zero imaginary part. Combining $(1)$ and $(2)$ results in

$$\frac{c}{s-p}+\frac{c^*}{s-p^*}=\frac{s(c+c^*)-(cp^*+c^*p)}{(s-p)(s-p^*)}\tag{3}$$

In order for $(3)$ to be a second-order bandpass transfer function, we have to make sure that

$$cp^*+c^*p=2\textrm{Re}\{cp^*\}=0\tag{4}$$

which is the case if

$$c_Rp_R+c_Ip_I=0\tag{5}$$

holds, where the subscripts $_R$ and $_I$ denote real and imaginary parts, respectively.

Using condition $(5)$, we can rewrite $(1)$ as

$$\frac{c_R-ic_R\frac{p_R}{p_I}}{s-p}\tag{6}$$

Of course, $(2)$ must be modified accordingly.

Note that the subproblems solved in each iteration of the algorithm are still linear in the parameters, because the pole locations are fixed per iteration.