Signal Processing Toolbox

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
  [b4,a4] = invfreqz(h,w,4,4)    % Model: n = 4, m = 4
  b4 =
      0.0466    0.1863    0.2795    0.1863    0.0466
  a4 =
      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

• [b4,a4] = invfreqz(h,w,3,3)   % Find third-order IIR
  b4 =
      0.0464    0.1785    0.2446    0.1276
  a4 =
      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 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). 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
  [b30,a30] = invfreqz(h,w,3,3,wt,30)  % 30 iterations
  b30 =
      0.0464    0.1829    0.2572    0.1549
  a30 =
      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:

• fvtool(b,a,b4,a4,b30,a30)
  
  
  

To verify the superiority of the fit numerically, type

• sum(abs(h-freqz(b4,a4,w)).^2)   % Total error, algorithm 1
  ans =
      0.0200
  sum(abs(h-freqz(b30,a30,w)).^2) % Total error, algorithm 2
  ans =
      0.0096
  

Time-Domain Based Modeling

Resampling