# Accurately calculate tick tock period of mechanical watches?

Posted 6 months ago
1335 Views
|
15 Replies
|
15 Total Likes
|
 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 chaika1601Blue 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 pointAccuracy 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,this (the last link does nearly what but alas it's bpm) 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] I have posted this over at stackexchange as well, incase this seems familiar to anyone. But I figured I'd try to hit two birds with two stones, should someone there post a solution I"ll gladly update it here. Thanks for the help! :)
15 Replies
Sort By:
Posted 6 months ago
 As additional information IAll 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. Some movement info and examples can be found: eta2412 Phenix 140SC As measured yesterday:Phenix 140SC 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. Thanks again for the help!
Posted 6 months ago
 Apparently my Linux machine does not like that format, so I cannot experiment. From the pictures, I wonder if you might use Clip to set values above some threshold to 1, and below to 0. The main frequency might then be more easily obtained from Fourier or similar.
Posted 6 months ago
 Hi Daniel! Ahh here are some links for the audio files, with .mp3 extension....the m4a are typically Apple. I attempted to something similar with MaxPeak, however it seemingly failed with the line plot only returning 1s (likely my poor MMA skills at work here)...The only issue I could see with it however, is that there isn't any nice cutoff for maximum peaks, between different movements things can change wildly. I'll look into trying clip, hope the mp3s are useful :) in the meantime.
Posted 6 months ago
 Mor,As vibration analysis is one of my interests I could not let this one go by! so here is my take on this!I attached the file at the bottom to make it easier to run.Start with your import. dir = NotebookDirectory[]; SetDirectory[dir]; file = "eta2412.m4a"; snd = Import[file]; sndSampleRate = Import[file, "SampleRate"]; Convert the audio to raw data: sndData = AudioData[snd][[1]]; Use peak finding to get the peaks -- I adjusted the values to do well on this type of data -- The goal is to get more points than you need around the peaks to make sure we get all peaks -- err on the side of overinclusiveness peaks = FindPeaks[sndData, 0.9, 0.0025, 0.002]; Length[peaks] ListLinePlot[sndData, Epilog -> {Red, PointSize[0.01], Point[peaks]}, PlotRange -> All, ImageSize -> Large] Now let's change to use binary (1's or zeros) for the peaks to make a nice time domain array binarypeaks = PeakDetect[sndData, 0.9, 0.0025, 0.002]; We need better frequency resolution because you only collected 11 seconds of data -- the longer the data, the better the resolution. We can make up for some of that if we pad the data with a ton of zeros: paddingAmount = 20; binarypeaksPad = ArrayPad[binarypeaks, paddingAmount*Length[binarypeaks]]; Use the approximate frequency to window the data to what we care about approxFreq = 6; bins = Length[binarypeaksPad]; approxbins = Round[approxFreq*(bins/sndSampleRate) *1.2] // N This will tell us our resolution we will get in the frequency domain -- make sure we are good enough resolution = 1/(bins/sndSampleRate) // N Compute a power spectral density of the binary bins and plot it -- we have a spike at DC and a spike at our desired frequency spec = Abs[Fourier[binarypeaksPad]^2]; ListLinePlot[spec[[1 ;; approxbins]], PlotRange -> All, ImageSize -> Large] Find our big peak in this region freqbin = Position[spec[[1 ;; approxbins]], Max[spec[[Round[approxbins/2] ;; approxbins]]]][[1, 1]] - 1 This is our actual frequency In[21]:= actualFreq = freqbin/(Length[spec]/sndSampleRate) // N Out[21]= 5.9996 This is our period In[22]:= 1/actualFreq // N Out[22]= 0.166678 It matches your hand processed results.Regards,Neil Attachments:
Posted 6 months ago
 Mor,One other thought -- I also realized that the technique above is insensitive to the parameters to the PeakDetect. For example, peaks = FindPeaks[sndData, 0, 0.002, 0.002]; Works fine so you should be able to generalize the algorithm to different audio clips -- you may only have to adjust the one number -- the amplitude cutoff for your clip loudness.Did you try what I posted? Is this what you wanted? If so, I will cross post it in StackExchange.Regards,Neil
Posted 6 months ago
Posted 6 months ago
 Mor, Firstly, is there a reason that you didn't filter the audio files? There is a lot of background noise in the recordings. Or was this simply not required, being FindPeaks ignored anything below 0.002? I only looked at the one data set so I did not see that the other sets have more noise in them. I agree with you that the filtering helps (especially the noisy sets). I was able to get away with leaving the noise on the eta2412 data so I saved the step. I applied your WeinerFilter and all data sets are improved by less noise and the PeakDetection works better.The key to the Peak detect is to not have peaks appear in between the ticks. I raised the threshold so only clusters of peaks form on the "tick". The values I used for both the chiaka and phenix was binarypeaks = PeakDetect[sndData, 0, 0.001, 0.0002]; Where sndData is filtered with your Weiner filter as shown in your post. What did you mean by DC? A resonsant frequency of some kind? DC refers to a steady state offset -- the data with a bunch of 1's and 0's has a mean that is not at zero. There are two ways to deal with that -- one way is to subtract the mean before taking the Power Spectrum. Second is to ignore the peak at the first data point (which is the zero frequency point or "DC") You mentioned the length of the data was low...the other files are around 60 seconds or so. If I use a longer file length, would the 0 padding be required? To be honest, I'm not sure I quite understand the additional padding at all. But you don't need to explain it. The resolution of the Fourier transform or Power spectrum is directly related to the time duration of your signal. If you have a 10 second signal you will have 6 times lower resolution in determining frequency information than if you have a 60 second signal. You can manipulate this a bit by padding the signal with zeros. You will get no more frequency information and your peaks will slowly shrink as you add more zeros, but you can increase the resolution this way. I changed the script to pad with 6 times your signal length to do the 60 second samples and that gave a good resolution. (even with this padding I only am getting 4 decimal places on the period -- is this good enough for you?). There are some clever ways to increase this without increasing the data size but I would need to understand your issues a bit more. (ie. anti-aliasing and downsampling -- you really do not need 44khz for this, taking longer samples instead of padding, and others)I find this problem interesting . I also have some ideas on how to significantly improve things. I'd recommend you email me directly at "NeilS" at the domain in my profile. We should probably talk directly.Regards,Neil
Posted 6 months ago
 I will write you an email today, as soon as I get a chance! Definitely looking forward to working with you on this problem :)!
Posted 6 months ago
 Mor,Even after fixing the algorithm I do not agree with your manual findings on the Phenix. Are you sure it is correct?Also, you have to balance the resolution of the FFT with the resolution of what you collect. For example, you have only 6 ticks per second for 60 seconds = 360 events. This fact limits your resolution as well.Regards,Neil
Posted 6 months ago
 I haven't heard a word why OP didn't say why he hadn't used a Tektronix Oscilliscope and whether the OS had given him the answer rather immediately (as is my guess it would do).OP has to name a range (a time range / sweep / horiz) certainly, volts (vertical if applicable), and also the sampling rate for an answer to be formed, i think.The question of if OS software is freeware and Mathematica has Tektronix software built in, i'm unsure of. Mathematica has many signal analysis functions built-in: but it is not arrange so to make it a drop in replacement for an OS front end.
Posted 6 months ago
Posted 6 months ago
Posted 12 days ago
 I feel like you didn't read i've read it again i see no mention that you already had a solution and were duplicating it for hobbythe OS i showed is obviously just one model - i did not say it was the least cost model needed
Posted 14 days ago
 @ Mor Bo (01/04/2019) : You provided "TickingClock.nb". Can you provide your "eta2412 data1" as another attachment available in this Mathematica Community post? Your internet link to this data file (and other data files above) are no longer working. Thank you! Gilmar Rodriguez-Pierluissi
Posted 14 days ago
 @ Mor Bo (01/04/2019) I see that Dr. Neil Singer was the author of "TickingClock.nb" but, my request about the eta2412 file still holds. Thank you! Gilmar Rodriguez-Pierluissi
Community posts can be styled and formatted using the Markdown syntax.