Signal Processing with Mathematica

Mariusz Jankowski

University of Southern Maine and Wolfram Research

What is signal processing?

The analysis, design, and verification of frequency-selective filters.

A fifth-order filter composed of two second-order Sallen-Key stages and a first-order input stage (left) and model of second-order digital filter (right).

Row[{,}]

Time-domain and frequency domain analysis of analog filters.

h=Signal`EquirippleFilterArray[{{0.,0.13,0.18,0.33,0.37,0.5} 2 π,{0,1,0},{1,1,10}},51];
{Manipulate[Plot[Evaluate[InverseLaplaceTransform[Signal`ButterworthFilterModel[n,s][s][[1,1]],s,t]//Chop],{t,0,30},PlotRange->{{0,30},All},AspectRatio->1/3],{{n,4},2,13,1}],Plot[Abs[Signal`ListDTFT[h,ω]//N],{ω,0,π},PlotRange->All,GridLines->Automatic],Plot[20 Log[10,Abs[Signal`ListDTFT[h,ω]//N]],{ω,0,π},PlotRange->All,GridLines->Automatic]}

Simulation of an analog filter: example circuits created in PSpice (left) and MathModelica (right).

Spectral analysis/visualization.

Signal processing foundations

Mathematica already has many of the essential foundational functionality to support advanced symbolic and numerical signal processing.

Time-domain modeling and analysis

Filters are represented by difference and differential equations.

Here is the step response of a (2^nd)-order Butterworth filter.

DSolve[{y''[t]+(Sqrt[2])y'[t]+y[t]==UnitStep[t],y'[0]==0,y[0]==0},y[t],t]

Plot[(E^(-(t/(Sqrt[2])))) ((E^(t/(Sqrt[2])))-Cos[(t/(Sqrt[2]))]-Sin[(t/(Sqrt[2]))]) UnitStep[t],{t,0,15},PlotRange->All,AspectRatio->1/3]

Impulse response of a discrete-time model of a (2^nd)-order Butterworth filter.

h=RecurrenceTable[{7.83 y[n]-6. y[n-1]+2.17 y[n-2]==DiscreteDelta[n]+2. DiscreteDelta[n-1]+DiscreteDelta[n-2],y[-1]==0,y[-2]==0},y[n],{n,0,16}]

Convolution. This shows the response of our (2^nd)-order Butterworth filter to data taken from a row of an example image. This filter has the effect of smoothing the data.

tmp=ListConvolve[h,row=ImageData[][[120]],1,0];

ListLinePlot[{row,%},PlotRange->All,AspectRatio->1/3]

Frequency-domain modeling and analysis

Classical frequency domain representations of filters are based on a variety of integrals and sums: Laplace, z, and Fourier (both continuous- and discrete-time).

Here is the system transfer function of a Butterworth filter.

B[s_]:=(1/(1+(Sqrt[2])s+(s^2)))

This gives the impulse response of the filter.

InverseLaplaceTransform[B[s],s,t]

This computes and plots the continuous-time Fourier transform of the impulse response, namely, the frequency response of the filter.

FourierTransform[% UnitStep[t],t,ω,FourierParameters->{1,-1}]

GraphicsRow[{Plot[((Abs[((Sqrt[2])/((Sqrt[2])+2 I ω-(Sqrt[2]) (ω^2)))])^2),{ω,0,4},PlotRange->All,AxesLabel->{ω,None},ImageSize->Small],Plot[20Log[10,Abs[((Sqrt[2])/((Sqrt[2])+2 I ω-(Sqrt[2]) (ω^2)))]],{ω,0,4},PlotRange->All,GridLines->Automatic,AxesLabel->{ω,None},ImageSize->Small]}]

Here are the pole locations of an n-th order Butterworth filter - it is well-known that they lie on a circle.

Manipulate[ListPlot[{Re[#],Im[#]}&/@First@First@TransferFunctionPoles[Signal`ButterworthFilterModel[n,s]],PlotRange->{{-1.5,1.5},{-1.5,1.5}},PlotStyle->{Red,PointSize[Large]},AspectRatio->1,AxesLabel->{"Re","Im"},Epilog->{Pink,Circle[{0,0},1]},ImageSize->Small],{{n,11},2,31,1}]

System transfer function of a discrete-time filter.

G[z_]:=((((1+z))^2)/(2-6z+8 (z^2)))

Impulse response of the filter.

InverseZTransform[G[z],z,n]UnitStep[n]//Simplify

Frequency response of the filter.

FourierSequenceTransform[%,n,ω]

Plot[((Abs[%])^2),{ω,-π,2 π},PlotRange->All,AspectRatio->1/3,AxesLabel->{ω,None}]

Other related functionality - wavelets, image processing, statistics, and control systems

Signal processing is leveraging some of the functionality already available in control systems

TransferFunctionModel

OutputResponse

Conversion from analog to digital filters using ToDiscreteTime

Visualization with BodePlot, SingularValuePlot, PoleZeroPlot, more ...

Project goals

To build the necessary specialized functionality to support design and analysis of analog and digital filters.

Analog filter design

ButterworthFilterModel, ChebyshevIFilterModel, ChebyshevIIFilterModel, EllipticFilterModel

Digital filter design

FrequencySamplingFilterArray, LeastSquaresFilterArray, EquirippleFilterArray

Filtering

RecursiveFilter

Transforms

ListDTFT, ListZTransform, more

Visualization

SpectrogramArrayPlot, SpectrogramArray

Filter design - analog filters

Filters are vital in modern electronic, they are used in a wide variety of applications. In the field of telecommunication, band-pass filters are used in the audio frequency range (0 kHz to 20 kHz) for modems and speech processing. High-frequency band-pass filters (several hundred MHz) are used for channel selection in telephone central offices. Data acquisition systems usually require anti-aliasing low-pass filters as well as low-pass noise filters in their preceding signal conditioning stages. System power supplies often use band-rejection filters to suppress the 60-Hz line frequency and high frequency transients.

Four filter types are important: Butterworth, Chebyshev I, Chebyshev II, and elliptic. They all share a common squared magnitude response fromulation:

| H(j ω) (|^2)= (1/(1 + ((Subscript[A, n])^2)(ω)))

where ((Subscript[A, n])^2)(ω) is a polynomial of order 2n which is unique to each of the analog filters and takes the following forms,

(Subscript[A, n])(ω) = {({{(ω^( n)), Butterworth,}, {ϵ (Subscript[T, n])(ω), Chebyshev I,}, {(1/(ϵ (Subscript[T, n])((1/ω)))), Chebyshev  II,}, {ϵ (Subscript[R, n])(ω), elliptic,}})

where (Subscript[T, n])(ω) is the Chebyshev polynomial of the first kind of order n and (Subscript[R, n])(ω) is the Chebyshev rational function.

This shows four realizations of a lowpass filter with cutoff frequency at 1 rad/sec.

Plot[Evaluate[Abs[#[{"Lowpass",{1.,1.3},{1.,20.}},s][I ω]//First//First]],{ω,0,4},PlotRange->{0,1},PlotLabel->#,AspectRatio->1/3,AxesLabel->{ω,None},Epilog->{Opacity[0.3],Pink,Rectangle[{0,0},{1.,0.89}],Rectangle[{1.3,0.1},{4,1}]}]&/@{Signal`ButterworthFilterModel,Signal`ChebyshevIFilterModel,Signal`ChebyshevIIFilterModel,Signal`EllipticFilterModel}

First@BodePlot[#[{"Lowpass",{1.,3.},{3.,90.}},s],{0.1,100},PlotRange->{{Full,{-150,5}},Automatic},PlotLabel->#,PlotLayout->List,GridLines->Automatic]&/@{Signal`ButterworthFilterModel,Signal`ChebyshevIFilterModel,Signal`ChebyshevIIFilterModel,Signal`EllipticFilterModel}

Function FrequencyTransformation can be used to convert any lowpass filter prototype into corresponding highpass, bandpass or bandstop filters. Here are four filter types all related to the same lowpass prototype.

tmp=({#,Signal`FrequencyTransformation[#,s->1/s],Signal`FrequencyTransformation[#,s->(1/10)(((s^2)+1)/s)],Signal`FrequencyTransformation[#,s->((10 s)/((s^2)+1))]}&)[Signal`ChebyshevIFilterModel[{"Lowpass",{1.,3.},{3.,60.}},s]];
MapThread[First@BodePlot[#1,{0.01,100},PlotRange->{{Full,{-150,5}},Automatic},PlotLabel->#2,PlotLayout->List,GridLines->Automatic]&,{tmp,{"Lowpass","Highpass","Bandpass","Bandstop"}}]

Filter design - analog filters

Now that we can obtain the transfer functions of any of the classic analog filters we can proceed with the design and/or verification of real active or passive circuits.

Unity gain Sallen-Key circuit configuration:

Standard nodal analysis gives the following transfer function for this circuit (assuming ideal operational amplifier):

((Subscript[v, out])/(Subscript[v, in]))->(1/(1+s (Subscript[C, 1]) (Subscript[R, 1])+s (Subscript[C, 1]) (Subscript[R, 2])+(s^2) (Subscript[C, 1]) (Subscript[C, 2]) (Subscript[R, 1]) (Subscript[R, 2]))) = (1/(1+(Sqrt[2])s+(s^2)))

Element values (assume (Subscript[R, 1])=(Subscript[R, 2])=1Ω)

Solve[{(Subscript[C, 1])((Subscript[R, 1])+(Subscript[R, 2]))==(Sqrt[2]),(Subscript[C, 1]) (Subscript[C, 2]) (Subscript[R, 1]) (Subscript[R, 2])==1,(Subscript[R, 1])==1,(Subscript[R, 2])==1},{(Subscript[C, 1]),(Subscript[C, 2]),(Subscript[R, 1]),(Subscript[R, 2])}]

Mathematica simulation. Filter response to sinusoidal excitation over a range of frequencies. (See also demonstration)

Manipulate[Plot[Evaluate[{Sin[ω t] UnitStep[t],y=First@OutputResponse[(TransferFunctionModel[(1/(1+(Sqrt[2]) s+(s^2))), ]),Sin[ω t] UnitStep[t],{t,0,40/ ω}]}],{t,0,40/ ω},ImageSize->240],{{ω,1.},0.1,5.,Appearance->"Labeled"}]

MathModelica simulation

Using MathModelica the following Sallen-Key Circuit can be designed by drag&drop

If we load MathModelica in Mathematica we can run a simulation.

Needs["MathModelicaLink`"];

sk1=ModelicaSimulate["Amplifiers.SallenKey1", {0, 100},ParameterValues->{"C1"->(1/(Sqrt[2])),"C2"->(Sqrt[2]),"f1"->{1.,3.}/(2 π)}];

{ModelicaPlot[sk1[[1]], {"sineVoltage1.v","idealOpAmp3Pin1.out.v"},{0,30}],ModelicaPlot[sk1[[2]], {"sineVoltage1.v","idealOpAmp3Pin1.out.v"},{0,30}]}

Filter design - analog-to-digital conversion

A common method of creating digital filters is to transform analog prototypes to their digital equivalents using one of several mapping methods.The most popular is the bilinear transformation.

Consider an analog lowpass filter of Chebyshev type with a cutoff frequency at (Subscript[ω, c])= 1 rad/sec.

tfa=Signal`ChebyshevIFilterModel[{"Lowpass",{1.,3.},{1.,40.}},s]

{Plot[Abs[tfa[I ω]],{ω,0,6},AxesLabel->{ω,None},GridLines->Automatic,AspectRatio->1/3],BodePlot[tfa,GridLines->Automatic,PlotLayout->List]//First}

This transforms the analog to filter to a digital one using bilinear transformation. The resulting frequency response, on the normalized domain 0≤ω≤2π, is shown below.

tfd=ToDiscreteTimeModel[tfa,1,z,Method->{"BilinearTransform"}]//Chop

{Plot[Abs[tfd[(E^(I ω))]],{ω,0,2π},AxesLabel->{ω,None},GridLines->Automatic,AspectRatio->1/3],BodePlot[tfd,{0.1,N[2 π]},GridLines->Automatic,PlotLayout->List]//First}

Filter design - least-squares method (aka window method)

The method of approximating a desired frequency response using the Fourier series and applying smoothing windows to diminish the inevitable ripple in the least squares approximation is the oldest direct digital filter design technique available.

Consider a lowpass filter with a cutoff frequency at (Subscript[ω, c])=(π/3). The impulse response is given by the inverse discrete-time Fourier transform.

(1/(2 π))(∫_(-(π/3))^(π/3))(E^(-I ω n))\[DifferentialD]ω

DiscretePlot[If[n==0,(1/3),((Sin[((n π)/3)])/(n π))],{n,-20,20},AxesLabel->{n,Style["h[n]",Italic]},PlotRange->All,AspectRatio->1/3]

The least-squares approximation is of finite length - here a FIR filter of length 21 has been selected.

Signal`LeastSquaresFilterArray[{"Lowpass",(π/3)},21]

Frequency response (red) of least squares filter over a range filter lengths.

Manipulate[tmp=Abs[Signal`ListDTFT[Signal`LeastSquaresFilterArray[{"Lowpass",(π/3.)},n],ω]];
If[type=="absolute",Plot[{UnitBox[((3 ω)/(2 π))],tmp},{ω,0,π},PlotRange->All,Exclusions->None,AxesLabel->{ω,None}],Plot[{Null,20. Log[10.,tmp]},{ω,0,π},PlotRange->{-60,0},GridLines->Automatic]],{{type,"absolute","Magnitude scaling"},{"absolute","logarithmic"}},{{n,21},5,81,1,Appearance->"Labeled"},ContentSize->Automatic]

Filter design - least-squares method (aka window method)

The ripple in the frequency response limits the stopband attenuation to a frequently unacceptable value of -20 dB. To diminish the magnitude of the ripple and therefore improve attenuation, smoothing windows can be applied to the filter.

A large selection of windows is available.

windows = {Signal`BartlettWindow,Signal`BartlettHannWindow,Signal`BlackmanWindow,Signal`BlackmanHarrisWindow,Signal`BlackmanNuttallWindow,Signal`BohmanWindow, Signal`ConnesWindow,Signal`CosineWindow,Signal`DirichletWindow,Signal`ExactBlackmanWindow,Signal`FlatTopWindow,Signal`GaussianWindow,Signal`HammingWindow,Signal`HannWindow,Signal`KaiserWindow,Signal`LanczosWindow,Signal`NuttallWindow, Signal`TukeyWindow, Signal`WelchWindow};
ListAnimate[Plot[#[x],{x,-1,1},PlotLabel->#,Filling->Axis,PlotRange->Full,ImageSize->Small]&/@windows,AnimationRunning->False,AnimationRate->2,AnimationRepetitions->1]

Sampled versions of window functions are needed for filter design. This gives the sampled version of the popular Blackman window function.

Signal`WindowArray[Signal`BlackmanWindow,21]

Show[{Plot[Signal`BlackmanWindow[x],{x,-1/2,1/2},AspectRatio->1/3],ListPlot[Signal`WindowArray[Signal`BlackmanWindow,21],Filling->Axis,PlotStyle->{Red,PointSize[Medium]},DataRange->{-1/2,1/2},AspectRatio->1/3]}]

Here are the results of applying four different windows to our least-squares example filter. Note the improvement in stopband attenuation.

win={Signal`DirichletWindow,Signal`BlackmanWindow,Signal`HannWindow,Signal`KaiserWindow};
Plot[Evaluate[Tooltip[20Log[10,Abs[Signal`ListDTFT[Signal`LeastSquaresFilterArray[{"Lowpass",(π/3)},21] * Signal`WindowArray[#,21],ω]]],#]&/@win],{ω,0,π},PlotRange->{-100,0},GridLines->Automatic]

Filter design - frequency sampling method

This method allows the user to specify the desired filter as a list of amplitude values uniformly spanning the normalized domain 0≤ω≤π.

Consider a simple list of ones and zeros defining the amplitude values. This returns a lowpass filter of length 11.

Signal`FrequencySamplingFilterArray[{1,1,1,0,0,0},11]

Due to symmetry and sampling options the same list of amplitudes can return different filters.

Module[{a,n},a={1,1,1,0,0,0};n=Length[a];
Grid[{{Null,Text["Odd length"],Text["Even length"]},
{Text["Even symmetry"],ListPlot[Signal`FrequencySamplingFilterArray[a,2n-1,"EvenSymmetry"->True,"Shifted"->False],Filling->Axis,PlotRange->All,Ticks->False,Axes->{True,False}],ListPlot[Signal`FrequencySamplingFilterArray[a,2n,"EvenSymmetry"->True,"Shifted"->True],Filling->Axis,PlotRange->All,Ticks->False,Axes->{True,False}]},
{Text["Odd symmetry"],ListPlot[Signal`FrequencySamplingFilterArray[a,2n+1,"EvenSymmetry"->False,"Shifted"->False],Filling->Axis,PlotRange->All,Ticks->False,Axes->{True,False}],ListPlot[Signal`FrequencySamplingFilterArray[a,2n,"EvenSymmetry"->False,"Shifted"->True],Filling->Axis,PlotRange->All,Ticks->False,Axes->{True,False}]}},Spacings->{3,3},Dividers->Center]]

Each filter has different frequency characteristics.

Module[{a,n},a={1,1,1,0,0,0};n=Length[a];
Grid[{{Null,Text["Odd length"],Text["Even length"]},
{Text["Even symmetry"],Show[{ListPlot[Thread[{Table[i ((2π)/(2Length[a]-1)),{i,0,Length[a]-1}],a}],PlotRange->{{0,π},All},PlotStyle->{Red,PointSize[Medium]},Filling->Axis,AxesLabel->{ω,None}],Plot[Evaluate[Abs[Signal`ListDTFT[Signal`FrequencySamplingFilterArray[a,2n-1,"EvenSymmetry"->True,"Shifted"->False],ω]]],{ω,0,π},PlotRange->All]}],Show[{ListPlot[Thread[{Table[(i-(1/2)) ((2π)/(2Length[a])),{i,1,Length[a]}],a}],PlotRange->{{0,π},All},PlotStyle->{Red,PointSize[Medium]},Filling->Axis,AxesLabel->{ω,None}],Plot[Evaluate[Abs[Signal`ListDTFT[Signal`FrequencySamplingFilterArray[a,2n,"EvenSymmetry"->True,"Shifted"->True],ω]]],{ω,0,π},PlotRange->All]}]},
{Text["Odd symmetry"],Show[{ListPlot[Thread[{Table[i ((2π)/(2Length[a]+1)),{i,0,Length[a]}],Prepend[a,0]}],PlotRange->{{0,π},All},PlotStyle->{Red,PointSize[Medium]},Filling->Axis,AxesLabel->{ω,None}],Plot[Evaluate[Abs[Signal`ListDTFT[Signal`FrequencySamplingFilterArray[a,2n+1,"EvenSymmetry"->False,"Shifted"->False],ω]]],{ω,0,π},PlotRange->All]}],Show[{ListPlot[Thread[{Table[(i-(1/2)) ((2π)/(2Length[a])),{i,1,Length[a]}],a}],PlotRange->{{0,π},All},PlotStyle->{Red,PointSize[Medium]},Filling->Axis,AxesLabel->{ω,None}],Plot[Evaluate[Abs[Signal`ListDTFT[Signal`FrequencySamplingFilterArray[a,2n,"EvenSymmetry"->False,"Shifted"->True],ω]]],{ω,0,π},PlotRange->All]}]}},Spacings->{3,3},Dividers->Center]]

Filter design - equiripple method

The Parks, McClellan, Rabiner algorithm is the most recent and most popular of thee three method of designing linear-phase FIR filters presented here. It minimizes the maximum deviation from the desired ideal frequency response and results in filters with equal ripples in each of the bands of the filter.

Here is a multiband filter of length 55 (from McClellan, Parks, and Rabiner (1979)).

Signal`EquirippleFilterArray[{{0.,0.05,.1,0.15,0.18,0.25,0.3,0.36,0.41,0.5} 2 π,{0.,1.,0.,1.,0.},{10.,1.,3.,1.,20.}},55]

Plot[20 Log[10,Abs[Signal`ListDTFT[%,ω]]],{ω,0,π},PlotRange->All,GridLines->Automatic]

Filtering with recursive filters

The response of a recursive digital filter given its model may be computed using several methods: RSolve, RecurrenceTable, and OutputResponse. However, none of these methods are suitable if the input is a vector or array of samples of a signal. Use function RecursiveFilter if you need to apply a recursive filter to sampled data.

Impulse response of a discrete-time model of a (2^nd)-order Butterworth filter.

Signal`RecursiveFilter[{{7.83,-6.,2.17 },{1.,2.,1.}},ArrayPad[{1},{2,20}]]

ListPlot[%[[3;;]],PlotRange->All,AspectRatio->1/3,Filling->Axis,PlotStyle->PointSize[Medium]]

Image smoothing using a 5th-order Butterworth filter with variable cutoff frequency.

Manipulate[Module[{tf,tfZ},tf=Signal`ButterworthFilterModel[{5,w},s];tfZ=ToDiscreteTimeModel[tf,1.,z]//Chop;{Plot[Abs[tfZ[(E^(I ω))]],{ω,0,π},PlotRange->All,GridLines->Automatic,AxesLabel->{ω,None}],Image[Signal`RecursiveFilter[{Reverse@CoefficientList[Denominator[#],z],Reverse@CoefficientList[Numerator[#],z]}&[tfZ[z][[1,1]]],ImageData[]]]}],{{w,1.,"frequency"},0.1,5.}]

Filtering with non recursive filters

Use function ListConvolve (ImageConvolve) when filtering with finite impulse response (non recursive) filters.

Simple equalizer

Manipulate[Module[{snd,data,rate,w,h},
snd=Import["C:\\Users\\Mariusz\\Documents\\WTC2011\\Data\\d3.wav"];
data=snd[[1,1,1]];
rate=snd[[1,2]];
w=Signal`WindowArray[Signal`HannWindow,71];
h=Signal`LeastSquaresFilterArray[{{(1/8.),(1/4.),(1/2.)}π,{a1,a2,a3,a4} }, 71];
ListPlay[Join[data,res=ListConvolve[w * h,data,1,0]],PlayRange->All,SampleRate->rate]],Row[{Control[{{a1,1.},0.,1.,VerticalSlider,ImageSize->{20,100}}],Spacer[10],Control[{{a2,0.5},0.,1.,VerticalSlider,ImageSize->{20,100}}],Spacer[20],Control[{{a3,0.5},0.,1.,VerticalSlider,ImageSize->{20,100}}],Spacer[60],Control[{{a4,1.},0.,1.,VerticalSlider,ImageSize->{20,100}}]}],ContinuousAction->False]

Here is a comparison of the spectra of the original and filtered signals.

{Signal`SpectrogramArrayPlot[data,SampleRate->rate,ImageSize->300],Signal`SpectrogramArrayPlot[res,SampleRate->rate,ImageSize->300]}

Filtering an image

Manipulate[Module[{i,h,win,fir},i=;
h=Signal`LeastSquaresFilterArray[{type,w }, 21];
win=If[window=="None",Signal`WindowArray[Signal`DirichletWindow,21],Signal`WindowArray[Signal`BlackmanWindow,21]];
fir=If[adj,((#)/(Total[#])),#]&[h win];GraphicsRow[{Plot[{Abs[Signal`ListDTFT[ fir,ω]],If[type=="Lowpass",UnitBox[(ω/(2w))],1-UnitBox[(ω/(2w))]]},{ω,0,π},Exclusions->None,PlotRange->All],ImageConvolve[i,Outer[Times,fir,fir]]}]],{{window,"Blackman","window type"},{"None","Blackman"}},{{type,"Lowpass","filter type"},{"Lowpass","Highpass"}},{{adj,False,"normalize"},{False,True}},{{w,N[(π/3)],"cutoff frequency"},0.1,0.9π,Appearance->"Labeled"}]

Thank you

Acknowledgments: Piotr Wendykier, Shadi Ashnai, and Jan Brugard