Essays, Posts & Presentations

Deploying Bandpass Filters Using the Wolfram Language Microcontroller Kit

Suba Thomas

Author

Suba Thomas

Title

Deploying Bandpass Filters Using the Wolfram Language Microcontroller Kit

Description

How to compute, analyze and deploy real-time filters to Arduinos directly from the Wolfram Language. Complete code provided for the entire workflow.

Deploying Bandpass Filters Using the Wolfram Language Microcontroller Kit

September 24, 2019

Suba Thomas, Software Engineer, Algorithms R&D

Real-time filters work like magic. Usually out of sight, they clean data to make it useful for the larger system they are part of, and sometimes even for human consumption. A fascinating thing about these filters is that they don’t have a big-picture perspective. They work wonders with only a small window into the data that is streaming in. On the other hand, if I had a stream of numbers flying across my screen, I would at the very least need to plot it to make sense of the data. These types of filters are very simple as well.

Take a basic lowpass filter based on two weights, 49/50 and 1/50. The filter computes its next output as a weighted sum of the current output and input. Specifically,

output

k+1

=49/50

output

+1/50

input

. So with a few elementary operations, it essentially strips away the high-frequency component on a properly sampled input signal:

In[]:=

With{u=Table[Sin[t/10]+Sin[10t],{t,0,150,1/10.}]},ListLinePlotu,RecurrenceFilter1,-

,

,u,PlotLegends{"original signal","filtered signal"}

Out[]=

	original signal
	filtered signal

In hindsight, this output makes sense. The filter relies only a bit on the current input, and hence misses the high-frequency portion of the signal, effectively getting rid of it. Any slowly changing component has already been captured in the output.

The filter weights need not be conjured up or arrived at by trial and error. They can be obtained in a systematic fashion using the signal processing functions in the Wolfram Language. (And that’s what I will be doing in the next section for a couple of nontrivial filters.) I can also compute the frequency response of the filter, which will tell me how it will shift and smother or amplify a signal. This creates a ballpark of what to expect from the real-time responses of the filter. I can simulate its real-time response to various signals as well. But how does the filter perform in real time? Can I see the shifting and smothering? And how does this compare to what is predicted by the frequency response?

To investigate all this, I will pick some slightly more involved filters—namely a bandpass Butterworth and a Chebyshev 1 filter—and deploy them to an Arduino Nano. To go about this, I need several things. Luckily, they are all now available in the Wolfram Language.

The functionality to analyze and design filters has existed in the Wolfram Language for quite some time now. The next link in the workflow is the C or Arduino filter code that needs to be deployed to the microcontroller. Thanks to the Microcontroller Kit, which was released with Version 12, I can generate and deploy the code directly from the Wolfram Language. It will automatically generate the needed microcontroller code, sparing me the wearisome task of having to write the code and make sure that I have the coefficients correct. And finally, I use the Wolfram Device Framework to transmit and receive the data from the Nano.

In the first part of this blog post, I will compute, analyze and deploy the filters. In the second part, I will acquire and analyze the filtered data to visualize the responses and evaluate the performance of the filters.

Compute, Analyze and Deploy the Filter I start off by creating a function to compute a bandpass filter that will pass signals with frequencies between 2 Hz and 5 Hz and attenuate signals of 1 Hz and 10 Hz by -30 dB:

In[]:=

tfm[filter_,s_:s]:=filter[{"Bandpass",{1,2,5,10},{30,1}},s]//N//TransferFunctionExpand//Chop

With this, I obtain the corresponding Butterworth and Chebyshev 1 bandpass filters:

In[]:=

{brw,cbs}=tfm/@{ButterworthFilterModel,Chebyshev1FilterModel};

Out[]=

Butterworth filter

81.

5088.47+4723.05s+4227.32

+2012.81

+824.694

+201.281

+42.2732

+4.72305

+0.508847

Chebyshev 1 filter

13.2653

1000.+296.502s+411.457

+72.5658

+41.1457

+2.96502

The Bode magnitude plot attests that frequencies between 2 Hz and 5 Hz will be passed through, and frequencies outside this range will be attenuated, with an attenuation of around -35 dB at 1 Hz and 10 Hz:

In[]:=

{bMagPlot,bPhasePlot}=BodePlot[{brw,cbs},{0.5,18},PlotLayout"List",GridLinesTable[{{1,2,5,10},Automatic},{2}],FrameLabel{{"frequency (Hz)","magnitude (dB)"},{"frequency (Hz)","phase (deg)"}},FrameTicksTable[{{Automatic,Automatic},{{1,2,5,10},Automatic}},{2}],ImageSizeMedium,PlotStyle{cB,cC,cI},PlotTheme"Detailed",PlotLegends{"Butterworth","Chebyshev 1"}];

In[]:=

bMagPlot

Out[]=

	Butterworth
	Chebyshev 1

The phase plot shows that at around 3 Hz for the Butterworth filter and 2 Hz for the Chebyshev, there will be no phase shift. For all other frequencies, there will be varying amounts of phase shifts:

In[]:=

bPhasePlot

Out[]=

	Butterworth
	Chebyshev 1

Later, I will put these frequency responses side by side with the filtered response from the microcontroller to check if they add up.

For now, I will simulate the response of the filters to a signal with three frequency components. The second component lies in the bandpass range, while the other two lie outside it:

In[]:=

inpC=Sin[0.75t]+Sin[2.5t]+Sin[10t];

The responses have a single frequency component, and the two frequences outside the bandpass have in fact been stripped away:

In[]:=

outC=Table[OutputResponse[sys,inpC,{t,0,15}],{sys,{brw,cbs}}];Plot[{outC,inpC},{t,0,15},PlotLegends{"Butterworth","Chebyshev 1","Input"},PlotStyle{cB,cC,cI},PlotTheme"Detailed"]

Out[]=

	Butterworth
	Chebyshev 1
	Input

I then discretize the filters and put them together into one model:

In[]:=

sp=0.1;

In[]:=

StateSpaceModel[ToDiscreteTimeModel[#,sp]]&/@{brw,cbs};sysD=NonlinearStateSpaceModel[SystemsModelParallelConnect[Sequence@@%,{{1,1}},{}]]/.Times[1.`,v_]v//Chop

Out[]=

x. 1	x. 2
x. 2	x. 3
x. 3	x. 4
x. 4	x. 5
x. 5	x. 6
x. 6	x. 7
x. 7	x. 8
x. 8	u. 1 -0.405815 x. 1 +3.43349 x. 2 -12.9245 x. 3 +28.257 x. 4 -39.2343 x. 5 +35.4225 x. 6 -20.3071 x. 7 +6.75879 x. 8
x. 9	x. 10
x. 10	x. 11
x. 11	x. 12
x. 12	x. 13
x. 13	x. 14
x. 14	u. 1 -0.750703 x. 9 +4.44058 x. 10 -11.2287 x. 11 +15.5191 x. 12 -12.3602 x. 13 +5.37913 x. 14
	0.000574178 u. 1 +0.000341168 x. 1 +0.00197143 x. 2 -0.00971768 x. 3 +0.0162245 x. 4 -0.0190824 x. 5 +0.0203388 x. 6 -0.0139566 x. 7 +0.00388075 x. 8
	-0.00131301 u. 1 +0.0022987 x. 9 -0.00583054 x. 10 +0.0108043 x. 11 -0.0203768 x. 12 +0.0201682 x. 13 -0.00706288 x. 14

0.1

For good measure, I simulate and verify the response of the discretized model as well:

In[]:=

With[{tmax=15},With[{inpD=Table[inpC,{t,0,tmax,sp}]},ListLinePlot[Join[OutputResponse[sysD,inpD],{inpD}],PlotLegends{"Butterworth","Chebyshev 1","Input"},PlotStyle{cB,cC,cI},PlotTheme"Detailed",PlotMarkers{Automatic,Scaled[0.025]}]]]

Out[]=

	Butterworth
	Chebyshev 1
	Input

And I wrap up this first stage by deploying the filter and setting up the Arduino to send and receive the signals over the serial pins:

In[]:=

ℳ=MicrocontrollerEmbedCode[sysD,<|"Target""ArduinoNano","Inputs""Serial","Outputs"{"Serial","Serial"}|>,"/dev/cu.usbserial-A106PX6Q"]

Out[]=

MicrocontrollerCodeData

Status:

Target: ArduinoNano



The port name /dev/cu.usbserial-A106PX6Q will not work for you. If you are following along, you will have to change it to the correct value. You can figure it out using Device Manager on Windows, or by searching for file names of the form /dev/cu.usb* and /dev/ttyUSB* on Mac and Linux systems, respectively.

At this point, I can connect any other serial device to send and receive the data. I will use the Device Framework in Mathematica to do that, as its notebook interface provides a great way to visualize the data in real time.

Acquire and Display the Data To set up the data transfer, I begin by identifying the start, delimiter and end bytes:

In[]:=

{sB,dB,eB}=Lookup[ℳ["Serial"],{"StartByte","DelimiterByte","EndByte"}]

Out[]=

{19,44,17}

Then I create a scheduled task that reads the filtered output signals and sends the input signal over to the Arduino, and runs at exactly the same sampling period as the discretized filters:

In[]:=

i1=i2=1;yRaw=y1=y2=u1=u2={};task1=ScheduledTask[If[DeviceExecute[dev,"SerialReadyQ"],If[i1>len1,i1=1];If[i2>len2,i2=1];AppendTo[yRaw,DeviceReadBuffer[dev,"ReadTerminator"eB]];AppendTo[u1,in1〚i1++〛];AppendTo[u2,in2〚i2++〛];DeviceWrite[dev,sB];DeviceWrite[dev,ToString[Last[u1]+Last[u2]]];DeviceWrite[dev,eB]],sp];

I also create a second and less frequent scheduled task that parses the data and discards the old values:

In[]:=

task2=ScheduledTask[yRaw1=yRaw;yRaw=Drop[yRaw,Length@yRaw1];{y1P,y2P}=Transpose[parseData/@yRaw1];y1=Join[y1,y1P];y2=Join[y2,y2P];{u1,u2,y1,y2}=drop/@{u1,u2,y1,y2};,3sp];

Now I am ready to actually send and receive the data, so I open a connection to the Arduino:

In[]:=

dev=DeviceOpen["Serial","/dev/cu.usbserial-A106PX6Q"]

Out[]=

DeviceObject

Class: Serial	ID: 1
Status: Not connected (/dev/cu.usbserial-A106PX6Q)



I then generate the input signals and submit the scheduled tasks to the kernel:

In[]:=

signals[w1,g1,w2,g2];{taskObj1,taskObj2}=SessionSubmit/@{task1,task2};

The input signals are of the form g1 Sin[w1 t] and g2 Sin[w1 t]. The utility function signals generates a cyclically padded list of the sampled values for each input. The code for this and other utility functions are available in the downloadable notebook.

At this point, the data is going back and forth between my Mac and the Arduino. To visualize the data and control the input signals, I create a panel. From the panel, I control the frequency and magnitude of the input signals. I plot the input and filtered signals, and also the frequency response of the filters. The frequency response plots have lines showing the expected magnitude and phase of the filtered signals, which I can verify on the signal plots:

In[]:=

DynamicModule

vars

panel



Since the Arduino needs to be up and running to see the panel update dynamically, I’ll show you some examples of my results here:

Out[]=

2 Hz

3 Hz

4 Hz

3 Hz & 7 Hz

A 2-Hz signal is expected to be attenuated slightly, and this can be seenin both filters' responses, although the Chebyshev appears to be attenuateda bit more.The filtered Butterworth response leads the input signal, while theChebyshev response has a small lag.

Finally, before disconnecting the Arduino, I remove the tasks and close the connection to the device:

In[]:=

TaskRemove/@{taskObj1,taskObj2};DeviceClose[dev]

It’s interesting to see in real time how the response of the deployed filter closely matches that of the analog version.

I think it’s extremely convenient to have the design and microcontroller code so closely coupled, the reason being that if I looked at just the microcontroller code, the coefficients of the filters are abstruse, especially compared to those of the simple first-order lowpass filter. However, the design specifications from which they were obtained are very concrete; thus, having them in one unified place makes it better suited for further analysis and experimentation.

The entire notebook with the filter design, analysis, code generation and visualization is available for download. Of course, you need an Arduino Nano or a similar microcontroller. And with just a few small tweaks—in some cases just changing the target and port name—you can replicate it on a whole host of other microcontrollers. Happy tinkering!

The Utility Functions There are several support functions I created for tasks ranging from setting plot styles, through generating the signals, to parsing the filtered data. I collected them in this separate initialized section to keep the essential workflow focused and smooth.

To deploy the code, the only thing I need is to load the paclet:

In[]:=

Needs["MicrocontrollerKit`"]

I set a color for each of the input signals and filters. Thus, for example, the plot of an input signal and its response on the Bode plot has the same color, and so does the frequency response of a filter and its time response:

In[]:=

{cI1,cI2,cI,cB,cC}=

;

And I created legends based on those colors:

In[]:=

With[{l1=List@@FullForm[ListLinePlot[RandomReal[{-1,1},{5,20}],PlotLegendsAutomatic,PlotStyle{cI1,cI2,cI,cB,cC},PlotTheme{"Business"}]]〚1,2,1〛},lM=Flatten[Cases[l1,HoldPattern[Rule][LegendMarkers,x_]x],1];l=l1〚1,All〛;]

In[]:=

Grid[{{lI=With[{i={1,2,3}},LineLegend[l〚i〛,{"input 1","input 2","input = input 1 + input 2"}〚i〛,LegendMarkers(lM〚i〛)]],lIO=With[{i={4,5,3}},LineLegend[l〚i〛,{"Butterworth filter output","Chebyshev 1 filter output","input = input 1 + input 2"},LegendMarkers(lM〚i〛)]]}}]

Out[]=

	input 1
	input 2
	input = input 1 + input 2

	Butterworth filter output
	Chebyshev 1 filter output
	input = input 1 + input 2

The next set of functions generate the lines showing the magnitude and phase of the filters at a specified frequency:

In[]:=

mLine[w_/;0.5≤w≤18,c_]:=Join[{c},Table[mLine1[w,20Log10[Abs[tf[w]〚1,1〛]]],{tf,{cbs,brw}}]]mLine[___]:={}mLine1[w_,mag_]:=Line[{{Log10[w],-60},{Log10[w],mag},{Log10[0.4],mag}}]pLine[w_/;0.5≤w≤18,c_]:=Join[{c},Table[pLine1[w,phase[w,i]],{i,2}]]pLine[___]:={}pLine1[w_,arg_]:=Line[{{Log10[w],-460},{Log10[w],arg},{Log10[0.4],arg}}]

To generate the lines, I need to compute the phase in degrees:

In[]:=

phase[w_,i_]/;AllTrue[Unevaluated[{cw1,cw2,cw3}],ValueQ]:=With[{arg=Arg[Switch[i,1,brw,_,cbs][w]〚1,1〛]/Degree},Which[i1&&w<cw1,arg+360,i1&&w>cw2,arg-360,i2&&w>cw3,arg-360,True,arg]]phase[w_,i_]:=({cw1,cw2}=getCrossOvers[brw];{cw3}=getCrossOvers[cbs];phase[w,i])

I also need to compute the frequencies at which the phase crosses the negative real axis (180° phase):

In[]:=

getCrossOvers[tfm_]:=Block[{re,im,w},{re,im}=Chop[ComplexExpand[ReIm[tfm[w]〚1,1〛]]];w/.NSolve[re<0&&im0&&w>0,w]]

The next set of functions generate the input signals:

In[]:=

signalGenerator[w_/;w≠0.,g_]:=With[{signal=Table[gSin[wt],{t,0,2.0π/w,sp}]},signal]signalGenerator[w_,g_]:=ConstantArray[g,{20}]

In[]:=

signals[w1_,g1_,w2_,g2_]:=(signal1[w1,g1];signal2[w2,g2];signalLengths;i1=i2=1;)

In[]:=

signal1[args__]:=in1=signalGenerator[args]signal2[args__]:=in2=signalGenerator[args]

I also created a function to determine the signal lengths:

In[]:=

signalLengths:=({len1,len2}=Length/@{in1,in2};lenM=2Max[len1,len2];)

I then use this function use to drop stale values to avoid the plot getting congested:

In[]:=

drop[ls_]:=With[{len=Length[ls]},Drop[ls,Max[len-2lenM,0]]]

And finally, I create a function to parse the incoming data:

In[]:=

parseData[{sB0_,sig1__,dB0_,sig2__}]/;{sB0,dB0}{sB,dB}:=ToExpression[FromCharacterCode/@{{sig1},{sig2}}]parseData[___]:=Nothing

Cite this as: Suba Thomas, "Deploying Bandpass Filters Using the Wolfram Language Microcontroller Kit" from the Notebook Archive (2019), https://notebookarchive.org/2019-09-b39ngmn