Frequency-Domain Based Modeling

The invfreqs and invfreqz functions implement the inverse operations of freqs and freqz; they find an analog or digital transfer function of a specified order that matches a given complex frequency response. Though the following examples demonstrate invfreqz, the discussion also applies to invfreqs.

To recover the original filter coefficients from the frequency response of a simple digital filter

[b,a] = butter(4,0.4)    % Design Butterworth lowpass
b =
    0.0466    0.1863    0.2795    0.1863    0.0466
a =
    1.0000   -0.7821    0.6800   -0.1827    0.0301
[h,w] = freqz(b,a,64);         % Compute frequency response
[bb,aa] = invfreqz(h,w,4,4)    % Model: n = 4, m = 4
bb =
    0.0466    0.1863    0.2795    0.1863    0.0466
aa =
    1.0000   -0.7821    0.6800   -0.1827    0.0301

The vector of frequencies w has the units in rad/sample, and the frequencies need not be equally spaced. invfreqz finds a filter of any order to fit the frequency data; a third-order example is

[bb,aa] = invfreqz(h,w,3,3)   % Find third-order IIR
bb =
    0.0464    0.1785    0.2446    0.1276
aa =
    1.0000   -0.9502    0.7382   -0.2006

Both invfreqs and invfreqz design filters with real coefficients; for a data point at positive frequency f, the functions fit the frequency response at both f and -f.

By default invfreqz 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 MATLAB's \ 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). wt(k) weights the error relative to the error at different frequencies. The syntax

invfreqz(h,w,n,m,wt)

includes a weighting vector. In this mode, the filter resulting from invfreqz is not guaranteed to be stable.

invfreqz provides a superior ("output-error") algorithm that solves the direct problem of minimizing the weighted sum of the squared error between the actual frequency response points and the desired response

To use this algorithm, specify a parameter for the iteration count after the weight vector parameter

wt = ones(size(w));    % Create unity weighting vector
[bbb,aaa] = invfreqz(h,w,3,3,wt,30)  % 30 iterations
bbb =
    0.0464    0.1829    0.2572    0.1549
aaa =
    1.0000   -0.8664    0.6630   -0.1614

The resulting filter is always stable.

Graphically compare the results of the first and second algorithms to the original Butterworth filter.

[H1,w] = freqz(b,a);
[H2] = freqz(bb,aa);
[H3] = freqz(bbb,aaa);
s.plot = 'mag';
s.xunits = 'rad/sample';
s.yunits = 'linear';
freqzplot([H1 H2 H3],w);
legend('Original','First Estimate','Second Estimate',3);
grid on

 

To verify the superiority of the fit numerically, type

sum(abs(h-freqz(bb,aa,w)).^2)   % Total error, algorithm 1
ans =
    0.0200
sum(abs(h-freqz(bbb,aaa,w)).^2) % Total error, algorithm 2
ans =
    0.0096