Main Content

Parametric modeling techniques find the parameters for a mathematical model describing a signal, system, or process. These techniques use known information about the system to determine the model. Applications for parametric modeling include speech and music synthesis, data compression, high-resolution spectral estimation, communications, manufacturing, and simulation.

The toolbox parametric modeling functions operate with the rational transfer function model. Given appropriate information about an unknown system (impulse or frequency response data, or input and output sequences), these functions find the coefficients of a linear system that models the system.

One important application of the parametric modeling functions is in the design of filters that have a prescribed time or frequency response.

Here is a summary of the parametric modeling functions in this toolbox.

Domain | Functions | Description |
---|---|---|

Time | Generate all-pole filter coefficients that model an input data sequence using the Levinson-Durbin algorithm. | |

Generate all-pole filter coefficients that model an input data sequence by minimizing the forward prediction error. | ||

Generate all-pole filter coefficients that model an input data sequence by minimizing the forward and backward prediction errors. | ||

Generate all-pole filter coefficients that model an input data sequence using an estimate of the autocorrelation function. | ||

Linear Predictive Coding. Generate all-pole recursive filter whose impulse response matches a given sequence. | ||

Generate IIR filter whose impulse response matches a given sequence. | ||

Find IIR filter whose output, given a specified input sequence, matches a given output sequence. | ||

Frequency | Generate digital or analog filter coefficients given complex frequency response data. |

The `lpc`

, `prony`

, and `stmcb`

functions find the coefficients of a digital rational transfer
function that approximates a given time-domain impulse response. The algorithms
differ in complexity and accuracy of the resulting model.

Linear prediction modeling assumes that each output sample of a signal,
`x(k)`

, is a linear combination of the past
`n`

outputs (that is, it can be linearly predicted from
these outputs), and that the coefficients are constant from sample to
sample:

$$x\left(k\right)=-a\left(2\right)x\left(k-1\right)-a\left(3\right)x\left(k-2\right)-\cdots -a\left(n+1\right)x\left(k-n\right).$$

An `n`

th-order all-pole model of a signal
`x`

is

a = lpc(x,n)

To illustrate `lpc`

, create a sample signal that is the
impulse response of an all-pole filter with additive white noise:

x = impz(1,[1 0.1 0.1 0.1 0.1],10) + randn(10,1)/10;

The coefficients for a fourth-order all-pole filter that models the system are

a = lpc(x,4)

`lpc`

first calls `xcorr`

to find a biased estimate of the correlation function
of `x`

, and then uses the Levinson-Durbin recursion,
implemented in the `levinson`

function, to find the
model coefficients `a`

. The Levinson-Durbin recursion is a fast
algorithm for solving a system of symmetric Toeplitz linear equations.
`lpc`

's entire algorithm for `n`

=
`4`

is

```
r = xcorr(x);
r(1:length(x)-1) = []; % Remove corr. at negative lags
a = levinson(r,4)
```

You could form the linear prediction coefficients with other assumptions by
passing a different correlation estimate to `levinson`

, such as
the biased correlation estimate:

r = xcorr(x,'biased'); r(1:length(x)-1) = []; % Remove corr. at negative lags a = levinson(r,4)

The `prony`

function models a signal
using a specified number of poles and zeros. Given a sequence
`x`

and numerator and denominator orders
`n`

and `m`

, respectively, the
statement

[b,a] = prony(x,n,m)

finds the numerator and denominator coefficients of an IIR filter whose
impulse response approximates the sequence `x`

.

The `prony`

function implements the method described in
[4] Parks and
Burrus. This method uses a variation of the covariance method of AR modeling to
find the denominator coefficients `a`

, and then finds the
numerator coefficients `b`

for which the resulting filter's
impulse response matches exactly the first `n`

+
`1`

samples of `x`

. The filter is not
necessarily stable, but it can potentially recover the coefficients exactly if
the data sequence is truly an autoregressive moving-average (ARMA) process of
the correct order.

**Note**

The functions `prony`

and `stmcb`

(described next) are more accurately described as ARX models in system
identification terminology. ARMA modeling assumes noise only at the inputs,
while ARX assumes an external input. `prony`

and
`stmcb`

know the input signal: it is an impulse for
`prony`

and is arbitrary for
`stmcb`

.

A model for the test sequence `x`

(from the earlier
`lpc`

example) using a third-order IIR filter is

[b,a] = prony(x,3,3)

The `impz`

command shows how well
this filter's impulse response matches the original sequence:

```
format long
[x impz(b,a,10)]
```

Notice that the first four samples match exactly. For an example of exact recovery, recover the coefficients of a Butterworth filter from its impulse response:

[b,a] = butter(4,.2); h = impz(b,a,26); [bb,aa] = prony(h,4,4);

Try this example; you'll see that `bb`

and
`aa`

match the original filter coefficients to within a
tolerance of 10^{-13}.

The `stmcb`

function determines the
coefficients for the system
*b*(*z*)/*a*(*z*)
given an approximate impulse response `x`

, as well as the
desired number of zeros and poles. This function identifies an unknown system
based on both input and output sequences that describe the system's behavior, or
just the impulse response of the system. In its default mode,
`stmcb`

works like `prony`

.

[b,a] = stmcb(x,3,3)

`stmcb`

also finds systems that match given input and
output sequences:

```
y = filter(1,[1 1],x); % Create an output signal.
[b,a] = stmcb(y,x,0,1)
```

In this example, `stmcb`

correctly identifies the system
used to create `y`

from `x`

.

The Steiglitz-McBride method is a fast iterative algorithm that solves for the
numerator and denominator coefficients simultaneously in an attempt to minimize
the signal error between the filter output and the given output signal. This
algorithm usually converges rapidly, but might not converge if the model order
is too large. As for `prony`

, `stmcb`

's
resulting filter is not necessarily stable due to its exact modeling
approach.

`stmcb`

provides control over several important algorithmic
parameters; modify these parameters if you are having trouble modeling the data.
To change the number of iterations from the default of five and provide an
initial estimate for the denominator coefficients:

n = 10; % Number of iterations a = lpc(x,3); % Initial estimates for denominator [b,a] = stmcb(x,3,3,n,a);

The function uses an all-pole model created with `prony`

as
an initial estimate when you do not provide one of your own.

To compare the functions `lpc`

, `prony`

,
and `stmcb`

, compute the signal error in each case:

a1 = lpc(x,3); [b2,a2] = prony(x,3,3); [b3,a3] = stmcb(x,3,3); [x-impz(1,a1,10) x-impz(b2,a2,10) x-impz(b3,a3,10)]

In comparing modeling capabilities for a given order IIR model, the last
result shows that for this example, `stmcb`

performs best,
followed by `prony`

, then `lpc`

. This
relative performance is typical of the modeling functions.

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 [h,w] = freqz(b,a,64); % Compute frequency response [b4,a4] = invfreqz(h,w,4,4) % Model: n = 4, m = 4

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
```

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 unit weighting vector [b30,a30] = invfreqz(h,w,3,3,wt,30) % 30 iterations

The resulting filter is always stable.

Graphically compare the results of the first and second algorithms to the original Butterworth filter with FVTool (and select the Magnitude and Phase Responses):

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 sum(abs(h-freqz(b30,a30,w)).^2) % Total error, algorithm 2