Identify continuous-time filter parameters from frequency response data
Transfer Function to Frequency Response Conversion
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
aa are equivalent to
a, respectively. However, the system is unstable because
aa has poles with positive real part. View the poles of
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.
Continuous-Time Transfer Function
Generate two vectors,
phase, that simulate magnitude and phase data gathered in a laboratory. Also generate a vector,
w, of frequencies.
rng('default') 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))';
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); freqs(b,a)
h — Frequency response
Frequency response, specified as a vector.
m — Desired order
positive integer scalar
Desired order of the numerator and denominator polynomials, specified as positive integer scalars.
wt — Weighting factors
Weighting factors, specified as a vector.
wt is a vector of
weighting factors that is the same length as
iter — Number of iterations in the search algorithm
positive real scalar
Number of iterations in the search algorithm, specified as a positive real scalar.
iter parameter tells
invfreqs to end the
iteration when the algorithm has converged to a solution, or after
iter iterations, whichever occurs first.
tol — Tolerance
0.01 (default) | scalar
Tolerance, specified as a scalar.
invfreqs defines convergence
as occurring when the norm of the (modified) gradient vector is less than
To obtain a weight vector of all ones, use
a — Transfer function coefficients
Transfer function coefficients, returned as vectors. Express the transfer function
in terms of
b = [1 3 3 1]/6 and
a = [3 0 1 0]/3
specify a third-order Butterworth filter with normalized 3 dB frequency
Complex Number Support: Yes
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
b. This corresponds to a rescaling of time.
invfreqs uses an equation error method to identify the best
model from the data. This finds
by creating a system of linear equations and solving them with the MATLAB®
\ operator. Here
B(w(k)) are the Fourier transforms
of the polynomials
b, respectively, at the
frequency w(k), and n is the number
of frequency points (the length of
algorithm is based on Levi . Several variants have
been suggested in the literature, where the weighting function
less attention to high frequencies.
The superior (“output-error”) algorithm uses the damped Gauss-Newton method for iterative search , 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.
 Levi, E. C. “Complex-Curve Fitting.” IRE Trans. on Automatic Control. Vol. AC-4, 1959, pp. 37–44.
 Dennis, J. E., Jr., and R. B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Englewood Cliffs, NJ: Prentice-Hall, 1983.