I’ve been working on some instrument design projects and have hit a brick wall of sorts. My prototypes are riddled with noise, most likely 60 hz. My thought here was to learn a bit more about signal processing to (a) see if I can get a better understanding of what’s going on and (b) see if this is a possible project for students.

So the setup is as follows. I’ve got an Arduino microcontroller that does one of two things, it either reads the signal from a noisy light detector (in this case, an LED connected to an op amp in a current-to-voltage configuration) or – for debugging purposes – outputs a fixed signal frequency by printing $A0 + A cos(2 Pi f millis()/1000)$ where $A0$ and $A$ are amplitude offset and signal amplitude, respectively, $f$ is the frequency and since millis() returns a value in milliseconds, it is divided by 1000. To enact a sampling rate, I set a delay(dt) in the loop routine where dt is the delay time in milliseconds.

On the *Mathematica* side, it’s pretty easy to read the serial data from the Arduino with d = DeviceOpen["Serial", {<port>, "BaudRate"-><baudrate>}] with replacing <port> and <baudrate> with your values. The code below is a tad clunky, but works well at grabbing data and converting it into a format that Mathematica wants.

dt = 0.075;
data = ToExpression /@ 
   Flatten@ImportString[
     FromCharacterCode@
      DeleteCases[DeviceReadBuffer[d], Alternatives @@ {13}], "CSV"];
time = Table[(n - 1) dt, {n, Length@data}];

Parsing the 2nd line that sets data, we see the command reads all of the data in the serial buffer with DeviceReadBuffer, which returns a bunch of ASCII codes. We get rid of ASCII code 13 to simplify double spacing problems in our data and then convert the remaining data into numeric strings with FromCharacterCode. Then, we use ImportString to convert the string of numbers into a list, Flatten it because there’s a bunch of nested lists in there, and then finally convert those numeric strings into values with ToExpression.

*Note*: When bouncing between Mathematica and the Arduino IDE, it’s important to disconnect the Arduino from Mathematica via DeviceClose.

I start with a very simple situation where my delay (dt) is set to 100 ms and I output a 2 hz signal.

What I want to do now is find the frequency of this signal using a Fourier transform. This question on Stack Exchange led me through much of the process. Here’s the code for processing the data.

working = data;
ft = Fourier[working, FourierParameters -> {-1, -1}];
freq = Table[N@(n - 1)/dt/Length@ft, {n, Length@ft}];
dtp = Transpose[{freq, Abs[ft]}][[1 ;; Round[Length@freq/2]]];
ListPlot[{dtp, 
  Transpose[{#1, PeakDetect[#2, 1] #2}] & @@ (Transpose@dtp)}, 
 PlotRange -> {{0, Max@freq/2}, All}, Joined -> {True, False}, 
 Filling -> {2 -> 0}, FillingStyle -> Opacity[1]]

In the code above ft is the Fourier transform with parameters that give frequency amplitudes that are 1/2 the signal amplitude. The data to plot (dtp) is just the first 1/2 of the transform to keep things simple and I’ve added a second trace in the plot to pull out the peaks just for fun.

The main frequency observed is pretty close to 2 hz and the first point is 129, which is pretty close to the offset of the amplitude (128). I played around with a couple frequencies (including fractional ones) and the results were as expected. I did stay under the Nyquist frequency (1/2 the sampling rate or 5 hz) for these exercises, but wanted to see what would happen when the signal was *higher* than the Nyquist frequency and aliasing would occur. So I set the frequency to 22 hz and here’s the FT:

I wanted to know if the 2 hz response is expected, so after a little bit of web searching, I found this site (among others) that shows how to estimate the alias frequency $f_{alias}$ given the signal frequency $f_{signal}$ and the sampling frequency $f_{sampling}$. I rearranged the equation a bit to get this: $ 2 \times Abs(f_{sampling} \times N – f_{signal}) < f_{sampling}$. That doesn’t give us the alias frequency, but the integer $N$ that satisfies that inequality will give us the alias frequency with $f_{alias} = Abs(f_{sampling} \times N – f_{signal})$. With the numbers in the above figure, we get $N = 2$ and $f_{alias} = 2$ which is reasonably close. This system can be solved using Mathematica fairly easily:

Module[{fsamp = 10, fsig = 22, tmp, n}, 
 tmp = FindInstance[2 Abs[fsamp n - fsig] < fsamp, n, Integers];
 Abs[fsamp n - fsig] /. First@tmp]

We can check the code by changing the sampling frequency. I set the delay to 75 ms (so $f_{samp} = 1/0.075 = 13.33 hz$) and kept the signal the same. The aliased frequency is predicted to be 4.6 hz.

The frequency of the peak in the FT spectrum is 4.5 hz, which is pretty close to the predicted value.

Now I wanted to move on to real data. Again, what I have here is an LED configured as a light sensor, and in its current configuration, the signal is dominated by noise. I’m fairly certain that the noise is 60 hz, but in order to detect that signal with the Arduino/Serial setup, I would need a sampling frequency of at least 120 hz, corresponding to a delay of about 8 ms. I am concerned that the analog read and the serial communication will add too much uncertainty to the timing. Keeping the sampling frequency at 13.3 hz and expecting the noise to be 60 hz, I anticipate the aliased signal will have a frequency of 6.5 hz.

…and sure enough, it does. Note that there’s an additional signal there at 0.5 hz, and at present, I don’t know its source. Is it an overtone? Is it another (aliased) signal in the data? Is it real? More investigations to follow, but if you happen to have an idea, feel free to mention it in the comments.