6

I need an algorithm to detect frequency and phase of a pure sine signal. The frequency of the input signal changes between 0 and 100 Hz.

The value of the signal gets captured with a frequency of 20 kHz (so I get 20000 values per second) - this is given and cannot be changed. I need to detect frequency and phase of this input signal and using PWM generate MCU interrupts with the same frequency as the input signal.

Can anyone suggest what algorithm to use to do this simple and efficient? Maybe Goertzel algorithm?

lennon310
  • 3,590
  • 19
  • 24
  • 27
jurij
  • 325
  • 3
  • 4
  • 8

4 Answers4

7

Note: I originally posted this answer for the Stack Overflow copy of this question, before realizing that it had also been asked here. It somewhat duplicates pichenettes' answer, but I felt it still worth (re)posting here, since it includes some extra details. (Whether those details are useful or not, I'll leave for you and the OP to judge.)

If you know your signal is a pure sine wave, you can just use zero crossing detection. Each cycle of the sine wave will have two zero crossings: one from negative to positive, and one from positive to negative. This will be a lot simpler and more efficient than trying to do something fancy like Fourier transforms.

There are a few details to keep in mind, though:

  • It's OK for the signal to be slightly biased, but if the bias might exceed the amplitude of the signal, you'll need to correct it somehow.

    You can either do this before sampling with an analog high-pass filter, or you can track a moving average of your sampled signal and use it as the "zero level" to compare to. Or you can, instead, look at zero crossings of the difference between successive samples (which correspond to the maximum and minimum of each cycle in the signal) to avoid any bias issues.

  • If your input (or ADC) is noisy, your samples might randomly fluctuate around the zero level when the signal is close to it. In such cases, naïvely comparing successive samples might detect multiple zero crossings where there's really just one.

    One way to fix this issue is to smooth your signal before processing it, but it may be easier and more efficient to implement hysteresis in your zero-crossing detector. That is, you'd only detect a positive-to-negative crossing when the signal dips below some pre-set threshold level −ε, and a negative-to-positive crossing only when it rises above +ε, like this pseudo-C code:

    boolean isPositive = (firstSample > 0);
    while (running) {
        int signalLevel = getNextSample();
        if (isPositive && signalLevel < -threshold) {
            isPositive = false;
            handleZeroCrossing();
        }
        else if (!isPositive && signalLevel > +threshold) {
            isPositive = true;
            handleZeroCrossing();
        }
    }
    

In fact, you may not even need an algorithm for this at all, since, as noted in the comments below, you can implement zero-crossing detection with hysteresis in hardware with a simple Schmitt trigger, which basically converts the sine wave input into a square wave signal with the same frequency and (almost) the same phase, which you can then read as a simple digital input. You might even be able to use the output of the Schmitt trigger to drive the MCU interrupt pin directly.

Ilmari Karonen
  • 245
  • 1
  • 4
  • 4
    If you feed your signal into a high gain op-amp's positive input, and the negative input is at zero volts, the output will be a square wave that goes to the positive rail when the input crosses above the zero volt point. It goes negative when input goes below zero. Count the time between positive-going edges and you have the frequency. If there's a DC bias, the duty cycle of the square wave will be wrong but the overall period will remain good. Having a sharp square wave improves phase measurement circuits, too. – Alan Campbell Nov 08 '14 at 23:34
  • 2
    @AlanCampbell: True. For a noisy signal, a Schmitt trigger would work even better. I was sort of assuming that the OP wanted to do this in software, since they asked specifically for an algorithm, but you're correct that using a DAC for this is really overkill, when all you really need is a hardware comparator and a digital input pin. (For all I know, the OP might even be able to drive the interrupts directly from the Schmitt trigger output.) – Ilmari Karonen Nov 09 '14 at 00:19
  • +1 on hardware op-amp / schmitt trigger + digital acquisition, rather than going through an ADC. – pichenettes Nov 09 '14 at 01:15
  • 1
    I need to do this in software, as I don't have any access to (already determined) hardware. – jurij Nov 10 '14 at 10:09
  • Also, to add, it has been suggested to use Phase Locked Loop algorithm here, but I don't understand neither how it works, nor how to use it. – jurij Nov 10 '14 at 10:12
  • I'd recomment to measure zero cross for both half wave and full period. Then use half mean between positive/negative parts to adjust phase value. – Vitaly Aug 05 '18 at 06:28
3

If the signal is pure sinusoidal and noise-free, you can simply count the number of samples between positive zero crossings - and use this as an estimate of the period of your signal, from which you deduce the frequency. If necessary, wait and average over several periods to get a more robust estimate. There's a good chance this can be achieved without having to sample your signal - directly through an input pin of your microcontroller and through one of its timers.

As for the phase, well, phase relative to what?

pichenettes
  • 19,413
  • 1
  • 50
  • 69
2

Algorithm removed by consensus of replies.

  • On a synthetic noiseless pure signal, I can get 10^(-16) accuracy on a handful of samples. https://www.dsprelated.com/showarticle/1284.php The limit is the precision of the variable, it can do much better as it is theoretically exact. – Cedron Dawg Jul 08 '20 at 11:19
  • @ Cedron Dawg Thanks for the link to your papers. I actually did also estimate the frequency by fitting FFT bins with a parabola (it looks like that is what you are suggesting) which I learned to do a decade ago from this CERN article:https://cds.cern.ch/record/738182/files/ab-2004-023.pdf. The accuracy was similar to the wave-shape fit I described above, but need to study your idea more to know which is better. By the FFT-CERN bin interpolation method I got 95996.90696 Hz sample (unknown error amount) rate instead of 95996.9088 and 95996.9090 Hz by sine shape fit. – fischertranscripts Jul 08 '20 at 11:47
  • The parabolic peak fit is also known as Jacobsen's estimator. It is an approximation and inherently limited requiring a large number of samples to get accurate. 1e-8 is not unreasonable for such a large number of samples. My methods (plural) are much more accurate. For an independent comparison using my original formulas see http://www.tsdconseil.fr/log/scriptscilab/festim/index-en.html The newer formulas are better yet. I also have time domain formulas which are significantly better than zero crossing interpolations starting with https://www.dsprelated.com/showarticle/1051.php – Cedron Dawg Jul 08 '20 at 12:07
  • While OP hasn't specified, your answer deals with a recorded data, but that might be a live input. That this might be the case is due to the requirement of "simple and efficient" which might sound ambiguous, but the rest on this page understood the meaning behind it. In this case, the other answers are very much on point. I the data would be available beforehand, in full, then your answer is not "simple", but it is "efficient". – a concerned citizen Jul 08 '20 at 13:24
  • @aconcernedcitizen The OP is long gone. I was working with the premise of fischertranscripts that it need not be simple. For real time application of a single pure tone, for a quick an dirty estimate, I don't think you can do better than the time domain formulas I referenced. They are quite simple. I use the term "near instantaneous" in the titles very purposefully. My two bin solution only requires the calculation of two bins, not a full DFT, so in that sense it is aslo efficient, but you have to calculate the right two for the best results. You are correct, it is far from simple. – Cedron Dawg Jul 08 '20 at 13:36
  • @CedronDawg I didn't use @ because I thought it would reply to fischertranscripts. Only now I realize it was beneath your reply. An unfortunate slip. :-) It's my impression that OP hoped for a simple if(a,b,c) formula, or similar, when even a "simple comparator" is not that simple. – a concerned citizen Jul 08 '20 at 19:17
  • @aconcernedcitizen Sorry, I did think you were replying to my comment. fischertranscripts approach could be implemented real time as well. Any processing is going to require some sort of lag. As for the OP's question, the lack of desired accuracy makes the question harder to answer. I do think the other answers do a good job. I would advise if your are counting zero crossings that you only pick the ones in the same direction as that will avoid slight DC offset problems. Picking three points, one at a peak, and two flanking near the crossings then applying "Clay's formula" is a one way. – Cedron Dawg Jul 08 '20 at 19:29
  • @CedronDawg "the lack of desired accuracy makes the question harder to answer" -- agreed, it's the main cause for the bandwidth of this whole thread. – a concerned citizen Jul 08 '20 at 19:31
  • "simple and efficient" which might sound ambiguous, but the rest on this page understood the meaning behind it. In this case, the other answers are very much on point.< You have stated that I fail to understand you, or the unstated things you and your pals are implying. I think what is learned from such a reply is: It is not worth posting information here, or reading the non-logic that prevails at this website.

    – fischertranscripts Jul 12 '20 at 08:51
  • @ To the two angry poster above. I have remove my algorithm that bothered you. One of your error reports about me read: >"simple and efficient" which might sound ambiguous, but the rest on this page understood the meaning behind it. In this case, the other answers are very much on point.< Thus you have stated that I fail to understand you, or the unstated things you and your pals are implying, and feel content leaving us that way. I think what is learned from such a reply is: It is not worth posting information here, or reading the non-logic that prevails at this website. – fischertranscripts Jul 12 '20 at 08:57
  • "too long by 472 characters" <<<-----that message sent by https://dsp.stackexchange.com/ censors here who wrote that. Doesn't seem correct, since 1000 characters is a page. So I will try instead to post it in pieces in the following posts (they also suddenly wrote me "no edits allowed", seems like non-sense): – fischertranscripts Jul 22 '20 at 17:52
  • The "original poster" asked about the "Goertzel algorithm", but if you search this page you will find no other occurrence of that word (up till now). – fischertranscripts Jul 22 '20 at 18:01
  • @CedronDawg I find reading such advertising boring, and I call ignoring the original poster's question about "Goertzel", "being off topic", especially since he only had those two short sentences with question marks that he wanted answered (only that first sentence had a problem of being vague). – fischertranscripts Jul 22 '20 at 18:04
  • For example the following statement of your's from above is what I meant by "advertising": "I don't think you can do better than the time domain formulas I referenced. They are quite simple. I use the term 'near instantaneous' in the titles very purposefully." – fischertranscripts Jul 22 '20 at 18:11
  • ??????? I haven't looked at this since making comments, so I am clueless about any thing going on here. I don't make any money off my blog so calling a reference "advertising" is a little off. As I recall, your answer, requiring an atan at every point is not going to be efficient, but as I mentioned to concernedcitizen, it should work. FYI, if you concentrate at zero crossing, your best fit linear on the phase should be just as accurate as when you do the full set of interior points. My comment was addressed to concernedcitizen. – Cedron Dawg Jul 22 '20 at 18:15
  • As to the question asked by lennon310 or jurij: I would not dare post any more info on it here myself, of course; I have learned. I like the question though, and it is so important for (seismology, detection of nuclear tests, Earth core and mantel structure investigation, etc.) to be able to measure the frequencies of sound below 100Hz that the first FFT algorithm was invented to solve it. The more noise the harder it is, if no noise then it is just a trivial curve fit, if tons of noise very hard (e.g. the IERS has to use "wavelet theory" for LOD analysis). – fischertranscripts Jul 22 '20 at 18:35
  • @CedronDawg a sub-100Hz sine wave sampled (say at 8bit resolution) at 20kHz as, specified in his request above, is not linear around the zero crossing as you gallantly profess. And when he reaches 1Hz there is going to be a huge staircase from the digital circuits.The tip of those stairs will shift due to the interacting phases of the digitalizers.... continued – fischertranscripts Jul 22 '20 at 21:53
  • @CedronDawg continued: ... The tiniest rumble of the the DC offset will yank the zero crossing all around the low count bit resolution of many of the 20k per second sample points. Looking at this situation reminds a scientist that there are several laws in Physics that state that the more information you through away, the less accurate you can declare for a measurement, e.g. telescopes with tiny obstructed apertures can not see details. .... continued – fischertranscripts Jul 22 '20 at 21:55
  • @CedronDawg continued(finale): ....That is why your ideas; of throwing away all FFT bins except two, and throwing away all captured data except 4 points; must be harmful for accuracy of measure. – fischertranscripts Jul 22 '20 at 22:00
  • Please read my commentary here, talking about Hz is meaningless without the other parameters specified. https://dsp.stackexchange.com/questions/69186/dft-exercise-in-the-book-understanding-digital-signal-processing-3-ed/69233#69233 The four values in two DFT bins represent an average of all signal values across the frame. The bin indexes represent the "unwinding rate" of the sinusoid in complex space. Thus, doing my calculations is the equivalent of having a regular function that is level, except for noise, and calculating its average height using all the samples. – Cedron Dawg Jul 22 '20 at 22:10
  • (above I meant "The tips of those stairs will shift..." not "tip" (singular) that I wrote by mistake.) – fischertranscripts Jul 22 '20 at 22:13
  • In terms of your approach. you will get the most accurate readings near the zero crossing, so if you took a level function that was noisier between integer domain points and wanted its average, you are better just using the higher SNR values near the integers and discarding the lower SNR values in between. If your frame spans over many integers, you are much better off finding the over all level using your higher quality readings. – Cedron Dawg Jul 22 '20 at 22:13
  • Alternatively, do think about a set of stairs. I can measure the slope by simply measuring the top step and bottom step. If that is a noisy measure, the next step up and down, included in linear regression will get be a better read. The values in the are less accurate and center aren't going to contribute, so why include them? – Cedron Dawg Jul 22 '20 at 22:18
  • @CedronDawg If this is important to you we should compare our self-invented techniques on a common sound file. I have 2015 recording of six tuning forks (8Bit, 3seconds for each fork, 11029.55 samples per second) that I measured with two of my (very different methods) (on was the algorithm I posted here). I got (all in Hz): 261.563,329.663,329.639,523.695,522.954 in 2015. My new way I get:261.5627,329.6635,329.6395,523.6967,522,9492. Can you match that consistency? Lets see. – fischertranscripts Jul 22 '20 at 22:27
  • Can you get 5 digit accuracy on 8 points of data? I am sure I can beat your results easily. Real world signals tend to waver more than the accuracy of my techniques. With synthesized pure tones I get $10^{-16}$ precision as my equations are exact. It's not important I prove it even more. I recommend you read my blog articles if you are curious how it works. https://www.dsprelated.com/blogs-1/nf/Cedron_Dawg.php – Cedron Dawg Jul 22 '20 at 22:35
  • (obviously that last note signals "End of Discussion") Bye to all the readers. – fischertranscripts Jul 23 '20 at 00:09
0

Frequency estimation of a sinusoid (in software) can be done robustly using an autocorrelation.

The center peak (lag of 0) will always be a (or the) maximum value. For a periodic signal, you'll get additional peaks every n lags. This second peak location tells you the period of the signal, in samples. However, you get no phase information.

You can then construct your own sinusoid of the same frequency, with a known phase reference. Then compute the cross-correlation function between your original signal and the reference signal. Now, the first peak value location will tell you the phase difference between the signals. Where the phase = 2*pi*(lag in samples)/(sinusoid period in samples)

wrjohns
  • 101
  • I should add if the sinusoid period isn't a nice integral multiple of the sample period, then you can use the k-th peak in the autocorrelation function, and divide the resulting lag value by k to get a better estimate. – wrjohns Mar 23 '16 at 18:08