Copyright © The Octave Project Developers
Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies.
Permission is granted to copy and distribute modified versions of this manual under the conditions for verbatim copying, provided that the entire resulting derived work is distributed under the terms of a permission notice identical to this one.
Permission is granted to copy and distribute translations of this manual into another language, under the above conditions for modified versions.
The Signal Toolkit contains signal processing tools, including filtering, windowing and display functions.
The Signal Toolkit must be installed and then loaded to be used.
It can be installed in GNU Octave directly from octave-forge,
If running in Windows, the package may already be installed, to check run:
pkg list signal
With an internet connection available, the Signal package can be installed from octave-forge using the following command within GNU Octave:
pkg install -forge signal
The latest released version of the toolkit will be downloaded and installed.
Otherwise, if the package file has already been downloaded it can be installed using the follwoing command in GNU Octave:
pkg install signal-1.4.8.tar.gz
Regardless of the method of installing the toolkit, in order to use its functions, the toolkit must be loaded using the pkg load command:
pkg load signal
The toolkit must be loaded on each GNU Octave session.
y = buffer (x, n, p, opt) ¶[y, z, opt] = buffer (…) ¶Buffer a signal into a data frame. The arguments to buffer are
The data to be buffered.
The number of rows in the produced data buffer. This is an positive integer value and must be supplied.
An integer less than n that specifies the under- or overlap between column in the data frame. The default value of p is 0.
In the case of an overlap, opt can be either a vector of length p or the string ’nodelay’. If opt is a vector, then the first p entries in y will be filled with these values. If opt is the string ’nodelay’, then the first value of y corresponds to the first value of x.
In the case of an underlap, opt must be an integer between 0 and
-p. The represents the initial underlap of the first
column of y.
The default value for opt the vector zeros (1, p)
in the case of an overlap, or 0 otherwise.
In the case of a single output argument, y will be padded with zeros to fill the missing values in the data frame. With two output arguments z is the remaining data that has not been used in the current data frame.
Likewise, the output opt is the overlap, or underlap that might
be used for a future call to code to allow continuous buffering.
(t) ¶(t, f0) ¶(t, f0, t1) ¶(t, f0, t1, f1) ¶(t, f0, t1, f1, shape) ¶(t, f0, t1, f1, shape, phase) ¶Evaluate a chirp signal at time t. A chirp signal is a frequency swept cosine wave.
vector of times to evaluate the chirp signal
frequency at time t=0 [ 0 Hz ]
time t1 [ 1 sec ]
frequency at time t=t1 [ 100 Hz ]
shape of frequency sweep ’linear’ f(t) = (f1-f0)*(t/t1) + f0 ’quadratic’ f(t) = (f1-f0)*(t/t1)^2 + f0 ’logarithmic’ f(t) = (f1/f0)^(t/t1) * f0
phase shift at t=0
For example:
specgram (chirp ([0:0.001:5])); # default linear chirp of 0-100Hz in 1 sec specgram (chirp ([-2:0.001:15], 400, 10, 100, "quadratic")); soundsc (chirp ([0:1/8000:5], 200, 2, 500, "logarithmic"), 8000);
If you want a different sweep shape f(t), use the following:
y = cos (2 * pi * integral (f(t)) + phase);
[psi, x] = cmorwavf (lb, ub, n, fb, fc) ¶Compute the Complex Morlet wavelet.
y = diric (x,n) ¶Compute the dirichlet function.
See also: sinc, gauspuls, sawtooth.
yi = gauspuls (t) ¶yi = gauspuls (t, fc) ¶yi = gauspuls (t, fc, bw) ¶yi = gauspuls (t, fc, bw, bwr) ¶[yi, yq] = gauspuls (…) ¶[yi, yq, ye] = gauspuls (…) ¶Generate a Gaussian modulated sinusoidal pulse sampled at times t.
The input arguments are:
10^(bwr/20) times the peak amplitude.
The output arguments are:
sqrt (yi.^2 + yq.^2)).
See also: pulstran, rectpuls, tripuls.
y = pulstran (t, d, func, …) ¶y = pulstran (t, d, p) ¶y = pulstran (t, d, p, Fs) ¶y = pulstran (t, d, p, Fs, method) ¶Generate the signal y=sum(func(t+d,...)) for each d. If d is a matrix of two columns, the first column is the delay d and the second column is the amplitude a, and y=sum(a*func(t+d)) for each d,a. Clearly, func must be a function which accepts a vector of times. Any extra arguments needed for the function must be tagged on the end.
Example:
fs = 11025; # arbitrary sample rate f0 = 100; # pulse train sample rate w = 0.001; # pulse width of 1 millisecond auplot (pulstran (0:1/fs:0.1, 0:1/f0:0.1, "rectpuls", w), fs);
If instead of a function name you supply a pulse shape sampled at frequency Fs (default 1 Hz), an interpolated version of the pulse is added at each delay d. The interpolation stays within the time range of the delayed pulse. The interpolation method defaults to linear, but it can be any interpolation method accepted by the function interp1.
Example:
fs = 11025; # arbitrary sample rate f0 = 100; # pulse train sample rate w = boxcar(10); # pulse width of 1 millisecond at 10 kHz auplot (pulstran (0:1/fs:0.1, 0:1/f0:0.1, w, 10000), fs);
y = rectpuls (t) ¶y = rectpuls (t, w) ¶Generate a rectangular pulse over the interval [-w/2,w/2),
sampled at times t. This is useful with the function pulstran
for generating a series of pulses.
Example:
fs = 11025; # arbitrary sample rate f0 = 100; # pulse train sample rate w = 0.3/f0; # pulse width 3/10th the distance between pulses plot (pulstran (0:1/fs:4/f0, 0:1/f0:4/f0, "rectpuls", w));
See also: gauspuls, pulstran, tripuls.
y = sawtooth (t) ¶y = sawtooth (t, width) ¶Generates a sawtooth wave of period 2 * pi with limits +1/-1
for the elements of t.
width is a real number between 0 and 1 which specifies
the point between 0 and 2 * pi where the maximum is. The
function increases linearly from -1 to 1 in [0, 2 *
pi * width] interval, and decreases linearly from 1 to
-1 in the interval [2 * pi * width, 2 * pi].
If width is 0.5, the function generates a standard triangular wave.
If width is not specified, it takes a value of 1, which is a standard sawtooth function.
[psi, x] = shanwavf (lb, ub, n, fb, fc) ¶Compute the Complex Shannon wavelet.
[y s] = sigmoid_train (t, ranges, rc) ¶Evaluate a train of sigmoid functions at t.
The number and duration of each sigmoid is determined from ranges.
Each row of ranges represents a real interval, e.g. if sigmoid
i starts at t=0.1 and ends at t=0.5, then
ranges(i,:) = [0.1 0.5].
The input rc is an array that defines the rising and falling time
constants of each sigmoid. Its size must equal the size of ranges.
The individual sigmoids are returned in s. The combined sigmoid train
is returned in the vector y of length equal to t, and such that
Y = max (S).
Run demo sigmoid_train to some examples of the use of this function.
(x) ¶(x, n) ¶(x, n, Fs) ¶(x, n, Fs, window) ¶(x, n, Fs, window, overlap) ¶[S, f, t] = specgram (…) ¶Generate a spectrogram for the signal x. The signal is chopped into
overlapping segments of length n, and each segment is windowed and
transformed into the frequency domain using the FFT. The default segment
size is 256. If fs is given, it specifies the sampling rate of the
input signal. The argument window specifies an alternate window to
apply rather than the default of hanning (n). The argument
overlap specifies the number of samples overlap between successive
segments of the input signal. The default overlap is
length (window)/2.
If no output arguments are given, the spectrogram is displayed. Otherwise, S is the complex output of the FFT, one row per slice, f is the frequency indices corresponding to the rows of S, and t is the time indices corresponding to the columns of S.
Example:
x = chirp([0:0.001:2],0,2,500); # freq. sweep from 0-500 over 2 sec.
Fs=1000; # sampled every 0.001 sec so rate is 1 kHz
step=ceil(20*Fs/1000); # one spectral slice every 20 ms
window=ceil(100*Fs/1000); # 100 ms data window
specgram(x, 2^nextpow2(window), Fs, window, window-step);
## Speech spectrogram
[x, Fs] = auload(file_in_loadpath("sample.wav")); # audio file
step = fix(5*Fs/1000); # one spectral slice every 5 ms
window = fix(40*Fs/1000); # 40 ms data window
fftn = 2^nextpow2(window); # next highest power of 2
[S, f, t] = specgram(x, fftn, Fs, window, window-step);
S = abs(S(2:fftn*4000/Fs,:)); # magnitude in range 0<f<=4000 Hz.
S = S/max(S(:)); # normalize magnitude so that max is 0 dB.
S = max(S, 10^(-40/10)); # clip below -40 dB.
S = min(S, 10^(-3/10)); # clip above -3 dB.
imagesc (t, f, log(S)); # display in log scale
set (gca, "ydir", "normal"); # put the 'y' direction in the correct direction
The choice of window defines the time-frequency resolution. In speech for example, a wide window shows more harmonic detail while a narrow window averages over the harmonic detail and shows more formant structure. The shape of the window is not so critical so long as it goes gradually to zero on the ends.
Step size (which is window length minus overlap) controls the horizontal scale of the spectrogram. Decrease it to stretch, or increase it to compress. Increasing step size will reduce time resolution, but decreasing it will not improve it much beyond the limits imposed by the window size (you do gain a little bit, depending on the shape of your window, as the peak of the window slides over peaks in the signal energy). The range 1-5 msec is good for speech.
FFT length controls the vertical scale. Selecting an FFT length greater than the window length does not add any information to the spectrum, but it is a good way to interpolate between frequency points which can make for prettier spectrograms.
After you have generated the spectral slices, there are a number of decisions for displaying them. First the phase information is discarded and the energy normalized:
S = abs(S); S = S/max(S(:));
Then the dynamic range of the signal is chosen. Since information in speech is well above the noise floor, it makes sense to eliminate any dynamic range at the bottom end. This is done by taking the max of the magnitude and some minimum energy such as minE=-40dB. Similarly, there is not much information in the very top of the range, so clipping to a maximum energy such as maxE=-3dB makes sense:
S = max(S, 10^(minE/10)); S = min(S, 10^(maxE/10));
The frequency range of the FFT is from 0 to the Nyquist frequency of one half the sampling rate. If the signal of interest is band limited, you do not need to display the entire frequency range. In speech for example, most of the signal is below 4 kHz, so there is no reason to display up to the Nyquist frequency of 10 kHz for a 20 kHz sampling rate. In this case you will want to keep only the first 40% of the rows of the returned S and f. More generally, to display the frequency range [minF, maxF], you could use the following row index:
idx = (f >= minF & f <= maxF);
Then there is the choice of colormap. A brightness varying colormap such as copper or bone gives good shape to the ridges and valleys. A hue varying colormap such as jet or hsv gives an indication of the steepness of the slopes. The final spectrogram is displayed in log energy scale and by convention has low frequencies on the bottom of the image:
imagesc(t, f, flipud(log(S(idx,:))));
s = square (t, duty) ¶s = square (t) ¶Generate a square wave of period 2 pi with limits +1/-1.
If duty is specified, it is the percentage of time the square wave is "on". The square wave is +1 for that portion of the time.
on time * 100
duty cycle = ------------------
on time + off time
See also: cos, sawtooth, sin, tripuls.
y = tripuls (t) ¶y = tripuls (t, w) ¶y = tripuls (t, w, skew) ¶Generate a triangular pulse over the interval [-w/2,w/2),
sampled at times t. This is useful with the function pulstran
for generating a series of pulses.
skew is a value between -1 and 1, indicating the relative placement of the peak within the width. -1 indicates that the peak should be at -w/2, and 1 indicates that the peak should be at w/2. The default value is 0.
Example:
fs = 11025; # arbitrary sample rate f0 = 100; # pulse train sample rate w = 0.3/f0; # pulse width 3/10th the distance between pulses plot (pulstran (0:1/fs:4/f0, 0:1/f0:4/f0, "tripuls", w));
See also: gauspuls, pulstran, rectpuls.
[out] = unshiftdata (in, perm, shifts) ¶Reverse what is done by shiftdata.
See also: shiftdata.
y = vco (x, fc, fs) ¶y = vco (x, [fmin, fmax], fs) ¶Creates a signal that oscillates at a frequency determined by input x with a sampling frequency fs.
Inputs:
Outputs:
[pks, loc, extra] = findpeaks (data) ¶… = findpeaks (…, property, value) ¶… = findpeaks (…, "DoubleSided") ¶Finds peaks on data.
Peaks of a positive array of data are defined as local maxima. For double-sided data, they are maxima of the positive part and minima of the negative part. data is expected to be a single column vector.
The function returns the value of data at the peaks in pks. The index indicating their position is returned in loc.
The third output argument is a structure with additional information:
A structure containing the parabola fitted to each returned peak. The structure has two fields, "x" and "pp". The field "pp" contains the coefficients of the 2nd degree polynomial and "x" the extrema of the interval where it was fitted.
The estimated height of the returned peaks (in units of data).
The height at which the roots of the returned peaks were calculated (in units of data).
The abscissa values (in index units) at which the parabola fitted to each of the returned peaks realizes its width as defined below.
This function accepts property-value pair given in the list below:
Minimum peak height (non-negative scalar). Only peaks that exceed this
value will be returned. For data taking positive and negative values
use the option "DoubleSided". Default value eps.
Minimum separation between (positive integer). Peaks separated by less than this distance are considered a single peak. This distance is also used to fit a second order polynomial to the peaks to estimate their width, therefore it acts as a smoothing parameter. The neighborhood size is equal to the value of "MinPeakDistance". Default value 1.
Minimum width of peaks (non-negative scalar). The width of the peaks is estimated using a parabola fitted to the neighborhood of each peak. The width is calculated with the formula
a * width^2 = 1
where a is the concavity of the parabola.
Default value eps.
Maximum width of peaks (positive integer).
Default value Inf.
Tells the function that data takes positive and negative values. The base-line for the peaks is taken as the mean value of the function. This is equivalent as passing the absolute value of the data after removing the mean.
Run demo findpeaks to see some examples.
y = peak2peak (x) ¶y = peak2peak (x, dim) ¶Compute the difference between the maximum and minimum values in the vector x.
If x is a matrix, compute the difference for each column and return them in a row vector.
If the optional argument dim is given, operate along this dimension.
See also: max, min, peak2rms, rms, rssq.
y = peak2rms (x) ¶y = peak2rms (x, dim) ¶Compute the ratio of the largest absolute value to the root-mean-square (RMS) value of the vector x.
If x is a matrix, compute the peak-magnitude-to-RMS ratio for each column and return them in a row vector.
If the optional argument dim is given, operate along this dimension.
See also: max, min, peak2peak, rms, rssq.
y = rssq (x) ¶y = rssq (x, dim) ¶Compute the root-sum-of-squares (RSS) of the vector x.
The root-sum-of-squares is defined as
rssq (x) = SQRT (SUM_i x(i)^2)
If x is a matrix, compute the root-sum-of-squares for each column and return them in a row vector.
If the optional argument dim is given, operate along this dimension.
See also: mean, meansq, sumsq, rms.
c = cconv (a, b, n) ¶c = cconv (a, b) ¶Compute the modulo-N circular convolution.
a and b are input vectors and c is the modolo-n
convolution of a and b. If n is not provided,
its assumed default value is length(a) + length(b) - 1,
which provides the same result as a linear convolution.
Examples:
cconv (1:2, 1:4)
⇒ 1 4 7 10 8
cconv (1:2, 1:4, 2)
⇒ 16 14
cconv (1:2, 1:4, 4)
⇒ 9 4 7 10
See also: conv, circshift.
(a, n) ¶If a is a column vector and x is a column vector of length n, then
convmtx(a, n) * x
gives the convolution of of a and x and is the
same as conv(a, x). The difference is if
many vectors are to be convolved with the same vector, then
this technique is possibly faster.
Similarly, if a is a row vector and x is a row vector of length n, then
x * convmtx(a, n)
is the same as conv(x, a).
See also: conv.
H = corrmtx (x, m) ¶H = corrmtx (x, m, method) ¶[H, R] = corrmtx (…) ¶Build a data matrix for autocorrelation matrix estimation.
Given a vector x of length N and a model order m, compute the rectangular Toeplitz matrix H such that H’*H is a biased estimate of the autocorrelation matrix. The size of H depends on the selected method:
Uses both prewindowed and postwindowed data.
Uses prewindowed data only.
Uses postwindowed data only.
Uses nonwindowed data.
Uses forward and backward prediction error estimates.
The optional second output R is the biased autocorrelation matrix estimate H’*H.
See also: aryule.
[R, lag] = xcorr ( X ) ¶… = xcorr ( X, Y ) ¶… = xcorr ( …, maxlag) ¶… = xcorr ( …, scale) ¶Estimates the cross-correlation.
Estimate the cross correlation R_xy(k) of vector arguments X and Y or, if Y is omitted, estimate autocorrelation R_xx(k) of vector X, for a range of lags k specified by argument "maxlag". If X is a matrix, each column of X is correlated with itself and every other column.
The cross-correlation estimate between vectors "x" and "y" (of length N) for lag "k" is given by
N
R_xy(k) = sum x_{i+k} conj(y_i),
i=1
where data not provided (for example x(-1), y(N+1)) is zero. Note the
definition of cross-correlation given above. To compute a
cross-correlation consistent with the field of statistics, see xcov.
ARGUMENTS
[non-empty; real or complex; vector or matrix] data
[real or complex vector] data
If X is a matrix (not a vector), Y must be omitted. Y may be omitted if X is a vector; in this case xcorr estimates the autocorrelation of X.
[integer scalar] maximum correlation lag If omitted, the default value is N-1, where N is the greater of the lengths of X and Y or, if X is a matrix, the number of rows in X.
[character string] specifies the type of scaling applied to the correlation vector (or matrix). is one of:
return the unscaled correlation, R,
return the biased average, R/N,
return the unbiased average, R(k)/(N-|k|),
return the correlation coefficient, R/(rms(x).rms(y)), where "k" is the lag, and "N" is the length of X. If omitted, the default value is "none". If Y is supplied but does not have the same length as X, scale must be "none".
RETURNED VARIABLES
array of correlation estimates
row vector of correlation lags [-maxlag:maxlag]
The array of correlation estimates has one of the following forms: (1) Cross-correlation estimate if X and Y are vectors.
(2) Autocorrelation estimate if is a vector and Y is omitted.
(3) If X is a matrix, R is an matrix containing the cross-correlation estimate of each column with every other column. Lag varies with the first index so that R has 2*maxlag+1 rows and P^2 columns where P is the number of columns in X.
If Rij(k) is the correlation between columns i and j of X
R(k+maxlag+1,P*(i-1)+j) == Rij(k)
for lag k in [-maxlag:maxlag], or
R(:,P*(i-1)+j) == xcorr(X(:,i),X(:,j)).
reshape(R(k,:),P,P) is the cross-correlation matrix for X(k,:).
See also: xcov.
(a) ¶(a, b) ¶(…, scale) ¶Compute the 2D cross-correlation of matrices a and b.
If b is not specified, computes autocorrelation of a, i.e.,
same as xcorr (a, a).
The optional argument scale, defines the type of scaling applied to the cross-correlation matrix. Possible values are:
No scaling.
Scales the raw cross-correlation by the maximum number of elements of a and b involved in the generation of any element of c.
Scales the raw correlation by dividing each element in the cross-correlation matrix by the number of products a and b used to generate that element.
Scales the normalized cross-correlation on the range of [0 1] so that a value of 1 corresponds to a correlation coefficient of 1.
See also: conv2, corr2, xcorr.
[R, lag] = xcov ( X ) ¶… = xcov ( X, Y ) ¶… = xcov ( …, maxlag) ¶… = xcov ( …, scale) ¶Compute covariance at various lags [=correlation(x-mean(x),y-mean(y))].
input vector
if specified, compute cross-covariance between X and Y, otherwise compute autocovariance of X.
is specified, use lag range [-maxlag:maxlag], otherwise use range [-n+1:n-1].
for covariance=raw/N,
for covariance=raw/(N-|lag|),
for covariance=raw/(covariance at lag 0),
for covariance=raw
is the default.
Returns the covariance for each lag in the range, plus an optional vector of lags.
See also: xcorr.
y = filtfilt (b, a, x) ¶Forward and reverse filter the signal. This corrects for phase distortion introduced by a one-pass filter, though it does square the magnitude response in the process. That’s the theory at least. In practice the phase correction is not perfect, and magnitude response is distorted, particularly in the stop band.
Example
[b, a]=butter(3, 0.1); # 5 Hz low-pass filter t = 0:0.01:1.0; # 1 second sample x=sin(2*pi*t*2.3)+0.25*randn(size(t)); # 2.3 Hz sinusoid+noise y = filtfilt(b,a,x); z = filter(b,a,x); # apply filter plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;')
zf = filtic (b, a, y) ¶zf = filtic (b, a, y, x) ¶Set initial condition vector for filter function The vector zf has the same values that would be obtained from function filter given past inputs x and outputs y
The vectors x and y contain the most recent inputs and outputs respectively, with the newest values first:
x = [x(-1) x(-2) ... x(-nb)], nb = length(b)-1 y = [y(-1) y(-2) ... y(-na)], na = length(a)-a
If length(x)<nb then it is zero padded If length(y)<na then it is zero padded
zf = filtic(b, a, y) Initial conditions for filter with coefficients a and b and output vector y, assuming input vector x is zero
zf = filtic(b, a, y, x) Initial conditions for filter with coefficients a and b input vector x and output vector y
y = medfilt1 (x, n) ¶y = medfilt1 (x, n, [], dim) ¶y = medfilt1 (..., NaN_flag, padding) ¶Apply a one dimensional median filter with a window size of n to the data x, which must be real, double and full. For n = 2m+1, y(i) is the median of x(i-m:i+m). For n = 2m, y(i) is the median of x(i-m:i+m-1).
The calculation is performed over the first non-singleton dimension, or over dimension dim if that is specified as the fourth argument. (The third argument is ignored; Matlab used to use it to tune its algorithm.)
NaN_flag may be omitnan or includenan (the default).
If it is omitnan then any NaN values are removed from the window
before the median is taken.
Otherwise, any window containing an NaN returns a median of NaN.
padding determines how the partial windows at the start and end of
x are treated.
It may be truncate or zeropad (the default).
If it is truncate then the window for y(i) is
the intersection of the window stated above with 1:length(x).
If it is zeropad, then partial windows have additional zeros
to bring them up to size n.
See also: filter, medfilt2.
[rmsx,w] = movingrms (x,w,rc,Fs=1) ¶Calculate moving RMS value of the signal in x.
The signal is convoluted against a sigmoid window of width w and risetime rc. The units of these parameters are relative to the value of the sampling frequency given in Fs (Default value = 1).
Run demo movingrms to see an example.
See also: sigmoid_train.
y = sgolayfilt (x) ¶y = sgolayfilt (x, p) ¶y = sgolayfilt (x, p, n) ¶y = sgolayfilt (x, p, n, m) ¶y = sgolayfilt (x, p, n, m, ts) ¶y = sgolayfilt (x, p, n, m, ts) ¶y = sgolayfilt (x, f) ¶Smooth the data in x with a Savitsky-Golay smoothing filter of polynomial order p and length n, n odd, n > p. By default, p=3 and n=p+2 or n=p+3 if p is even.
If f is given as a matrix, it is expected to be a filter as
computed by sgolay.
These filters are particularly good at preserving lineshape while removing high frequency squiggles. Particularly, compare a 5 sample averager, an order 5 butterworth lowpass filter (cutoff 1/3) and sgolayfilt(x, 3, 5), the best cubic estimated from 5 points:
[b, a] = butter (5, 1/3);
x = [zeros(1,15), 10*ones(1,10), zeros(1,15)];
plot (sgolayfilt (x), "r;sgolayfilt;", ...
filtfilt (ones (1,5)/5, 1, x), "g;5 sample average;", ...
filtfilt (b, a, x), "c;order 5 butterworth;", ...
x, "+b;original data;");
See also: sgolay.
y = sosfilt (sos, x) ¶Second order section IIR filtering of x. The second order section filter is described by the matrix sos with:
| [ B1 A1 ] | |
| sos = | [ … ], |
| [ BN AN ] |
where B1 = [b0 b1 b2] and A1 = [1 a1 a2] for
section 1, etc. The b0 entry must be nonzero for each section.
L = filternorm (b, a) ¶L = filternorm (b, a, pnorm) ¶L = filternorm (b, a, 2, tol) ¶Compute the 2-norm of a digital filter defined by the numerator coefficients, b, and the denominator coefficients, a. It is also possible to compute the infinity-norm by passing inf in the pnorm parameter. pnorm only accepts 2 or inf.
Example:
[b, a] = butter (8, 0.5); filternorm (b, a)
n = filtord (b, a) ¶n = filtord (sos) ¶Returns the filter order n for a filter defined by the numerator coefficients, b, and the denominator coefficients, a. It also accepts a filter defined by a matrix of second-order sections, sos.
Example:
[b, a] = butter (8, 0.5); filtord (b, a)
h = freqs (b, a, w) ¶[h, wout] = freqs (b, a, n) ¶(…) ¶Compute the s-plane frequency response of the analog filter B(s)/A(s).
The frequency response is evaluated at the angular frequencies specified by vector w (in rad/s). If the third argument is a scalar integer n, or if it is omitted, the frequency response is computed at n logarithmically spaced frequencies (default n = 200).
If no output argument is requested, the magnitude and phase are plotted.
Example:
b = [1 2]; a = [1 1]; w = linspace (0, 4, 128); freqs (b, a, w);
f = freqspace (N) ¶f = freqspace (N, "whole") ¶[f1, f2] = freqspace (n) ¶[f1, f2] = freqspace ([m, n]) ¶[x, y] = freqspace (…, "meshgrid") ¶Generate frequency spacing for frequency responses.
For one-dimensional usage, freqspace (N) returns N
equally-spaced points along the upper half of the unit circle, with
normalized frequency in [0, 1]. When "whole" is specified, returns points
along the entire unit circle, ranging from [0, 2).
For two-dimensional usage, [f1, f2] = freqspace (n)
returns two frequency vectors for an n-n matrix;
[f1, f2] = freqspace ([m, n]) returns
vectors for an m-n matrix. The 2-D frequency vectors are
equally spaced between -1 and 1, where 1 corresponds to the Nyquist
frequency.
When "meshgrid" is specified, it is equivalent to calling freqspace
followed by meshgrid on the results.
f = fwhm (y) ¶f = fwhm (x, y) ¶f = fwhm (…, "zero") ¶f = fwhm (…, "min") ¶f = fwhm (…, "alevel", level) ¶f = fwhm (…, "rlevel", level) ¶Compute peak full-width at half maximum (FWHM) or at another level of peak maximum for vector or matrix data y, optionally sampled as y(x). If y is a matrix, return FWHM for each column as a row vector.
The default option "zero" computes fwhm at half maximum, i.e. 0.5*max(y). The option "min" computes fwhm at the middle curve, i.e. 0.5*(min(y)+max(y)).
The option "rlevel" computes full-width at the given relative level of peak
profile, i.e. at rlevel*max(y) or rlevel*(min(y)+max(y)),
respectively. For example, fwhm (…, "rlevel", 0.1) computes
full width at 10 % of peak maximum with respect to zero or minimum; FWHM is
equivalent to fwhm(…, "rlevel", 0.5).
The option "alevel" computes full-width at the given absolute level of y.
Return 0 if FWHM does not exist (e.g. monotonous function or the function does not cut horizontal line at rlevel*max(y) or rlevel*(max(y)+min(y)) or alevel, respectively).
[g, w] = grpdelay (b) ¶[g, w] = grpdelay (b, a) ¶[g, w] = grpdelay (…, n) ¶[g, w] = grpdelay (…, n, "whole") ¶[g, f] = grpdelay (…, n, Fs) ¶[g, f] = grpdelay (…, n, "whole", Fs) ¶[g, w] = grpdelay (…, w) ¶[g, f] = grpdelay (…, f, Fs) ¶(…) ¶Compute the group delay of a filter.
[g, w] = grpdelay(b) returns the group delay g of the FIR filter with coefficients b. The response is evaluated at 512 angular frequencies between 0 and pi. w is a vector containing the 512 frequencies. The group delay is in units of samples. It can be converted to seconds by multiplying by the sampling period (or dividing by the sampling rate fs).
[g, w] = grpdelay(b,a) returns the group delay of the rational IIR filter whose numerator has coefficients b and denominator coefficients a.
[g, w] = grpdelay(b,a,n) returns the group delay evaluated at n angular frequencies. For fastest computation n should factor into a small number of small primes.
[g, w] = grpdelay(b,a,n,’whole’) evaluates the group delay at n frequencies between 0 and 2*pi.
[g, f] = grpdelay(b,a,n,Fs) evaluates the group delay at n frequencies between 0 and Fs/2.
[g, f] = grpdelay(b,a,n,’whole’,Fs) evaluates the group delay at n frequencies between 0 and Fs.
[g, w] = grpdelay(b,a,w) evaluates the group delay at frequencies w (radians per sample).
[g, f] = grpdelay(b,a,f,Fs) evaluates the group delay at frequencies f (in Hz).
grpdelay(...) plots the group delay vs. frequency.
If the denominator of the computation becomes too small, the group delay is set to zero. (The group delay approaches infinity when there are poles or zeros very close to the unit circle in the z plane.)
Theory: group delay, g(w) = -d/dw [arg{H(e^jw)}], is the rate of change of phase with respect to frequency. It can be computed as:
d/dw H(e^-jw)
g(w) = -------------
H(e^-jw)
where
H(z) = B(z)/A(z) = sum(b_k z^k)/sum(a_k z^k).
By the quotient rule,
A(z) d/dw B(z) - B(z) d/dw A(z)
d/dw H(z) = -------------------------------
A(z) A(z)
Substituting into the expression above yields:
A dB - B dA
g(w) = ----------- = dB/B - dA/A
A B
Note that,
d/dw B(e^-jw) = sum(k b_k e^-jwk)
d/dw A(e^-jw) = sum(k a_k e^-jwk)
which is just the FFT of the coefficients multiplied by a ramp.
As a further optimization when nfft>>length(a), the IIR filter (b,a) is converted to the FIR filter conv(b,fliplr(conj(a))). For further details, see http://ccrma.stanford.edu/~jos/filters/Numerical_Computation_Group_Delay.html
[x, t] = impz (b) ¶[x, t] = impz (b, a) ¶[x, t] = impz (b, a, n) ¶[x, t] = impz (b, a, n, fs) ¶(…) ¶Generate impulse-response characteristics of the filter.
The filter coefficients correspond to the z-plane rational function
with numerator b and denominator a. If a is not specified, it
defaults to 1. When n is a scalar, it specifies the number of
points to compute (default: determined by impzlength). When
n is a vector of non-negative integers, the response is
computed only at those sample indices. The sampling frequency
fs (default: 1) controls the time spacing in the output
t. With no output arguments, the result is plotted.
See also: freqz, zplane, stepz,impzlength.
len = impzlength (b) ¶len = impzlength (b, a) ¶len = impzlength (sos) ¶len = impzlength (…, tol) ¶Return the impulse response length of the specified filter.
For a finite impulse response (FIR) filter specified by the numerator coefficients b, the length is simply the number of coefficients in b.
For an infinite impulse response (IIR) filter specified by the numerator b and denominator a polynomials in z^-1, the function computes an effective impulse response sequence length.
The filter can also be specified by a K-by-6 second-order sections matrix sos, where K is the number of sections. In this case, the matrix is converted to transfer function form b and a before computing the length.
The algorithm proceeds as follows:
floor (M * log10 (tol) / log10 (maxpole)) + delay
where M is the multiplicity of the dominant pole and d is the initial delay (number of leading zeros in b).
floor (6 / log10 (maxpole))
The optional argument tol specifies the tolerance used to estimate the effective length of an IIR filter’s impulse response. The default tolerance is 5e-5. Increasing tol estimates a shorter effective length, while decreasing tol produces a longer effective length.
The returned value len is the effective impulse response length
of the specified filter. This function is used by impz and
stepz to determine the number of points to plot.
See also: impz, stepz.
L = isallpass (b, a) ¶L = isallpass (sos) ¶Determine whether a digital filter is allpass. The filter might be defined by the numerator coefficients, b, and the denominator coefficients, a, or, alternatively, by a matrix of second-order sections, sos.
Example:
a = [1 2 3]; b = [3 2 1]; isallpass (b, a)
Ref [1] Shyu, Jong-Jy, & Pei, Soo-Chang, A new approach to the design of complex all-pass IIR digital filters, Signal Processing, 40(2–3), 207–215, 1994. https://doi.org/10.1016/0165-1684(94)90068-x
Ref [2] Vaidyanathan, P. P. Multirate Systems and Filter Banks. 1st edition, Pearson College Div, 1992.
L = ismaxphase (b, a) ¶L = ismaxphase (sos) ¶L = ismaxphase (…, tol) ¶Determine whether a digital filter is maximum phase (maximum energy-delay). The filter might be defined by the numerator coefficients, b, and the denominator coefficients, a, or, alternatively, by a matrix of second-order sections, sos. A tolerance tol might be given to define when two numbers are close enough to be considered equal.
Example:
b = [1 2 4 4 2 1]; zplane (b); ismaxphase (b)
Ref [1] Oppenheim, Alan, and Ronald Schafer. Discrete-Time Signal Processing. 3rd edition, Pearson, 2009.
L = isminphase (b, a) ¶L = isminphase (sos) ¶L = isminphase (…, tol) ¶Determine whether a digital filter is minimum phase. The filter might be defined by the numerator coefficients, b, and the denominator coefficients, a, or, alternatively, by a matrix of second-order sections, sos. A toleranve tol might be given to define when two numbers are close enough to be considered equal.
Example:
a = [1 0.5]; b = [3 1]; isminphase (b, a)
Ref [1] Oppenheim, Alan, and Ronald Schafer. Discrete-Time Signal Processing. 3rd edition, Pearson, 2009.
FLAG = isstable (B, A) ¶FLAG = isstable (sos) ¶Returns a logical output equal to TRUE, if the filter is stable. This can be done with coeffients of the filer B and A. Alternatively by using a second order sections matrix (SOS).
Inputs:
Output:
Examples:
b = [1 2 3 4 5 5 1 2]; a = [4 5 6 7 9 10 4 6]; flag = isstable (b, a) flag = 0
Using SOS
[z, p, k] = butter (6, 0.7, 'high'); sos = zp2sos (z, p, k); flag = isstable (sos) flag = 1
[phi, w] = phasez (b, a, n) ¶[phi, w] = phasez (b, a) ¶[phi, w] = phasez (sos, n) ¶[phi, w] = phasez (sos) ¶[phi, w] = phasez (…, n, "whole") ¶[phi, w] = phasez (…, n, Fs) ¶(…) ¶Compute the phase response of digital filter defined either by its coefficients (b and a are the numerator and denominator coefficients respectively) or by its second-order sections representation, given by the matrix sos. The output phi is the phase response computed in a vector the vector of frequencies w.
The phase response is evaluated at n angular frequencies between 0 and pi.
If a is omitted, the denominator is assumed to be 1 (this corresponds to a simple FIR filter).
If n is omitted, a value of 512 is assumed.
If the third/forth argument, "whole", is given, the response is
evaluated at n angular frequencies between 0 and
2*pi.
It is possible also to pass the value "half", which will lead to
the default behaviour.
Example:
[b, a] = butter (2, [.15,.3]); phasez (b, a);
Ref [1] Oppenheim, Alan, and Ronald Schafer. Discrete-Time Signal Processing. 3rd edition, Pearson, 2009.
See also: freqz, phasedelay.
[x, t] = stepz (b) ¶[x, t] = stepz (b, a) ¶[x, t] = stepz (b, a, n) ¶[x, t] = stepz (b, a, n, fs) ¶(…) ¶Generate step-response characteristics of the filter.
The filter coefficients correspond to the z-plane rational function
with numerator b and denominator a. If a is not specified, it
defaults to 1. When n is a scalar, it specifies the number of
points to compute (default: determined by impzlength). When
n is a vector of non-negative integers, the response is
computed only at those sample indices. The sampling frequency
fs (default: 1) controls the time spacing in the output
t. With no output arguments, the result is plotted.
See also: freqz, zplane, impz, impzlength.
(z, p) ¶(b, a) ¶[hz, hp, ht] = zplane (___) ¶Plot the poles and zeros on a complex plane. If the arguments are column vectors z and p, the complex zeros z and poles p are displayed. If the arguments are row vectors b and a, the zeros and poles of the transfer function represented by these filter coefficients are displayed.
If z and p are matrices, the columns are distinct sets of zeros and poles and are displayed together in distinct colors.
Note that due to the nature of the roots function, poles and zeros
may be displayed as occurring around a circle rather than at a single
point.
The transfer function is
B(z) b0 + b1 z^(-1) + b2 z^(-2) + ... + bM z^(-M)
H(z) = ---- = --------------------------------------------
A(z) a0 + a1 z^(-1) + a2 z^(-2) + ... + aN z^(-N)
b0 (z - z1) (z - z2) ... (z - zM)
= -- z^(-M+N) ------------------------------
a0 (z - p1) (z - p2) ... (z - pN)
If called with only one argument, the poles p defaults to an empty vector, and the denominator coefficient vector a defaults to 1.
If output variables are provided, hz is the handle to the zero lines, hp is the handle to the pole lines of the pole-zero plot. ht is a vector of handles to the axes/unit circle line and to text objects. If there are no zeros or no poles, hz or hp is the empty matrix, [].
sos = cell2sos (cll) ¶[sos, g] = cell2sos (cll) ¶Convert a second-order-section cell array to matrix form.
Given a cell array cll representing second-order sections, return the equivalent matrix sos and an optional overall gain g.
A valid cll input is a cell array of 2-element cell arrays.
Each element consists of a numerator vector [b0, b1, b2] and a
denominator vector [1, a1, a2]. If the first element contains
two scalars {gn, gd}, it is treated as an overall gain.
Examples:
## Two sections, no gain:
cll = {{[2 4 2] [6 0 2]}, {[3 3 0] [6 0 0]}};
sos = cell2sos (cll)
⇒ sos =
2 4 2 6 0 2
3 3 0 6 0 0
## With gain:
cll = {{7 5}, {[2 4 2] [6 0 2]}, {[3 3 0] [6 0 0]}};
[sos, g] = cell2sos (cll)
⇒ sos =
2 4 2 6 0 2
3 3 0 6 0 0
⇒ g = 1.4000
See also: sos2cell.
[b, a] = eqtflength (num, den) ¶[b, a, n, m] = eqtflength (num, den) ¶Equalize numerator and denominator polynomial lengths.
Pads the shorter polynomial with trailing zeros so that the returned numerator b and denominator a have equal length while representing the same discrete-time transfer function.
Outputs n and m are the numerator and denominator orders, respectively, excluding trailing zeros.
Example:
[b, a] = eqtflength ([1 2], [1 0.5 0.25])
b = [1 2 0]
a = [1 0.5 0.25]
See also: tf2ss, tf2zp.
[r, p, f, m] = residued (b, a) ¶Compute the partial fraction expansion (PFE) of filter
H(z) = B(z)/A(z). In the usual PFE function residuez, the
IIR part (poles p and residues r) is driven in parallel
with the FIR part (f). In this variant, the IIR part is driven by
the output of the FIR part. This structure can be more accurate in
signal modeling applications.
INPUTS:
b and a are vectors specifying the digital filter
H(z) = B(z)/A(z). See help filter for documentation of the
b and a filter coefficients.
RETURNED:
EXAMPLES:
See test residued verbose to see a number of examples.
For the theory of operation, see
‘http://ccrma.stanford.edu/~jos/filters/residued.html’
See also: residue, residued.
[r, p, f, m] = residuez (b, a) ¶Compute the partial fraction expansion of filter H(z) = B(z)/A(z).
INPUTS:
b and a are vectors specifying the digital filter
H(z) = B(z)/A(z). See help filter for documentation of the
b and a filter coefficients.
RETURNED:
EXAMPLES:
See test residuez verbose to see a number of examples.
For the theory of operation, see
‘http://ccrma.stanford.edu/~jos/filters/residuez.html’
See also: residue, residued.
C = sos2cell (S) ¶C = sos2cell (S, G) ¶Convert a second-order-section matrix to a cell array.
S is an L-by-6 matrix, where each row represents a second-order section in the form:
S = [B1 A1;
B2 A2;
...
BL AL]
where Bi and Ai are the numerator and denominator coefficients of a linear or quadratic polynomial. The function converts this matrix into a cell array C with the following format:
C = { {B1, A1}, {B2, A2}, ..., {BL, AL} }
Each element of the cell array is a cell containing a pair of vectors: Bi and Ai.
If an additional gain term G is provided, the function returns:
C = { {G, 1}, {B1, A1}, {B2, A2}, ..., {BL, AL} }
where {G, 1} represents the constant gain term applied to the filter.
S = [ [1, 2, 3, 4, 5, 6];
[7, 8, 9, 10, 11, 12] ];
C = sos2cell(S);
[a, b, c, d] = sos2ss (sos) ¶Convert series second-order sections to state-space.
See also: sos2ss, ss2tf.
[b, a] = sos2tf (sos) ¶[b, a] = sos2tf (sos, g) ¶Convert series second-order sections to transfer function.
INPUTS:
sos = [B1.' A1.'; ...; BN.' AN.']
where
B1.' = [b0 b1 b2] and A1.' = [a0 a1 a2] for
section 1, etc.
a0 is usually equal to 1 because all 2nd order transfer functions
can be scaled so that a0 = 1.
However, this is not mandatory for this implementation, which supports
all kinds of transfer functions, including first order transfer functions.
See filter for documentation of the second-order direct-form filter
coefficients Bi and Ai.
RETURNED:
b and a are vectors specifying the analog or digital filter
H(s) = B(s)/A(s) or H(z) = B(z)/A(z).
See filter for further details.
See also: tf2sos, zp2sos, sos2pz, zp2tf, tf2zp.
[z, p, k] = sos2zp (sos) ¶[z, p, k] = sos2zp (sos, g) ¶Convert series second-order sections to zeros, poles, and gains (pole residues).
INPUTS:
sos = [B1.' A1.'; ...; BN.' AN.']
where
B1.' = [b0 b1 b2] and A1.' = [a0 a1 a2] for
section 1, etc.
a0 is usually equal to 1 because all 2nd order transfer functions can
be scaled so that a0 = 1.
However, this is not mandatory for this implementation, which supports
all kinds of transfer functions, including first order transfer functions.
See filter for documentation of the second-order direct-form filter
coefficients Bi and Ai.
RETURNED:
EXAMPLE:
[z, p, k] = sos2zp ([1 0 1, 1 0 -0.81; 1 0 0, 1 0 0.49])
⇒ z =
0 + 1i
0 - 1i
0 + 0i
0 + 0i
⇒ p =
-0.9000 + 0i
0.9000 + 0i
0 + 0.7000i
0 - 0.7000i
⇒ k = 1
See also: zp2sos, sos2tf, tf2sos, zp2tf, tf2zp.
[num, den] = ss2tf (a, b, c, d) ¶Conversion from state-space to transfer function representation. The state space system:
.
x = Ax + Bu
y = Cx + Du
is converted to a transfer function:
num(s)
G(s)=-------
den(s)
[z, p, k] = ss2zp (a, b, c, d) ¶Converts a state space representation to a set of poles and zeros; k is a gain associated with the zeros.
[sos, g] = tf2sos (b, a) ¶sos = tf2sos (b, a) ¶Convert direct-form filter coefficients to series second-order sections.
INPUTS:
b and a are vectors specifying the digital filter
H(z) = B(z)/A(z). See filter for documentation of the b
and a filter coefficients.
RETURNED:
sos = [b1.' a1.'; ...; bn.' an.']
where
B1.' = [b0 b1 b2] and A1.' = [1 a1 a2] for
section 1, etc. The b0 entry must be nonzero for each section (zeros at
infinity not supported).
If called with only one output argument, the overall filter gain is applied to the first second-order section in the matrix sos.
EXAMPLE:
B = [1 0 0 0 0 1];
A = [1 0 0 0 0 .9];
[sos, g] = tf2sos (B, A)
sos =
1.00000 0.61803 1.00000 1.00000 0.60515 0.95873
1.00000 -1.61803 1.00000 1.00000 -1.58430 0.95873
1.00000 1.00000 -0.00000 1.00000 0.97915 -0.00000
g = 1
See also: sos2tf, zp2sos, sos2pz, zp2tf, tf2zp.
[a, b, c, d] = tf2ss (num, den) ¶Conversion from transfer function to state-space. The state space system:
.
x = Ax + Bu
y = Cx + Du
is obtained from a transfer function:
num(s)
G(s)=-------
den(s)
The state space system matrices obtained from this function will be in observable companion form as Wolovich’s Observable Structure Theorem is used.
[z, p, k] = tf2zp (num, den) ¶Convert transfer functions to poles-and-zero representations.
Returns the zeros and poles of the system defined by num/den. k is a gain associated with the system zeros.
[sos, g] = zp2sos (z) ¶[sos, g] = zp2sos (z, p) ¶[sos, g] = zp2sos (z, p, k) ¶sos = zp2sos (…) ¶Convert filter poles and zeros to second-order sections.
INPUTS:
RETURNED:
sos = [B1.' A1.'; ...; BN.' AN.']
where
B1.' = [b0 b1 b2] and A1.' = [a0 a1 a2] for
section 1, etc.
See filter for documentation of the second-order direct-form filter
coefficients Bi and Ai, i=1:N.
If called with only one output argument, the overall filter gain is applied to the first second-order section in the matrix sos.
EXAMPLE:
[z, p, k] = tf2zp ([1 0 0 0 0 1], [1 0 0 0 0 .9]);
[sos, g] = zp2sos (z, p, k)
sos =
1.0000 0.6180 1.0000 1.0000 0.6051 0.9587
1.0000 -1.6180 1.0000 1.0000 -1.5843 0.9587
1.0000 1.0000 0 1.0000 0.9791 0
g =
1
See also: sos2zp, sos2tf, tf2sos, zp2tf, tf2zp.
[a, b, c, d] = zp2ss (z, p, k) ¶Conversion from zero / pole to state space.
Inputs
Vectors of (possibly) complex poles and zeros of a transfer function. Complex values must come in conjugate pairs (i.e., x+jy in z means that x-jy is also in z).
Real scalar (leading coefficient).
Outputs
The state space system, in the form:
.
x = Ax + Bu
y = Cx + Du
[zero, pole, gain] = besselap (n) ¶Return bessel analog filter prototype.
References:
http://en.wikipedia.org/wiki/Bessel_polynomials
[b, a] = besself (n, w) ¶[b, a] = besself (n, w, "high") ¶[z, p, g] = besself (…) ¶[a, b, c, d] = besself (…) ¶[…] = besself (…, "z") ¶Generate a Bessel filter. Default is a Laplace space (s) filter.
[b,a] = besself(n, Wc) low pass filter with cutoff pi*Wc radians
[b,a] = besself(n, Wc, ’high’) high pass filter with cutoff pi*Wc radians
[z,p,g] = besself(...) return filter as zero-pole-gain rather than coefficients of the numerator and denominator polynomials.
[...] = besself(...,’z’) return a discrete space (Z) filter, W must be less than 1.
[a,b,c,d] = besself(...) return state-space matrices
References:
Proakis & Manolakis (1992). Digital Signal Processing. New York: Macmillan Publishing Company.
[Zb, Za] = bilinear (Sb, Sa, T) ¶[Zb, Za] = bilinear (Sz, Sp, Sg, T) ¶[Zz, Zp, Zg] = bilinear (…) ¶Transform a s-plane filter specification into a z-plane specification. Filters can be specified in either zero-pole-gain or transfer function form. The input form does not have to match the output form. 1/T is the sampling frequency represented in the z plane.
Note: this differs from the bilinear function in the signal processing toolbox, which uses 1/T rather than T.
Theory: Given a piecewise flat filter design, you can transform it from the s-plane to the z-plane while maintaining the band edges by means of the bilinear transform. This maps the left hand side of the s-plane into the interior of the unit circle. The mapping is highly non-linear, so you must design your filter with band edges in the s-plane positioned at 2/T tan(w*T/2) so that they will be positioned at w after the bilinear transform is complete.
The following table summarizes the transformation:
+---------------+-----------------------+----------------------+ | Transform | Zero at x | Pole at x | | H(S) | H(S) = S-x | H(S)=1/(S-x) | +---------------+-----------------------+----------------------+ | 2 z-1 | zero: (2+xT)/(2-xT) | zero: -1 | | S -> - --- | pole: -1 | pole: (2+xT)/(2-xT) | | T z+1 | gain: (2-xT)/T | gain: (2-xT)/T | +---------------+-----------------------+----------------------+
With tedious algebra, you can derive the above formulae yourself by substituting the transform for S into H(S)=S-x for a zero at x or H(S)=1/(S-x) for a pole at x, and converting the result into the form:
H(Z)=g prod(Z-Xi)/prod(Z-Xj)
Please note that a pole and a zero at the same place exactly cancel. This is significant since the bilinear transform creates numerous extra poles and zeros, most of which cancel. Those which do not cancel have a "fill-in" effect, extending the shorter of the sets to have the same number of as the longer of the sets of poles and zeros (or at least split the difference in the case of the band pass filter). There may be other opportunistic cancellations but I will not check for them.
Also note that any pole on the unit circle or beyond will result in an unstable filter. Because of cancellation, this will only happen if the number of poles is smaller than the number of zeros. The analytic design methods all yield more poles than zeros, so this will not be a problem.
References:
Proakis & Manolakis (1992). Digital Signal Processing. New York: Macmillan Publishing Company.
[z, p, g] = buttap (n) ¶Design lowpass analog Butterworth filter.
This function exists for MATLAB compatibility only, and is equivalent
to butter (n, 1, "s").
See also: butter.
[b, a] = butter (n, wc) ¶[b, a] = butter (n, wc, filter_type) ¶[z, p, g] = butter (…) ¶[a, b, c, d] = butter (…) ¶[…] = butter (…, "s") ¶Generate a Butterworth filter. Default is a discrete space (Z) filter.
The cutoff frequency, wc should be specified in radians for
analog filters. For digital filters, it must be a value between zero
and one. For bandpass filters, wc is a two-element vector
with w(1) < w(2).
The filter type must be one of "low", "high",
"bandpass", or "stop". The default is "low"
if wc is a scalar and "bandpass" if wc is a
two-element vector.
If the final input argument is "s" design an analog Laplace
space filter.
Low pass filter with cutoff pi*Wc radians:
[b, a] = butter (n, Wc)
High pass filter with cutoff pi*Wc radians:
[b, a] = butter (n, Wc, "high")
Band pass filter with edges pi*Wl and pi*Wh radians:
[b, a] = butter (n, [Wl, Wh])
Band reject filter with edges pi*Wl and pi*Wh radians:
[b, a] = butter (n, [Wl, Wh], "stop")
Return filter as zero-pole-gain rather than coefficients of the numerator and denominator polynomials:
[z, p, g] = butter (...)
Return a Laplace space filter, Wc can be larger than 1:
[...] = butter (..., "s")
Return state-space matrices:
[a, b, c, d] = butter (...)
References:
Proakis & Manolakis (1992). Digital Signal Processing. New York: Macmillan Publishing Company.
n = buttord (wp, ws, rp, rs) ¶n = buttord ([wp1, wp2], [ws1, ws2], rp, rs) ¶n = buttord ([wp1, wp2], [ws1, ws2], rp, rs, "s") ¶[n, wc_p] = buttord (…) ¶[n, wc_p, wc_s] = buttord (…) ¶Compute the minimum filter order of a Butterworth filter with the desired response characteristics. The filter frequency band edges are specified by the passband frequency wp and stopband frequency ws. Frequencies are normalized to the Nyquist frequency in the range [0,1]. rp is the allowable passband ripple measured in decibels, and rs is the minimum attenuation in the stop band, also in decibels.
The output arguments n and wc_p (or n and wc_n) can
be given as inputs to butter.
Using wc_p makes the filter characteristic touch at least one pass band
corner and using wc_s makes the characteristic touch at least one
stop band corner.
If wp and ws are scalars, then wp is the passband cutoff frequency and ws is the stopband edge frequency. If ws is greater than wp, the filter is a low-pass filter. If wp is greater than ws, the filter is a high-pass filter.
If wp and ws are vectors of length 2, then wp defines the passband interval and ws defines the stopband interval. If wp is contained within ws (ws1 < wp1 < wp2 < ws2), the filter is a band-pass filter. If ws is contained within wp (wp1 < ws1 < ws2 < wp2), the filter is a band-stop or band-reject filter.
If the optional argument "s" is given, the minimum order for an analog
elliptic filter is computed. All frequencies wp and ws are
specified in radians per second.
Theory: For Low pass filters, |H(W)|^2 = 1/[1+(W/Wc)^(2N)] = 10^(-R/10). With some algebra, you can solve simultaneously for Wc and N given Ws,Rs and Wp,Rp. Rounding N to the next greater integer, one can recalculate the allowable range for Wc (filter characteristic touching the pass band edge or the stop band edge).
For other types of filter, before making the above calculation, the requirements must be transformed to LP requirements. After calculation, Wc must be transformed back to original filter type.
See also: butter, cheb1ord, cheb2ord, ellipord.
(n, x) ¶Returns the value of the nth-order Chebyshev polynomial calculated at the point x. The Chebyshev polynomials are defined by the equations:
/ cos(n acos(x), |x| <= 1
Tn(x) = |
\ cosh(n acosh(x), |x| > 1
If x is a vector, the output is a vector of the same size, where each element is calculated as y(i) = Tn(x(i)).
[z, p, g] = cheb1ap (n, Rp) ¶Design lowpass analog Chebyshev type I filter.
This function exists for MATLAB compatibility only, and is equivalent
to cheby1 (n, Rp, 1, "s").
Input:
Output:
Example
[z, p, g] = cheb1ap (2, 1) z = [](0x1) p = -0.54887 - 0.89513i -0.54887 + 0.89513i g = 0.98261
See also: buttap, cheby1, cheb2ap, ellipap.
n = cheb1ord (wp, ws, rp, rs) ¶n = cheb1ord ([wp1, wp2], [ws1, ws2], rp, rs) ¶n = cheb1ord ([wp1, wp2], [ws1, ws2], rp, rs, "s") ¶[n, wc] = cheb1ord (…) ¶[n, wc_p, wc_s] = cheb1ord (…) ¶Compute the minimum filter order of a Chebyshev type I filter with the desired response characteristics. The filter frequency band edges are specified by the passband frequency wp and stopband frequency ws. Frequencies are normalized to the Nyquist frequency in the range [0,1]. rp is the allowable passband ripple measured in decibels, and rs is the minimum attenuation in the stop band, also in decibels.
The output arguments n and wc_p (or n and wc_s) can
be given as inputs to cheby1.
Using wc_p makes the filter characteristic touch at least one pass band
corner and using wc_s makes the characteristic touch at least one
stop band corner.
If wp and ws are scalars, then wp is the passband cutoff frequency and ws is the stopband edge frequency. If ws is greater than wp, the filter is a low-pass filter. If wp is greater than ws, the filter is a high-pass filter.
If wp and ws are vectors of length 2, then wp defines the passband interval and ws defines the stopband interval. If wp is contained within ws (ws1 < wp1 < wp2 < ws2), the filter is a band-pass filter. If ws is contained within wp (wp1 < ws1 < ws2 < wp2), the filter is a band-stop or band-reject filter.
If the optional argument "s" is given, the minimum order for an analog
elliptic filter is computed. All frequencies wp and ws are
specified in radians per second.
See also: buttord, cheby1, cheb2ord, ellipord.
[z, p, g] = cheb2ap (n, Rs) ¶Design lowpass analog Chebyshev type II filter.
This function exists for MATLAB compatibility only, and is equivalent
to cheby2 (n, Rs, 1, "s").
Demo
demo cheb2ap
See also: cheby2.
n = cheb2ord (wp, ws, rp, rs) ¶n = cheb2ord ([wp1, wp2], [ws1, ws2], rp, rs) ¶n = cheb2ord ([wp1, wp2], [ws1, ws2], rp, rs, "s") ¶[n, wc_s] = cheb2ord (…) ¶[n, wc_s, wc_p] = cheb2ord (…) ¶Compute the minimum filter order of a Chebyshev type II filter with the desired response characteristics. The filter frequency band edges are specified by the passband frequency wp and stopband frequency ws. Frequencies are normalized to the Nyquist frequency in the range [0,1]. rp is the allowable passband ripple measured in decibels, and rs is the minimum attenuation in the stop band, also in decibels.
The output arguments n and wc_p (or n and wc_s) can
be given as inputs to cheby2.
Using wc_p makes the filter characteristic touch at least one pass band
corner and using wc_s makes the characteristic touch at least one
stop band corner.
If wp and ws are scalars, then wp is the passband cutoff frequency and ws is the stopband edge frequency. If ws is greater than wp, the filter is a low-pass filter. If wp is greater than ws, the filter is a high-pass filter.
If wp and ws are vectors of length 2, then wp defines the passband interval and ws defines the stopband interval. If wp is contained within ws (ws1 < wp1 < wp2 < ws2), the filter is a band-pass filter. If ws is contained within wp (wp1 < ws1 < ws2 < wp2), the filter is a band-stop or band-reject filter.
If the optional argument "s" is given, the minimum order for an analog
elliptic filter is computed. All frequencies wp and ws are
specified in radians per second.
See also: buttord, cheb1ord, cheby2, ellipord.
[b, a] = cheby1 (n, rp, w) ¶[b, a] = cheby1 (n, rp, w, "high") ¶[b, a] = cheby1 (n, rp, [wl, wh]) ¶[b, a] = cheby1 (n, rp, [wl, wh], "stop") ¶[z, p, g] = cheby1 (…) ¶[a, b, c, d] = cheby1 (…) ¶[…] = cheby1 (…, "s") ¶Generate a Chebyshev type I filter with rp dB of passband ripple.
[b, a] = cheby1(n, Rp, Wc) low pass filter with cutoff pi*Wc radians
[b, a] = cheby1(n, Rp, Wc, ’high’) high pass filter with cutoff pi*Wc radians
[b, a] = cheby1(n, Rp, [Wl, Wh]) band pass filter with edges pi*Wl and pi*Wh radians
[b, a] = cheby1(n, Rp, [Wl, Wh], ’stop’) band reject filter with edges pi*Wl and pi*Wh radians
[z, p, g] = cheby1(...) return filter as zero-pole-gain rather than coefficients of the numerator and denominator polynomials.
[...] = cheby1(...,’s’) return a Laplace space filter, W can be larger than 1.
[a,b,c,d] = cheby1(...) return state-space matrices
References:
Parks & Burrus (1987). Digital Filter Design. New York: John Wiley & Sons, Inc.