FFT Signal Processing in Mathematica

Let’s create and plot some data.

We’ll first make 1024 evenly spaced time values, 120 μs apart, which is about the closest time spacing one can easily get sampling analog voltages with an Arduino.  We won’t simulate any noise in the data.  We’ll add two frequencies,  440 Hz and 1400 Hz, but the 1400 Hz part will have 10 times the amplitude as the 440 Hz part.  The first plot shows all the data, and the next plot shows the details in a smaller time window.  Note that there are several data points in any recognizable oscillation.  The last few commands in this section create and play a sound file of the data.  Since it’s such a short file, the code loops, and you control the loop with the button.  Note that what you mostly hear is the 440 Hz tone (concert A), although your ears may also hear some of the smaller 1.4 kHz part.

tData=Table[i 0.00012,{i,0,1023}];
yData=Table[ 0.1 Sin[2 Pi 440 i 0.00012]+Sin[2 Pi 1400 i 0.00012],{i,0,1023}];
plotData=Partition[Riffle[tData,yData],2];
ListPlot[plotData, Joined->True, Frame->True]
ListPlot[plotData, Joined->True, Frame->True,PlotRange->{{0,0.02},All},Epilog->Point[plotData]]
sound=ListPlay[yData,SampleRate->8333];
DynamicModule[{state="pause"},Grid[{{Toggler[Dynamic@state,{"pause"->Style["▸",24],"play"->Style["◼",16]}],Dynamic@If[state=="play",If[ScheduledTasks[]=={},RunScheduledTask[EmitSound[sound]]];"playing",If[ScheduledTasks[]≠{},RemoveScheduledTask[ScheduledTasks[]]];"paused"]}}]]

Now take the Fourier transform and analyze the results

ft=Fourier[yData];

Note that the Fourier transform returns complex numbers (the complex numbers contain the phase relationships between the various frequency contributions).  So in order to plot the Fourier components, use Abs[ ] to get the magnitude of each component.

ListPlot[Abs[ft],Joined->True, PlotRange->All]

Now find the indices that correspond to the peaks shown in the Fourier spectrum.  First, to decode the indices:  the first index corresponds to the “zero” frequency, which is the DC or average value of the signal being analyzed (which should be near zero in this case, but won’t be in the lab).  The next few indices correspond to (Subscript[f, 0])=1/(Subscript[t,  total]), then 2(Subscript[f, 0]), 3(Subscript[f, 0]), etc.  Near the middle, right after (Subscript[f,  max]), the indices correspond to negative frequencies, counting down to -(Subscript[f, 0])at the right end of the spectrum.  The next command finds the peaks.  The three parameters (Automatic, Automatic, 1) tell FindPeaks what kind of peaks to look for -- the “1” is the minimum height to look for.

peaks=FindPeaks[Abs[ft],Automatic,Automatic,1]

Now we calculate the lowest or fundamental frequency, then use the indices of the peaks to calculate the frequencies of the low and the high peak.

f0=1./tData[[Length[tData]]]

fLo=54 f0
fHigh=172 f0

Filtering the signal

What would the signal look or sound like if we removed the low frequencies (both positive and negative just for the sake of completeness)?  First, we zero the Fourier amplitudes near the low frequency peaks, then re-create the flitered signal.  Note that we plot the real part of the re-created function, since the InverseFourier command, like the Fourier command, produces a complex result.  Note the difference in the sound of the recreated file.

Do[ft[[i]]=0,{i,50,60}];
Do[ft[[i]]=0,{i,966,976}];
filteredSignal=InverseFourier[ft];
replotData=Partition[Riffle[tData,Re[filteredSignal]],2];
ListPlot[replotData, Joined->True, Frame->True,PlotRange->{{0,0.02},All},Epilog->Point[replotData]]

filteredSound=ListPlay[Re[filteredSignal],SampleRate->8333];
DynamicModule[{state="pause"},Grid[{{Toggler[Dynamic@state,{"pause"->Style["▸",24],"play"->Style["◼",16]}],Dynamic@If[state=="play",If[ScheduledTasks[]=={},RunScheduledTask[EmitSound[filteredSound]]];"playing",If[ScheduledTasks[]≠{},RemoveScheduledTask[ScheduledTasks[]]];"paused"]}}]]