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.
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.
Generate digital or analog filter coefficients given complex frequency response data.
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
n outputs (that is, it can be linearly predicted
from these outputs), and that the coefficients are constant from sample
nth-order all-pole model of a signal
a = lpc(x,n)
lpc, create a sample signal
that is the impulse response of an all-pole filter with additive white
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
find a biased estimate
of the correlation function of
and then uses the Levinson-Durbin recursion,
implemented in the
levinson function, to find the model
a. The Levinson-Durbin recursion is
a fast algorithm for solving a system of symmetric Toeplitz linear
lpc's entire algorithm for
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
such as the biased correlation estimate:
r = xcorr(x,'biased'); r(1:length(x)-1) = ; % Remove corr. at negative lags a = levinson(r,4)
[b,a] = prony(x,n,m)
finds the numerator and denominator coefficients of an IIR filter
whose impulse response approximates the sequence
implements the method described in  Parks and
Burrus (pgs. 226-228). This method uses a variation
of the covariance method of AR modeling to find the denominator coefficients
a, and then finds the numerator coefficients
which the resulting filter's impulse response matches exactly the
1 samples of
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.
A model for the test sequence
x (from the
lpc example) using a third-order IIR filter
[b,a] = prony(x,3,3)
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
the original filter coefficients to within a tolerance of 10-13.
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
[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
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
resulting filter is not necessarily stable due to its exact modeling
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
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
an initial estimate when you do not provide one of your own.
To compare the functions
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,
best, followed by
This relative performance is typical of the modeling functions.
invfreqz functions implement the
inverse operations of
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
the discussion also applies to
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.
a filter of any order to fit the frequency data; a third-order example
[b4,a4] = invfreqz(h,w,3,3) % Find third-order IIR
filters with real coefficients; for a data point at positive frequency
f, the functions fit the frequency
response at both
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
\ operator. Here A(w(k))
are the Fourier transforms of the polynomials
at the frequency w(k), and n is
the number of frequency points (the length of
weights the error relative to the error at different frequencies.
includes a weighting vector. In this mode, the filter resulting
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
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):
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