Frequency of ticks in mechanical watches

Question

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

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:

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"]

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]

Thanks for the help! :)

Answer 1

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]}
]

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:

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.]]

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.

Answer 2

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}]

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}]

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"]}]

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]

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]

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.