Vernier Sensors in the Wolfram Language

Armeen Mahdian and Ian Johnson

Join the Conversation #WolframTechConf

Finding and Initializing Devices

Find all Vernier devices.

FindDevices["Vernier"]

Open the first Vernier device found.

forceSensor=First[DeviceOpen[FindDevices["Vernier"]]]

We can use DeviceConfigure to modify the MeasurementPeriod of the device.

DeviceConfigure[forceSensor,"MeasurementPeriod"->(Quantity[1, "s", seconds, "Seconds"])]
DeviceConfigure[forceSensor,"MeasurementPeriod"->0.05]
DeviceConfigure[forceSensor,"MeasurementPeriod"->1]
DeviceConfigure[forceSensor,"MeasurementPeriod"->Automatic]
DeviceConfigure[forceSensor,"MeasurementPeriod"->Default]

Reading From Devices

Read a single value.

DeviceRead[forceSensor]

Read the entire buffer and plot it.

bufferValues=DeviceReadBuffer[forceSensor];
ListLinePlot[bufferValues]

Dynamically view values being read from the sensor.

DynamicModule[{append, readings},
append=False;
readings={};
Row[{Column[{
Button[Dynamic[If[append,"Stop","Start"]],
If[append == False,
readings={};baseTime=AbsoluteTime[];append=!append;,
append=!append;
];,
ImageSize->{100,180}
],
Button["Keep",keptValues=readings,ImageSize->{100,180}]
}]
,"   ",
Dynamic[If[append,readings = Append[readings,{AbsoluteTime[]-baseTime,DeviceRead[forceSensor]}]];myListLinePlot[readings]
]
}]
,
Initialization:>{
forceSensor = FindDevices["Vernier"][[1]];
DeviceConfigure[forceSensor,"MeasurementPeriod"->0.005];
DeviceOpen[forceSensor];
myListLinePlot[data_]:=ListLinePlot[
data,
AxesStyle->Thick,
PlotStyle->Directive[Darker[Green],Thick],
AxesLabel->{"s",If[Length[data]>0,data[[1,2,2]],"No Data"]},
ImageSize->{600,400},
LabelStyle->Large
]
}
]

View kept values.

keptValues

Manual Fitting

Reading

Empty the buffer 30 times to gather some data. Remove the units to speed up further calculations.

DeviceRead[forceSensor];
vals=Last[Reap[Do[Pause[1];Sow[DeviceReadBuffer[forceSensor]],{30}]]];
data=Flatten[QuantityMagnitude[vals]];

Manipulate and find what portion of the data you are interested in.

Manipulate[ListLinePlot@Take[data,{i,j}],{{i,7773-4500,"Lower"},1,j,1},{{j,7773-3000,"Upper"},i,7773,1},Button["Set Values",{selection,lower,upper}={Take[data,{i,j}],i,j}],SynchronousUpdating->False]

Finding Peaks and Troughs

Visualize the peaks on the graph.

Dynamic[ListLinePlot[
Take[selection], 
Epilog-> {Red,PointSize[0.02],Point[peaks=FindPeaks[selection]]}
]]

Then we can flip our selection to find the troughs.

ListLinePlot[{selection,-selection}]

Highlight the peaks on the flipped graph.

ListLinePlot[{selection,-selection},Epilog-> {Red,PointSize[0.02],Point[FindPeaks[-selection]]}]

Now highlight the troughs on the original graph.

ListLinePlot[{selection,-selection},Epilog-> {Red,PointSize[0.02],Point[troughs=Transpose[Transpose[FindPeaks[-selection]]*{1,-1}]]}]

Combine both for peaks and troughs on a single graph.

ListLinePlot[
selection, 
Epilog-> {Red,PointSize[0.02],Point[Join[troughs,peaks]]}
]

Frequency and Period

Find the periods by calculating the differences between the peaks.

periods=Differences[peaks][[All,1]]

Calculating the Mean of these gives us an average period for our selected portion.

N[averagedPeriod=Mean[periods]]

However, we still need to multiply this by the time between each data point in order to calculate the frequency from the period.

deltat=(1/200);
N[initialFrequency=(averagedPeriod*deltat)]

Now in order to convert this frequency into radians per second, we divide 2π by this number.

N[frequency=(2π/initialFrequency)]

Amplitude and Displacement

Now we can find the average amplitude by subtracting the average minimum Y value from the average maximum Y value, and dividing by 2.

peakAverage=Mean[peaks[[All,2]]]
troughAverage=Mean[(troughs)[[All,2]]]
amplitude=((peakAverage-troughAverage)/2)

Additionally, we can find the displacement from the X axis by finding the average between the peak average and the trough average.

displacement=((peakAverage+troughAverage)/2)

With all of this, we can combine everything into a single equation:

Plot[amplitude Sin[N[frequency] t]+displacement,{t,1,2000 deltat}]

Comparing and Modifying Phase Shift

Making sure that we have the same number of points as our original signal with sampling, we can plot our equation alongside the original.

range=upper-lower;
ListLinePlot[{Table[amplitude Sin[frequency t]+displacement,{t,0,range*deltat,deltat}],selection}]

This, however, does not account for the phase shift, so we can find the differences between the X values of the peaks.

FindPeaks[Table[amplitude Sin[frequency t]+displacement,{t,0,range*deltat,deltat}]][[All,1]]-peaks[[All,1]]

So by averaging these, we get the average time step difference between the two signals.

Mean[%]

Now we have to multiply that by the actual sampling frequency and then divide 2π by that value to return a phase shift in radians.

phaseShift=N[2π/(%)]

Accounting for the phase shift, we can graph the signals again:

ListLinePlot[{Table[amplitude Sin[frequency t+phaseShift]+displacement,{t,0,range*deltat,deltat}],selection}]

Closer, but not perfect. Let’s try using Manipulate to see if we can get any closer.

Manipulate[
ListLinePlot[{Table[amplitude Sin[frequency t+phaseShift+ϕ]+displacement,{t,0,range*deltat,deltat}],selection}],{{ϕ,-1},-2,2}]

Automatic Fitting

Now we can use built-in functions to fit the data that we got from the sensor.

First we assume we are at our maximum 200 readings per second value.

samplingRate=200;

maxRange=Floor[N[((Length[data])/samplingRate)]]

Now we have our underdamped oscillator formula for a better fit.

solx[a_,γ_,ω0_,δ_][t_] :=  a Exp[-γ t]Cos[ω0 t + δ]

We will select a range that we would like to try to fit.

myRange=samplingRate{0,maxRange}+1;
Length[d0=Transpose[{Range[0,-Subtract@@myRange]/samplingRate,QuantityMagnitude[data[[Span@@myRange]]-Mean[data]]}]]

nlm=NonlinearModelFit[d0,solx[a,γ,ω0,δ]''[t],{a,{ω0,12},δ,γ},t];
nlm[t]
nlm["ParameterTable"]

Show[{Plot[nlm[t],{t,0,11},ImageSize->800],ListPlot[d0, PlotStyle->Red]}]