27
$\begingroup$

I would like to measure how fast my 6-bladed fidget spinner spins. To do so, after giving it a hard spin, I gently touch the spinner with a wooden stirring stick to create a buzzing sound which usually last for a minute.

enter image description here

I have recorded and plotted the sound it generates:

enter image description here

How can I automatically generate a list of peak times for the above data? My final goal is to plot revolutions per second as a function of time to show spin decay.

Data

To hear the sound in your Mathematica notebook, run the following code:

Sound[SampledSoundList[
  Flatten@ImageData@Import["https://i.stack.imgur.com/qHpp6.png"], 
  22050]]

This will download the following image, turn it into an array, and finally, convert it to a sound object.

enter image description here

$\endgroup$
  • $\begingroup$ Can't you take the Fourier transform of the signal? $\endgroup$ – anderstood Dec 20 '17 at 4:29
  • $\begingroup$ I did, but I didn't detect a significant signal. I was mainly looking in the range of 0 to 200 Hz. $\endgroup$ – Miladiouss Dec 20 '17 at 4:31
  • $\begingroup$ Very interesting: I used two methods to measure the rotation frequency: (1) stroboscope effect: with the flickering "light" source being alternatively the screen of my PC and the image frequency of my mobile phone video camera (2) sound, similiar to your idea, but I used the App called "spectroid" on my mobile phone where I can see the frequency directly in graphic form. Hint: you can reach "very high" frequencies of your fidget spinner using a stream of compressed air from a cleaning device (a can similar to a pepper spray anti dog "weapon"). But be careful and use protective glasses. $\endgroup$ – Dr. Wolfgang Hintze Dec 26 '17 at 19:35
16
$\begingroup$

First, import the audio and extract usable data from it:

audio = Sound[
  SampledSoundList[
   Flatten@ImageData@Import["https://i.stack.imgur.com/qHpp6.png"], 
   22050]]

audioDuration = Duration[audio];
audioSampleRate = AudioSampleRate[audio];
data = AudioData[audio][[1]]; 

Second, use PeakDetect to see which points are peaks (= 1) and which points are not peaks (= 0). Find the location of peaks in seconds.

peaks = PeakDetect[data, 150, 0.0, 0.4];
peakPos = 1./audioSampleRate Position[peaks, 1] // Flatten;
Length[peakPos]

The period of the spinner is the separation between the beats (peaks) times the number of blades:

periods = 6 (peakPos[[2 ;; -1]] - peakPos[[1 ;; -2]])/1

Spin rate, that is revolutions per second, is reciprocal of the period:

spinRates = 1/periods;(* Revolutions per second *)

Convert the data into a list of {time, spin rate} and plot it:

spinRateVStime = 
  Table[{i audioDuration/Length[spinRates], spinRates[[i]]}, {i, 
    Length[spinRates]}];

enter image description here

As it can be seen, the spinner spins 6 times per second and eventually comes to a stop after 12 seconds.

Details

The parameters for PeakDetect needs to be adjusted. To do so, you need to reduce the amount of data to speed up the process, and plot PeakDetect on top of the data and look for a good agreement.

data = AudioData[audio][[1]][[800 ;; 11111]];
peaks = PeakDetect[data, 150, 0.0, 0.4];
ListLinePlot[{data , peaks}, PlotRange -> {All, {0, 1.1}}]

enter image description here

$\endgroup$
  • $\begingroup$ data is just 0 when I run your code. I think it should be data = AudioData[audio]; at least in MMA 11.0. $\endgroup$ – anderstood Dec 28 '17 at 1:02
  • 1
    $\begingroup$ Nice answer! I just think in calculating periods you are reinventing Differences. $\endgroup$ – The Vee Jan 7 '18 at 9:53
9
$\begingroup$

Using Fourier Discrete Transform.

Let's start by observing the Fourier content of the beginning of the signal:

data = Flatten@ImageData@Import["https://i.stack.imgur.com/qHpp6.png"];
fs = 22050; (* sampling frequency *)

data1 = data[[;;20000]];
fourierAbs = Abs[Fourier[data1 - Mean[data1]]];
ListLinePlot[fourierAbs, PlotRange -> {{0, 200}, Full}]

enter image description here

What you see is the left peak (the frequency where looking for) and its harmonics (of multiple frequency). To get the $x$ axis in frequency, you need to multiply the value by $\Delta f$, which is the ratio between the sampling frequency fs and the number of points (20000 here). Note that I substracted the mean of the signal to avoid having a peak in 0.

Now, it is easy to find the peak and its corresponding frequency. Let us wrap that in a function peak:

peak[data_] := Module[{},
  fourierAbs = Abs[Fourier[data - Mean[data]]];
  posPeak = First@First@Position[fourierAbs, 
      Max[fourierAbs[[1 ;; Round[Length@fourierAbs/2]]]]];
  N@posPeak*fs/Length@fourierAbs]

Then it suffices to use a moving map. Since I did not managed to use MovingMap with a step, I did it by myself.

freqs = Table[peak[data[[i ;; i + 20000]]], {i, 1, Length@data - 20000, 1000}];
ListLinePlot[freqs/6, DataRange -> Length@data/fs]

enter image description here

I divided by 6 to count in revolutions per second and compare to Miladiouss's answer. The peaks are due to high harmonics content that you can hear in the end. The curve steps corresond to the frequency resolution $\Delta f$.

Assuming the speed cannot increase, you can improve the above code by limiting the peak search in the range [0,previous peak]: simply initialize posPeak=Round[Length@data/2] and change to this line in peak:

posPeak = First@First@Position[fourierAbs, Max[fourierAbs[[1 ;; posPeak]]]]

$\endgroup$
5
$\begingroup$

You could use the simple noob solution:

snd = Sound[SampledSoundList[
   Flatten@ImageData@Import["https://i.stack.imgur.com/qHpp6.png"], 22050]];
Spectrogram[AudioTrim[snd, {0.1, .4}]]

Mathematica graphics

Here I count 10 hits in the very beginning of your sample. This means, for this short time interval of 0.3s we have

.3/10*6
(* 0.18 *)

about 180ms per round (6 blades) which makes about 330 rounds per minute.

Looking at a sample from around second 9 in your audio and you are down to 50 rpm.

$\endgroup$
  • 3
    $\begingroup$ Nice hair style.. $\endgroup$ – yode Dec 20 '17 at 8:18
  • 1
    $\begingroup$ You are counting. This needs to be automated. $\endgroup$ – Miladiouss Dec 20 '17 at 8:19
2
$\begingroup$

This seems to be a very nice application for wavelets, because the frequency content changes with time. Unfortunately I'm not a wavelet guru, but here is an example, using a lowpass filter, a decimation, an this answer.

data = Flatten@ImageData@Import["https://i.stack.imgur.com/qHpp6.png"];
fs = 22050;(*sampling frequency*)
data2 = LowpassFilter[data, .01];
data3 = data2[[1 ;; -1 ;; 10]];
sampleRate = fs/10;
cwd = ContinuousWaveletTransform[data3, DGaussianWavelet[7], {11, 12}, 
   SampleRate -> sampleRate];
freq = (cwd[
       "SampleRate"]/(#*
        cwd["Wavelet"]["FourierFactor"])) & /@ (Thread[{Range[
        cwd["Octaves"]], 1}] /. cwd["Scales"]);
ticks = Transpose[{Range[Length[freq]], 1/6*freq}];

WaveletScalogram[cwd, Frame -> True, 
 FrameTicks -> {{ticks, Automatic}, Automatic}, 
 FrameLabel -> {"Time", "Spin rate"}, ColorFunction -> "RustTones"]

The output displays the same curve as Miladiouss's answer. The ridge still has to be extracted though. Hopefully someone can improve this answer and mark it as a community wiki.

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.