This is FilterSupport.m in view mode; [Download] [Up]
(* :Title: One Dimensional Filters *)
(* :Authors: Wally McClure, Brian Evans, James McClellan *)
(* :Summary: Supporting code for designing one-dimensional filters *)
(* :Context: SignalProcessing`Support`FilterSupport` *)
(* :PackageVersion: 2.6 *)
(*
:Copyright: Copyright 1989-1991 by Brian L. Evans
Georgia Tech Research Corporation
Permission to use, copy, modify, and distribute this software
and its documentation for any purpose and without fee is
hereby granted, provided that the above copyright notice
appear in all copies and that both that copyright notice and
this permission notice appear in supporting documentation,
and that the name of the Georgia Tech Research Corporation,
Georgia Tech, or Georgia Institute of Technology not be used
in advertising or publicity pertaining to distribution of the
software without specific, written prior permission. Georgia
Tech makes no representations about the suitability of this
software for any purpose. It is provided "as is" without
express or implied warranty.
*)
(*
:History: date begun -- December 24, 1989
cleaned up for premier release -- October 29, 1990
added functional definitions for FIR/IIR -- March, 1991
*)
(* :Keywords: filter, design *)
(*
:Source: {Theory and Application of Digital Signal Processing}
by Lawrence Rabiner and Bernard Gold, 1975
*)
(* :Warning: *)
(* :Mathematica Version: 1.2 or 2.0 *)
(*
:Functions: BilinearTransformation
FilterOrder
FilterParameters
IIRFunction
*)
If [ TrueQ[ $VersionNumber >= 2.0 ],
Off[ General::spell ];
Off[ General::spell1 ] ]
(* B E G I N P A C K A G E *)
BeginPackage[ "SignalProcessing`Support`FilterSupport`",
"SignalProcessing`Support`SupCode`",
"SignalProcessing`Support`SigProc`" ]
(* U S A G E I N F O R M A T I O N *)
Analog::usage =
"Analog is continuous-time signal processing."
Bandpass::usage =
"Bandpass describes a filter that passes a band of frequencies \
and rejects all others."
Bandstop::usage =
"Bandstop describes a filter that rejects a band of frequencies \
and passes all others."
Bessel::usage =
"Bessel is a type of infinite impulse response (IIR) filter. \
The step response of a lowpass Bessel function has low overshoot. \
The transfer function is of the form d0 / Bn(s), where d0 is a \
constant and Bn(s) is the nth-order Bessel polynomial. \
The magnitude response of a lowpass Butterworth filter is smooth \
(no ripples) and monotonically decreases with respect to frequency."
BilinearTransformation::usage =
"BilinearTransformation[omega, w, k] assumes that two of three \
arguments have value and will set the value of the unknown parameter. \
It returns a list of the three arguments. \
BilinearTransformation[omega, w, k, h, s, z, t] will transform \
the analog transfer function h = f(s) into a digital transfer \
function g(z). \
Now, if h is instead a function of t and not of s, \
then the Bilinear transformation will be applied to the Laplace \
transform of h. If h is an analog filter (see CIIR and CFIR), \
then the result will be a digital filter (see IIR and FIR)."
Butterworth::usage =
"Butterworth is a type of infinite impulse response (IIR) filter \
whose magnitude squared response is of the form 1/(1 + (w/wc)^(2 n)) \
where w is the frequency variable, wc is the cutoff frequency, \
and n is the filter order. \
The magnitude response of a lowpass Butterworth filter is smooth \
(no ripples) and monotonically decreases with respect to frequency."
CFIR::usage =
"CFIR[t, {h0, h1, h2, ...}] represents the finite impulse response \
of an all-zero analog filter. \
Input to the filter is specified as FIR[t, h] [ x[t] ], \
where x[t] is some symbolic signal processing expression. \
For the generalized N-dimensional FIR structure, \
the argument t becomes a list of indices (like {t1, t2} for N=2) \
and the feedback coefficients are in an N x N matrix. \
Supports option Roots, which for 1-D filters is a list of zeroes."
CIIR::usage =
"CIIR[t, {a0, a1, a2, ...}] represents the infinite impulse response \
of an all-pole analog filter with gain 1, where a1, a2, ..., are the \
feedback coefficients and a0 is the scaling factor of the output. \
The Laplace transform of such an expression is \
1 / ( a0 + a1 s + a2 s^2 + ... ). \
By default, the input to the filter is a Dirac delta function. \
An alternate input h can be specified by the code \
CIIR[t, {a0, a1, a2, ...}][h]. \
For the generalized N-dimensional IIR structure, the argument t \
becomes a list of indices (like {t1, t2} for N=2) and the feedback \
coefficients are specified by an N x N matrix. \
Supports option Roots, which for 1-D filters is a list of poles."
ChebyshevI::usage =
"ChebyshevI is a type of infinite impulse response (IIR) filter \
whose magnitude squared response is of the form \
1/(1 + e^2 Tn(w/wc)^2) where e controls passband ripple, \
w is the frequency variable, and wc is the cutoff frequency. \
The term Tn(s) is the nth-order Chebyshev polynomial. \
The magnitude response of the stopband monotonically decreases \
with respect to frequency."
ChebyshevII::usage =
"ChebyshevII is a type of infinite impulse response (IIR) filter \
whose magnitude squared response is of the form \
1/(1 + e^2 Xn(w/wc)^2) where e controls stopband ripple, \
w is the frequency variable, and wc is the cutoff frequency. \
Here, Xn(w) = Tn(wr) / Tn(wr / w), where Tn(w) is the nth-order \
Chebyshev polynomial. \
The magnitude response of the passband \
monotonically decreases with respect to frequency."
DesignFilter::usage =
"DesignFilter[domain, class, type, time-variable, deltap, \
deltas, w1, w2, w3, w4] returns a filter expression (involving \
IIR/FIR or CIIR/CFIR filter sections). \
Domain is either Analog or Digital. \
Possibilities for class include Butterworth, \
ChebyshevI, ChebyshevII, and Elliptic. \
Possible filter types are Lowpass, Highpass, Bandpass, and Bandstop. \
For lowpass filters, the magnitude response of the passband resides \
on [1, 1 - deltap] and the magnitude response of the stopband \
occupies [0, deltas]. \
The last two arguments, w3 and w4, are only necessary for \
bandpass and bandstop filters. \
Frequencies w1, w2, w3, and w4 are measured in rad/sec and \
should be given in ascending order."
Digital::usage =
"Digital is discrete-time signal processing."
Elliptic::usage =
"Elliptic is a type of infinite impulse response (IIR) filter \
whose magnitude squared response is of the form \
1/(1 + e^2 Rn(w,L)^2) where e controls ripple, w is the frequency \
variable, and wc is the cutoff frequency. \
Here, Rn(w, L) is called a Chebyshev rational function \
such that L is a parameter describing ripple properties. \
Elliptic filters have equiripple in the passbands and stopbands."
FilterOrder::usage =
"FilterOrder[filter-type, transition-ratio, nu] returns the filter \
order for the lowpass filter that meets the transition ratio \
(for digital filters, see BilinearTransformation) \
and the filter parameter nu (which is a function of the passband \
and stopband ripple parameters-- c.f. FilterParameters)."
FilterParameters::usage =
"FilterParameters[delta1, delta2, epsilon, A, Eripple, ATT, nu] \
fills in the values of the missing parameters. \
It assumes that either the passband ripple (delta1) and \
the stopband ripple (delta2) or epsilon (controls ripple) and \
A (1 / delta2) are given. \
The values of Eripple, ATT, and nu will be computed and assigned \
to the symbols passed into those slots. \
Note that the passband takes on values between 1 and (1 - delta1) \
and that the stopband takes on values between 0 and delta2. \
The function returns a list of all seven parameters."
FIR::usage =
"FIR[n, {h0, h1, h2, ...}] represents the finite impulse response \
of an all-zero digital filter. \
As a function of n, the response is actually equal to \
h0, h1, h2, ..., for n = 0, 1, 2, ....
Note that FIR[n, 2] will be simplified to 2 Impulse[n]. \
FIR can also represent an all-zero filter. \
In this case, the input x[n] to the filter is specified as
FIR[n, h] [ x[n] ], where x[n] is some symbolic \
signal processing expression. \
For the generalized N-dimensional FIR structure, \
the argument n becomes a list of indices (like {n1, n2} for N=2)
and the feedback coefficients are in an N x N matrix."
Highpass::usage =
"Highpass describes a filter that passes all frequencies above \
a cutoff frequency and rejects all others."
IIR::usage =
"IIR[n, {a0, a1, a2, ...}] represents the infinite impulse response \
of an all-pole filter with gain 1, where a1, a2, ..., are the \
feedback coefficients and a0 is the scaling factor of the output. \
In discrete time, the output of IIR is the solution to y[n] \
in the difference equation a0 y[n] + a1 y[n-1] + ... + aN y[n - N] = \
impulse[n], with y[n] = 0 for n < 0. \
By default, the input into the filter (a.k.a. the forcing function) \
is an impulse. \
An alternate input h can be specified by the code \
IIR[n, {a0, a1, a2, ...}][h]. \
For the generalized N-dimensional IIR structure, the argument n \
becomes a list of indices (like {n1, n2} for N=2) and the feedback \
coefficients are specified by an N x N matrix."
IIRFunction::usage =
"IIRFunction[n, coeffs] and IIRFunction[n, coeffs, input] generate \
the object IIRFunction[n, coeffs, input, TheFunction -> hiir] that \
represents a digital IIR filter with feedback coefficients coeffs \
(see IIR for meaning) and input function input. \
When n is an integer or real-valued integer, IIRFunction will return \
the value of the filter computed at that sample. \
For negative n, the value of IIRFunction is zero but the \
value of hiir is the value of the input function (this allows \
initial conditions to be specified in the input function). \
Here, hiir is the actual recursive object that computes \
values of the IIR filter only as needed. \
To see what values have already been computed, type ??<hiir_slot>."
Lowpass::usage =
"Lowpass describes a filter that passes all frequencies below \
a cutoff frequency and rejects all others."
(* E N D U S A G E I N F O R M A T I O N *)
CFIR::invalid =
"Mismatch between the number of roots and the number of coefficients."
CIIR::invalid =
"Mismatch between the number of roots and the number of coefficients."
Begin[ "`Private`" ]
(* Augment convolution operator *)
Unprotect[Convolve]
Convolve[n_][FIR[n_, coeffs_, rest___], x_] :=
FIR[n, coeffs, rest][x]
Convolve[n_][x_, FIR[n_, coeffs_, rest___]] :=
FIR[n, coeffs, rest][x]
Convolve[n_][FIR[n_, coeffs_, rest___], x_, right__] :=
FIR[n, coeffs, rest][Convolve[n][x, right]]
Convolve[n_][IIR[n_, coeffs_, rest___], x_] :=
IIR[n, coeffs, rest][x]
Convolve[n_][x_, IIR[n_, coeffs_, rest___]] :=
IIR[n, coeffs, rest][x]
Convolve[n_][IIR[n_, coeffs_, rest___], x_, right__] :=
IIR[n, coeffs, rest][Convolve[n][x, right]]
Protect[Convolve]
(* BilinearTransformation *)
SetAttributes[BilinearTransformation, HoldAll]
BilinearTransformation[omega_, w_, k_] :=
Block [ {},
Which [ ! ValueQ[omega], omega = k Tan[w/2],
! ValueQ[w], w = 2 ArcTan[omega / k],
! ValueQ[k], k = omega / Tan[w/2] ];
{ omega, w, k } ]
BilinearTransformation[omega_, w_, k_, h_, s_:Global`s, z_:Global`z] :=
Block [ {expr},
BilinearTransformation[omega, w, k];
expr = k ( 1 - z^-1 ) / ( 1 + z^-1 );
h /. s -> expr ]
(* FilterOrder -- see [Rabiner & Gold, 241] *)
FilterOrder[ Butterworth, k_, k1_ ] := Max[1, Ceiling[ N[ Log[k1] / Log[k]]]]
FilterOrder[ ChebyshevI, k_, k1_ ] :=
Max[1, Ceiling[N[ArcCosh[1 / k1] / Log[ (1 + Sqrt[1 - k^2]) / k ]]]]
FilterOrder[ ChebyshevII, k_, k1_ ] :=
Max[1, Ceiling[N[ArcCosh[1 / k1] / Log[ (1 + Sqrt[1 - k^2]) / k ]]]]
FilterOrder[ Elliptic, k_, k1_ ] :=
Max[ 1, Ceiling[ N[ EllipticK[k^2] EllipticK[1 - k1^2] /
( EllipticK[k1^2] EllipticK[1 - k^2] ) ] ] ]
(* FilterParameters -- see [Rabiner & Gold, 239-240] *)
SetAttributes[FilterParameters, HoldAll]
FilterParameters[ delta1_, delta2_, epsilon_, A_, Eripple_, ATT_, nu_ ] :=
Block [ {},
(* Parameters delta1 and epsilon are related *)
If [ NumberQ[N[delta1]] && SameQ[Head[epsilon], Symbol],
epsilon = N[Sqrt[2 - delta1] Sqrt[delta1] / (1 - delta1)]];
If [ NumberQ[N[epsilon]] && SameQ[Head[delta1], Symbol],
delta1 = N[1 - 1 / Sqrt[1 + epsilon^2]]];
(* Parameters delta2 and A are related *)
If [ NumberQ[N[delta2]] && SameQ[Head[A], Symbol],
A = N[1 / delta2]];
If [ NumberQ[N[A]] && SameQ[Head[delta2], Symbol],
delta2 = N[1 / A]];
(* Now, determine values of Eripple, ATT, nu *)
Eripple = N[ 20 Log[10, Sqrt[1 + epsilon]] ];
ATT = N[ 20 Log[10, A] ];
nu = N[ epsilon / Sqrt[A^2 - 1] ];
(* Return a list of all seven parameters *)
{ delta1, delta2, epsilon, A, Eripple, ATT, nu } ]
(* degenfilterq -- degenerate filter check *)
degenfilterq[var_, coeffs_] :=
(! ListQ[coeffs]) && (! PatternQ[var]) &&
(! PatternQ[coeffs]) && FreeQ[coeffs, var]
(* CFIR *)
CFIR/: CFIR[t_Symbol, h0_] := h0 Delta[t] /; degenfilterq[t, h0]
CFIR/: Format[ CFIR[t_, h_, Roots -> r_] ] := CFIR[t, h]
(* CIIR *)
CIIR/: CIIR[t_Symbol, a0_] := 1 / a0 /; degenfilterq[t, a0]
CIIR/: Format[ CIIR[t_, a_, Roots -> r_] ] := CIIR[t, a]
(* FIR *)
FIR/: FIR[n_Symbol, h0_] := h0 Impulse[n] /; degenfilterq[n, h0]
FIR/: TheFunction[ FIR[n_Symbol, h_][x_] ] :=
Block [ {i, input, length},
length = Length[h];
input = TheFunction[x] /. ( n -> n - i + 1 );
Dot[ h, Table[input, {i, 1, length}] ] ]
FIR/: TheFunction[ FIR[n_Symbol, h_] ] := SequenceToFunction[h, n] /;
FreeQ[h, n]
(* IIR *)
IIR/: IIR[n_Symbol, a0_] := 1 / a0 /; degenfilterq[n, a0]
IIR/: TheFunction[ IIR[n_Symbol, a_][x_] ] := IIRFunction[n, a, x]
IIR/: TheFunction[ IIR[n_Symbol, a_] ] := IIRFunction[n, a]
(* IIRFunction *)
(* two issues here -- initial values on filter *)
(* the computational function for non-negative n *)
(* IIRFunction computes values for all n, hiir for non-negative n *)
(* We represent hiir as an option for readability and so that *)
(* the object work correctly with GetVariables. *)
IIRFunction[n_Symbol, a_List] := IIRFunction[n, a, Impulse[n]]
IIRFunction[n_Symbol, a_List, in_] :=
Block [ {hiir, i, input, length, str, termstr},
input = TheFunction[in];
hiir = Unique["hiir"];
str = ToString[ StringForm[ "``[n_Integer] := `` /; n < 0",
hiir, InputForm[input] ] ];
GenerateCode[ str ];
str = ToString[ StringForm[ "``[n_Integer] := ``[n] = (``",
hiir, hiir, InputForm[input] ] ];
length = Length[a];
For [ i = 2, i <= length, i++,
termstr = ToString[ StringForm[ "+ (``) * ``[n - ``]",
InputForm[a[[i]]], hiir, i-1 ] ];
str = StringJoin[str, termstr] ];
str = StringJoin[ str,
ToString[StringForm[ ") / (``)",
InputForm[ a[[1]] ] ]] ];
GenerateCode[ str ];
IIRFunction[n, a, input, TheFunction -> hiir] ]
IIRFunction[n_, a_List, input_, op_] := 0 /; n < 0
IIRFunction[n_Integer, a_List, input_, TheFunction -> hiir_] := hiir[n]
IIRFunction[n_Integer, a_List, input_, TheFunction :> hiir_] := hiir[n]
IIRFunction[n_Real, a_List, input_, TheFunction -> hiir_] :=
hiir[Rationalize[n]] /;
IntegerQ[Rationalize[n]]
IIRFunction[n_Real, a_List, input_, TheFunction :> hiir_] :=
hiir[Rationalize[n]] /;
IntegerQ[Rationalize[n]]
(* E N D P A C K A G E *)
End[]
EndPackage[]
If [ TrueQ[ $VersionNumber >= 2.0 ],
On[ General::spell ];
On[ General::spell1 ] ]
(* H E L P I N F O R M A T I O N *)
Block [ {newfuns},
newfuns = {BilinearTransformation, FilterOrder,
FilterParameters, IIRFunction};
Combine[ SPfunctions, newfuns ];
Apply[ Protect, newfuns ] ]
Block [ {newoperators},
newoperators = {FIR, IIR};
Combine[ SPoperators, newoperators ];
Apply[ Protect, newoperators ] ]
Block [ {newsignals},
newsignals = {FIR, IIR};
Combine[ SPsignals, newsignals ];
Apply[ Protect, newsignals ] ]
(* E N D I N G M E S S A G E *)
Print[ "Supporting routines for filter design are loaded." ]
Null
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.