Identify continuous-time filter parameters from frequency response data


[b,a] = invfreqs(h,w,n,m)
[b,a] = invfreqs(h,w,n,m,wt)
[b,a] = invfreqs(h,w,n,m,wt,iter)
[b,a] = invfreqs(h,w,n,m,wt,iter,tol)
[b,a] = invfreqs(h,w,n,m,wt,iter,tol,'trace')
[b,a] = invfreqs(h,w,'complex',n,m,...)


invfreqs is the inverse operation of freqs. It finds a continuous-time transfer function that corresponds to a given complex frequency response. From a laboratory analysis standpoint, invfreqs is useful in converting magnitude and phase data into transfer functions.

[b,a] = invfreqs(h,w,n,m) returns the real numerator and denominator coefficient vectors b and a of the transfer function


whose complex frequency response is given in vector h at the frequency points specified in vector w. Scalars n and m specify the desired orders of the numerator and denominator polynomials.

The length of h must be the same as the length of w. invfreqs uses conj(h) at -w to ensure the proper frequency domain symmetry for a real filter.

[b,a] = invfreqs(h,w,n,m,wt) weights the fit-errors versus frequency, where wt is a vector of weighting factors the same length as w.

[b,a] = invfreqs(h,w,n,m,wt,iter) and

[b,a] = invfreqs(h,w,n,m,wt,iter,tol) provide a superior algorithm that guarantees stability of the resulting linear system and searches for the best fit using a numerical, iterative scheme. The iter parameter tells invfreqs to end the iteration when the solution has converged, or after iter iterations, whichever comes first. invfreqs defines convergence as occurring when the norm of the (modified) gradient vector is less than tol, where tol is an optional parameter that defaults to 0.01. To obtain a weight vector of all ones, use


[b,a] = invfreqs(h,w,n,m,wt,iter,tol,'trace') displays a textual progress report of the iteration.

[b,a] = invfreqs(h,w,'complex',n,m,...) creates a complex filter. In this case no symmetry is enforced, and the frequency is specified in radians between –π and π.


collapse all

Convert a simple transfer function to frequency-response data and then back to the original filter coefficients.

a = [1 2 3 2 1 4];
b = [1 2 3 2 3];

[h,w] = freqs(b,a,64);
[bb,aa] = invfreqs(h,w,4,5)
bb = 1×5

    1.0000    2.0000    3.0000    2.0000    3.0000

aa = 1×6

    1.0000    2.0000    3.0000    2.0000    1.0000    4.0000

bb and aa are equivalent to b and a, respectively. However, the system is unstable because aa has poles with positive real part. View the poles of bb and aa.


Use the iterative algorithm of invfreqs to find a stable approximation to the system.

[bbb,aaa] = invfreqs(h,w,4,5,[],30)
bbb = 1×5

    0.6816    2.1015    2.6694    0.9113   -0.1218

aaa = 1×6

    1.0000    3.4676    7.4060    6.2102    2.5413    0.0001

Verify that the system is stable by plotting the new poles.


Generate two vectors, mag and phase, that simulate magnitude and phase data gathered in a laboratory. Also generate a vector, w, of frequencies.


fs = 1000;
t = 0:1/fs:2;
mag = periodogram(sin(2*pi*100*t)+randn(size(t))/10,[],[],fs);
phase = randn(size(mag))/10;
w = linspace(0,fs/2,length(mag))';

Use invfreqs to convert the data into a continuous-time transfer function. Plot the result.

[b,a] = invfreqs(mag.*exp(1j*phase),w,2,2,[],4);



When building higher order models using high frequencies, it is important to scale the frequencies, dividing by a factor such as half the highest frequency present in w, so as to obtain well conditioned values of a and b. This corresponds to a rescaling of time.


By default, invfreqs uses an equation error method to identify the best model from the data. This finds b and a in


by creating a system of linear equations and solving them with the MATLAB® \ operator. Here A(w(k)) and B(w(k)) are the Fourier transforms of the polynomials a and b, respectively, at the frequency w(k), and n is the number of frequency points (the length of h and w). This algorithm is based on Levi [1]. Several variants have been suggested in the literature, where the weighting function wt gives less attention to high frequencies.

The superior (“output-error”) algorithm uses the damped Gauss-Newton method for iterative search [2], with the output of the first algorithm as the initial estimate. This solves the direct problem of minimizing the weighted sum of the squared error between the actual and the desired frequency response points.



[1] Levi, E. C. “Complex-Curve Fitting.” IRE Trans. on Automatic Control. Vol. AC-4, 1959, pp. 37–44.

[2] Dennis, J. E., Jr., and R. B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations.Englewood Cliffs, NJ: Prentice-Hall, 1983.

See Also

| | |

Introduced before R2006a