9
$\begingroup$

I'm looking to find the frequency of a relatively large dataset (or a relatively short audio file) of a repeating sound, like so:

sound_file Here are the original files, for those willing to play. :) eta2412 data1 Other examples: phenix140sc chaika1601

Blue is the original soundwave, and orange is a Highpass and Wiener filtered form at 44100Hz sample.

A finely zoomed in picture of each "tick" (unfiltered in an audio program) looks like so: enter image description here

It's extremely oscillatory, however, there is a clear and obvious local maximum (when zoomed in far enough anyways (the duration of a tick is 0.016s or so)

I'm not particularly interested per say in the amplitude, however what I'd like to do, is to be able to automate and count the freq. of the "beats/ticks". This is currently done 'by hand' for each individual file I have to deal with and when a person has MMA...well, you get the point

Accuracy is a must in this case, the possible periods of the beats are at 0.200 [s] and 0.166 [s] repeating (5hz, and 6hz respectfully). Any kind of drift more than +/-0.002 second isn't acceptable.

Through many readings of this link this link and many others, I've managed to cobble together something...that doesn't seem to get me anywhere.

dir = NotebookDirectory[];
SetDirectory[dir];
file = "eta2412.m4a";
snd = Import[file];
sndhpf = HighpassFilter[snd, 20000];
sndwf = WienerFilter[HighpassFilter[snd, 20000], 25];
sndSampleRate = Import[file, "SampleRate"];
p1 = AudioPlot[snd,  PlotRange -> {{0, 5}, {-1*10^-2, 1*10^-2}}];
p2 = AudioPlot[sndwf, PlotRange -> {{0, 5}, {-1*10^-2, 1*10^-2}}, PlotStyle -> Purple];
p3 = AudioPlot[sndhpf, PlotRange -> {{0, 5}, {-1*10^-2, 1*10^-2}}, PlotStyle -> Green];
data = Flatten[AudioData[sndwf]];
data1 = Drop[data, {50000, Length[data]}];

I've attempted to do audio processing itself,

getting to

res = AudioLocalMeasurements[sndwf,"RMSAmplitude",PartitionGranularity ->{Quantity[0.1, "Milliseconds"],Quantity[0.1, "Milliseconds"]}];
ListLinePlot[res, PlotStyle -> Red, PlotRange -> {{0, 2}, {0, 0.005}}]

Show[AudioPlot[sndwf, PlotRange -> {{0, 5}, {0, 0.005}}],ListLinePlot[res, PlotStyle -> Red, PlotRange -> {{0, 5}, {0, 0.005}}], AspectRatio -> 3/6, ImageSize -> "Large"]

show

Now, I'm unfortunately at a complete loss on how to get each ticks local peak and calculate their global frequency, or rather, period.

How would one go about this? I've tried in Fourier as suggested in other posts..but the noise seems too extreme to give me anything I can understand or use to get further.

ListLinePlot[p = PeriodogramArray[data1], PlotRange -> All, ImageSize -> Large]

fourier

Thanks for the help! :)

$\endgroup$
  • $\begingroup$ How noisy are other clips? Is there going to be more interference? $\endgroup$ – Carl Lange Jul 29 '18 at 23:37
  • 1
    $\begingroup$ Ahh, all recordings are in the same place, so the only extra noise/interference will come from the watch movement themselves, however, they mostly look all the same, I"ll add other examples! $\endgroup$ – morbo Jul 29 '18 at 23:49
  • $\begingroup$ Great! Also,you say it's either 5hz or 6hz, are those the only two options? (eg are you classifying watch ticks into those two boxes) Or are you trying to see if any of your watches are 8hz? Also can you post the manual tick count of your examples so we don't have to all do it :) $\endgroup$ – Carl Lange Jul 29 '18 at 23:55
  • 1
    $\begingroup$ Ahh, all the watches i deal with for work are 5hz or 6hz explicitly, to gauge their accuracy, or (where the research is) the degradation of them over time, I load up the recorded file, and measure the peaks by hand in sonic visualiser and average the peak period. ranfft.de/cgi-bin/… as an example or ranfft.de/cgi-bin/… $\endgroup$ – morbo Jul 30 '18 at 0:01
  • 1
    $\begingroup$ The Phenix has a period of 0.187, making it 96 seconds per day fast eta2412 has a period of 0.16668, making it -0.112 seconds slow. chaika1600 has a period of 0.1678, making it -9 second slow. As examples as of today hope that helps $\endgroup$ – morbo Jul 30 '18 at 0:10
6
$\begingroup$

Third, much better answer

This is a bit easier to understand graphically, so I put a Manipulate together.

a = filterAudio[First@audios];
Manipulate[
 ListLinePlot[Rescale[Ramp[First@AudioData[#]]], 
    PlotRange -> Full] & /@ RandomChoice[AudioPartition[a, x], 6],
 {x, N[1/7], N[1/4]}
]

Image of the Manipulate

What we're doing here is splitting the audio file into partitions of length x. We can see that when the peaks (the ticks) line up correctly, x is our tick period.

If we Rescale the data, we always get some point at 1. Thus we don't need to use FindPeaks at all.

Now that we've established some understanding of how that works, we can solidify it a little. What do we mean by "when the peaks line up correctly"? Essentially, we mean what the difference between the position of 1. is between each clip. So when the difference between each of those positions is close to 0, we've found a a period where the peaks line up.

Now, let's get a list of those values for different potential periods.

Table[{x, 
  Mean[Differences[
    Flatten[First@
        Position[#, 1.] & /@ (Rescale[Ramp[First@AudioData[#]]] & /@ 
        Take[AudioPartition[a, x], {3, 12}])]]]}, {x, 0.16, 0.17, 
  0.0001}]

(Modify the 0.0001 for more or less precision).

We can plot this also, as a helpful mental model of what's going on:

Plotted data

NOW! The great thing about this is that when there's no (or very little) difference, we must have a period where the ticks line up! That is, when the mean difference crosses the x-axis.

The icing on this cake is that we can use the points to generate a LinearModel to get the data, then find out when it crosses the axis!

solvePeriods[aud_] := Module[{a, model},


a = filterAudio[aud];
  model = 
   LinearModelFit[
    Table[{Mean[
       Differences[
        Flatten[First@
            Position[#, 
             1.] & /@ (Rescale[Ramp[First@AudioData[#]]] & /@ 
            Take[AudioPartition[a, x], {3, 12}])]]], x}, {x, 0.16, 
      0.17, 0.005}], v, v];
  model[0.]]

Solved 0

There we go! a (the audio file) in this context was eta2412, so 0.166776 is looking good to me. Hope this is helpful! You'll need to modify the bounds for x above, depending on your initial guess of the period.


Second bad answer

averageTicksPerSecond[aud_] := Mean[Module[{data, peaks},
 data = Rescale[Ramp[First@AudioData[#]]];
 peaks = FindPeaks[data, 150, 0.00, 0.3];
 Length@peaks
 ] & /@ AudioPartition[filterAudio@AudioTrim@aud, Quantity[1, "Seconds"]]]

We split the audio into 1 second partitions, use FindPeaks to find and count the peaks in that small partition and average the number of found peaks across the whole audio file to find our average ticks per second. Then convert that into the period like so:

N[1/averageTicksPerSecond[aud]]

This gives me a result of a period of 0.164179 for the eta2412 (which is still off from 0.16668 but it's not awful) and 0.159383 for the chaika (from 0.1678).

FindPeaks takes a lot of magic number arguments though, I don't like this answer much either. It does feel a little less hacky though. It would be nice to be able to inform FindPeaks of roughly the number of expected peaks (eg, the hz) to help it automatically guess parameters. Maybe someday :)

side note: you'll find this helpful if you're playing with FindPeaks:

Module[{data, peaks},
 data = Rescale[Ramp[First@AudioData[#]]];
 peaks = FindPeaks[data, 150, 0.00, 0.3];
 ListLinePlot[data, Epilog -> { Red, PointSize[0.04], Point[peaks]},
  PlotRange -> Full]
] & /@ AudioPartition[filterAudio@AudioTrim@aud, Quantity[1, "Seconds"]]

First, bad answer

OK, so the following is very hacky and doesn't work well, but it gives a crap estimate.

Let's set up a function to filter the audio (same way you did it).

filterAudio[aud_] := HighpassFilter[WienerFilter[aud, 10], 44100]

Now, let's write a function to count the number of ticks. (Some people may want to look at FindPeaks - I think that's probably very useful but I couldn't make it happen).

Let's take the first channel of the AudioData,Ramp it to only get >0 values, and Threshold it to throw away anything too low. 0.0004 is a magic number here and will want to change per audio file.

threshold[aud_] := 
  Module[{ad = First@AudioData[aud]}, 
    Threshold[Ramp[ad], Threshold[ad, {"Hard", 0.0004}]]]

Then, use SequenceSplit to get the number of sequences (ticks) surrounded by long lists of zeros (silence).

numberOfTicks[aud_] :=
  Length[SequenceSplit[threshold[aud], Table[0., 1000]]]]

Turning that number of ticks into the period is easy, so to tie it together...

period[aud_] := 1/(numberOfTicks[filterAudio[aud]]/Duration@aud)

And finally:

audios = Import /@ \
{"http://lizard-truth.com/wp-content/uploads/2018/07/eta2412.m4a", 
   "http://lizard-truth.com/wp-content/uploads/2018/07/phenix_140sc.\
m4a", "http://lizard-truth.com/wp-content/uploads/2018/07/chaika.m4a"}

UnitConvert[period /@ audios, "Seconds"]

(* {Quantity[0.164251, "Seconds"], Quantity[0.194472, "Seconds"], Quantity[0.164715, "Seconds"]} *)

That's not enormously far away from your numbers, but definitely far too far away for your needs. However, I thought it may be useful for someone.

$\endgroup$
  • 1
    $\begingroup$ Ahh Those are awesome solutions....though, unfortunately they don't give the desired accuracy. I had found this post this morning, mathematica.stackexchange.com/questions/4272/… Maybe it'll help for an idea. This person managed to find a half decent way to measure BPMs of a song, however, I couldn't apply it to anything useful. I'll play with your second equally bad answer, maybe I can find some better FindPeak numbers :) $\endgroup$ – morbo Jul 30 '18 at 21:13
  • 1
    $\begingroup$ Carl, your third algorithm seems to work great, Im going to do some testing and get back to you. $\endgroup$ – morbo Aug 1 '18 at 11:04
  • $\begingroup$ In case you haven't seen it, there's an answer on the Wolfram Community that looks like it's written by someone who knows what they're doing :) $\endgroup$ – Carl Lange Aug 1 '18 at 11:10
5
$\begingroup$

The data looks like periodic pulses of incoherent noise, familiar to an old pulsar astronomer. So, use a method sometimes used for pulsars. Start by working in intensity rather than amplitude. That's appropriate if the contents of the bursts are incoherent.

intensity[dat_Audio] := AudioData[dat][[1]]^2

Have a look at some intensity data.

eta = Import["eta2412.m4a"];
ListPlot[intensity[eta], {PlotRange -> All, Joined -> True}]

enter image description here

Looks like a nice regular train of irregular pulses, like a pulsar.

The Fourier power density looks good, too.

ListPlot[PeriodogramArray[intensity[eta]][[2 ;; 2000]], {PlotRange -> All, Joined -> True}]

enter image description here

I've plotted elements [[2;;2000]] because the whole plot is just a wad of unresolved stuff. Element 1 is the zero frequency term, large and uninteresting. Note the small peaks between the larger ones: there's really some "tick, tock" going on. There's some feature in the sound that alternates from pulse to pulse. Investigating this is left for the future.

To get an accurate period and a pulse profile, we'll fold and resample the data. First, choose how many resampled points to have per pulse period. This shouldn't be critical: it should be enough to resolve the persistent features in the pulses without oversampling them too much.

samplesPerTick = 1000;

The resampling function periodicDecimate filters and resamples the data so that there are exactly samplesPerTick samples per pulse period.

periodicDecimate[dat_Audio, frequency_Quantity] := ArrayResample[intensity[dat], 
   Scaled[frequency samplesPerTick/AudioSampleRate[dat]], Antialiasing -> True]

Now, take the average signal in modular time, with the modulus being the assumed pulse period. The resampling makes this a simple array manipulation.

fold[dat_Audio, frequency_Quantity] :=
    With[{dd = periodicDecimate[dat, frequency]}, 
       Mean /@ Transpose[ArrayReshape[dd, 
           {Floor[Length[dd]/samplesPerTick], samplesPerTick}]]]

Now, a couple of parameters to restrict the frequency search space to minimize the chance of a spurious result.

fLow = Quantity[5.99, "Hz"];
fHigh = Quantity[6.01, "Hz"];

The folded pulse profile will be weak and broad if the assumed frequency is wrong. It will be strong and narrow if the frequency is right. Its variance is a good measure. Annoyingly, we have to dequantify and quantify to plot it.

Plot[Variance[fold[eta, f Quantity["Hz"]]], 
  {f, fLow/Quantity["Hz"], fHigh/Quantity["Hz"]}]

enter image description here

Note the very sharp peak. Unlike the usual periodogram, this incorporates information from the detailed pulse shape, or, if you prefer to think that way, from the train of harmonics in the spectrum. Find the peak:

fMax[dat_Audio] := 
   f Quantity["Hz"] /. NMaximize[Variance[fold[dat, f Quantity["Hz"]]], 
     {f, fLow/Quantity["Hz"], fHigh/Quantity["Hz"]}][[2]]

fMax[eta]
(* Quantity[5.9991, "Hertz"] *)

Here's the mean pulse profile:

ListPlot[fold[eta, %], PlotRange -> All, Joined -> True]

enter image description here

Now, the Phenix. We have to change the frequency range.

phenix = Import["phenix_140sc.m4a"];
fLow = Quantity[4.99, "Hz"];
fHigh = Quantity[5.01, "Hz"];
fMax[phenix]
(* Quantity[5.00396, "Hertz"] *)
ListPlot[fold[phenix, %], PlotRange -> All, Joined -> True]

enter image description here

Warning: this code uses the sample rate in the audio file. Clock oscillators in audio equipment are often specified to only a part in 10^4. This precision is more than adequate for audio, but may not be enough for your purposes. Frequency drift within spec, especially with temperature, is common, so a single calibration may not suffice.

$\endgroup$
  • $\begingroup$ Very educational, what a great answer! $\endgroup$ – Carl Lange Aug 2 '18 at 23:39
  • $\begingroup$ This is very nice. For what it's worth, the string length method of Lafler and Kinman (1965 paper) shows a drop at 0.166678, corresponding to 5.9995 Hz (thus in reasonable agreement with this analysis). I doubt string lengths is the best method to employ here-- the fact that it gives something in reasonable agreement is encouraging though. $\endgroup$ – Daniel Lichtblau Aug 3 '18 at 16:08

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.