Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 2)
Introduction
This is an article that is a continuation of a digression from trying to give a better understanding of the Discrete Fourier Transform (DFT). It is recommended that my previous article "Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 1)"[1] be read first as many sections of this article are directly dependent upon it.
A second family of formulas for calculating the frequency of a single pure tone in a short interval in the time domain is presented. It is very similar to the first family presented in the previous article. The formulas also work for complex tones. The formulas are still exact, meaning no approximations were taken in the derivation. The advantage of the second family over the first is that it offers a more robust solution over the same sampling interval. It also requires marginally more calculations to achieve this. Like the first family, this one should be evaluated at peaks of real tones and the formulas become indeterminate and noise sensitive at zero crossings. In the complex case, there is no such limitation.
Pure Tone Signal Definitions
This is the equation used in all my real tone blog articles to describe a discrete single pure tone: $$ S_n = M \cdot \cos( \alpha n + \phi ) \tag {1} $$ This is the corresponding equation for a pure complex tone. $$ S_n = M \cdot e^{i(\alpha n + \phi )} \tag {2} $$ In both types, the parameter being solved for with these families of formulas is $\alpha$, which is the one pertaining to the frequency.
One of the Family
Here is a member of the second family of formulas: $$ \begin{aligned} P_{n,1} &= S_{n+2} + S_{n-2} \\ P_{n,2} &= S_{n+4} + S_{n-4} \\ P_{n,3} &= S_{n+6} + S_{n-6} \\ P_{n,4} &= S_{n+8} + S_{n-8} \end{aligned} \tag {3} $$ $$ \alpha = \frac{1}{2} \cos^{-1}\left( \frac{ 30 S_n + 26 P_{n,1} + 16 P_{n,2} + 6 P_{n,3} + P_{n,4} } { 40 S_n + 30 P_{n,1} + 12 P_{n,2} + 2 P_{n,3} } \right) \tag {4} $$ Like in the previous article, this one is the $d=2,k=4$ case. The difference is that this family uses all the neighbor pairs in both the numerator and denominator. This makes the formulas more robust, that is noise resistant, though it does make for a few more calculations.
See Appendix A for a list of the formulas up to $k=9$.
Neighbor Pairs of Signal Values
The two neighbor sample values, on either side of the sample center by the same interval, are added together to form a neighbor pair sum. For both types of signals the formula is the same. In this article, the sample spacing will be incorporated from the start. $$ P_{n,m} = S_{n+md} + S_{n-md} = 2 S_n \cos[ (\alpha d) m ] \tag {5} $$ $$ S_n \cos[ (\alpha d) m ] = \frac{P_{n,m}}{2} \tag {6} $$
Powers of a Cosine Binomial
Instead of taking the powers of the cosine function, this family is based on taking the powers of a binomial containing the cosine. The calculations are a little more complicated and are shown in Appendix B. $$ \left[1 + \cos(\theta)\right]^0 = 1 \tag {7} $$ $$ \left[1 + \cos(\theta)\right]^1 = 1 + \cos(\theta) \tag {8} $$ $$ \left[1 + \cos(\theta)\right]^2 = \frac{3}{2} + 2 \cos(\theta) + \frac{1}{2} \cos(2\theta) \tag {9} $$ $$ \left[1 + \cos(\theta)\right]^3 = \frac{5}{2} + \frac{15}{4} \cos(\theta) + \frac{3}{2} \cos(2\theta) + \frac{1}{4} \cos(3\theta) \tag {10} $$ $$ \left[1 + \cos(\theta)\right]^4 = \frac{35}{8} + 7 \cos(\theta) + \frac{7}{2} \cos(2\theta) + \cos(3\theta) + \frac{1}{8} \cos(4\theta) \tag {11} $$ The main difference is that all the power terms appear instead of skipping every other one like the first family does.
Signal Calculations
The signal calculation takes the general form: $$ W_{n,k} = S_n \left[1 + \cos[ (\alpha d) m ]\right]^k \tag {12} $$ Where $n$ is the center sample index and $k$ is the degree of the exponentiation.
The first few examples can be found from the cosine power series in the previous section. $$ W_{n,0} = S_n \tag {13} $$ $$ W_{n,1} = S_n + S_n \cos[ (\alpha d) m ] = S_n + \frac{1}{2} P_{n,1} \tag {14} $$ $$ \begin{aligned} W_{n,2} &= \frac{3}{2} S_n + 2 S_n \cos[ (\alpha d) m ] + \frac{1}{2} S_n \cos(2\alpha) \\ &= \frac{3}{2} S_n + P_{n,1} + \frac{1}{4} P_{n,2} \end{aligned} \tag {15} $$ $$ W_{n,3} = \frac{5}{2} S_n + \frac{15}{8} P_{n,1} + \frac{3}{4} P_{n,2} + \frac{1}{8} P_{n,3} \tag {16} $$ $$ W_{n,4} = \frac{35}{8} S_n + \frac{7}{2} P_{n,1} + \frac{7}{4} P_{n,2} + \frac{1}{2} P_{n,3} + \frac{1}{16} P_{n,4} \tag {17} $$ The coefficients for the first nine formulas can be found in Appendix C.
Frequency Calculation
Finding the frequency value is similar to how it was done with the first family. Taking the quotient of two adjacent signal calculations yields an expression for the frequency term. $$ q = \frac{ W_{n,k} }{ W_{n,k-1} } = \frac{S_n \left[1 + \cos[ (\alpha d) m ]\right]^{k}}{S_n \left[1 + \cos[ (\alpha d) m ]\right]^{k-1}} =1 + \cos[ (\alpha d) m ] \tag {18} $$ Solving for the frequency is then very straighforward: $$ \alpha = \cos^{-1}(q-1) = \cos^{-1}\left(\frac{ W_{n,k} }{ W_{n,k-1} }-1\right) \tag {19} $$
Base Case
The base case occurs when $k$ is set to zero and $d$ is set to one. $$ q = \frac{ W_{n,1} }{ W_{n,0} } = \frac{S_n + \frac{1}{2} P_{n,1}}{S_n} = 1 + \frac{ P_{n,1} }{2 S_n} \tag {20} $$ Applying a little algebra shows that the base case of this family is the same as the base case of the first family. $$ \alpha = \cos^{-1}\left( \frac{ P_{n,1} }{2 S_n} \right) = \cos^{-1}\left( \frac{ S_{n+1}+S_{n-1} }{2 S_n} \right) \tag {21} $$
Better Signal Value
Similar as what was done in the previous article, (12) can be rearranged to give a formula for an estimated value of the signal at a given sample point. $$ G_n = \frac{ W_{n,k} }{ \left[1 + \cos(\alpha d)\right]^k } = \frac{ W_{n,k} }{ q^k } \tag {22} $$ In the noiseless case, this guess will match the true signal value. In a noisy case it is likely to be more accurate.
Simplifying the Formula
The fractions in the coefficients can be eliminated by multiply the numerator and denominator of the $q$ calculation by the appropriate power of two. $$ q = \frac{ W_{n,k} }{ W_{n,k-1} } = \frac{ 2^k W_{n,k} }{ 2^k W_{n,k-1} } \tag {23} $$ Subtracting the one from $q$ to get $r$ can also be precalculated, resulting is simpler formulas for the frequency calculation. $$ r = q - 1 = = \frac{ 2^k W_{n,k} }{ 2^k W_{n,k-1} } - 1 = \frac{ 2^k W_{n,k} - 2^k W_{n,k-1} }{ 2^k W_{n,k-1} } \tag {24} $$ The actual calculation for the frequency also is simplified. $$ \alpha = \frac{1}{d} \cos^{-1}(r) \tag {25} $$ The coefficients for the first nine formulas in the family can be found in Appendix A.
Numerical Example
Here is a sample calculation of the formula given in (3). Since it is a noiseless case, the accuracy is good to the precision of the numerical variables, which is higher than the values displayed.
alpha = 0.0626894 S[144] = 2.6701126 S[145] = 2.7086362 S[146] = 2.7365186 S[147] = 2.7536500 S[148] = 2.7599633 S[149] = 2.7554336 S[150] = 2.7400787 S[151] = 2.7139589 S[152] = 2.6771768 Pn1 = 5.5090836 Pn2 = 5.4765972 Pn3 = 5.4225951 Pn4 = 5.3472894 Vn3 = 22.0147123 Vn4 = 43.9861803 q = 1.9980357 alpha_n = 0.0626894 Gn = 2.7599633
The Larger Family
Replacing the one in (12) with a variable $x$ makes the two families of equations covered in this and my previous article merely two points on a much larger set. $$ W_{n,k}(x) = S_n \left[x + \cos[ (\alpha d) m ]\right]^k \tag {26} $$ The first family from the previous article corresponds to $x=0$ and the second family, covered in this article, corresponds to $x=1$. The effect of varying x is changing the weighting of the sample points around the center.
People who are good with filters (and I'm not really one) will recognize that the numerator and denominator of the frequency calculations are equivalent to finite impulse response (FIR) filters. I have steered away from employing this terminology because the derivation of these formulas is inherently center based and FIRs are generally defined relative to the last sample. Therefore calling them FIRs may lead to more confusion than benefit.
Limited Applicability
The formulas in this article and the previous article are really only applicable to single tone low noise signals. This makes their usefulness rather limited except for those specific applications. The mathematics are pretty neat though and in those few cases where these formulas are useful, they are probably going to be the best ones. I have not done any extensive testing to see what settings are optimal for any given situation. One reason is the limited applicability. Another reason is because of the wide range of parameter settings would make the testing large and cumbersome, and the results would correspondingly be large and cumbersome also. Having said that, it is still on my todo list, but readers are welcome to beat me to it.
Conclusion
A second family of exact time domain frequency formulas for single tones has been derived. It turns out just to be one of infinitely many such families. A lot more testing is needed to determine what the optimal parameter values are for any given situation.
Reference
[1] Dawg, Cedron, Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 1)
Appendix A. Simplified Formulas
These are the coefficients for the first nine simplified formulas given by (23). The numerator and denominator need to be dotted with the vector $ \begin{bmatrix} S_n & P_{n,1} & P_{n,2} & P_{n,3} & ... \end{bmatrix} $ in order to be evaluated. $$ 1:\frac{ \begin{bmatrix} 0 & 1 \end{bmatrix} }{ \begin{bmatrix} 2 \end{bmatrix} } \tag {27} $$ $$ 2:\frac{ \begin{bmatrix} 2 & 2 & 1 \end{bmatrix} }{ \begin{bmatrix} 4 & 2 \end{bmatrix} } \tag {28} $$ $$ 3:\frac{ \begin{bmatrix} 8 & 7 & 4 & 1 \end{bmatrix} }{ \begin{bmatrix} 12 & 8 & 2 \end{bmatrix} } \tag {29} $$ $$ 4:\frac{ \begin{bmatrix} 30 & 26 & 16 & 6 & 1 \end{bmatrix} }{ \begin{bmatrix} 40 & 30 & 12 & 2 \end{bmatrix} } \tag {30} $$ $$ 5:\frac{ \begin{bmatrix} 112 & 98 & 64 & 29 & 8 & 1 \end{bmatrix} }{ \begin{bmatrix} 140 & 112 & 56 & 16 & 2 \end{bmatrix} } \tag {31} $$ $$ 6:\frac{ \begin{bmatrix} 420 & 372 & 255 & 130 & 46 & 10 & 1 \end{bmatrix} }{ \begin{bmatrix} 504 & 420 & 240 & 90 & 20 & 2 \end{bmatrix} } \tag {32} $$ $$ 7:\frac{ \begin{bmatrix} 1584 & 1419 & 1012 & 561 & 232 & 67 & 12 & 1 \end{bmatrix} }{ \begin{bmatrix} 1848 & 1584 & 990 & 440 & 132 & 24 & 2 \end{bmatrix} } \tag {33} $$ $$ 8:\frac{ \begin{bmatrix} 6006 & 5434 & 4004 & 2366 & 1092 & 378 & 92 & 14 & 1 \end{bmatrix} }{ \begin{bmatrix} 6864 & 6006 & 4004 & 2002 & 728 & 182 & 28 & 2 \end{bmatrix} } \tag {34} $$ $$ 9:\frac{ \begin{bmatrix} 22880 & 20878 & 15808 & 9828 & 4928 & 1940 & 576 & 121 & 16 & 1 \end{bmatrix} }{ \begin{bmatrix} 25740 & 22880 & 16016 & 8736 & 3640 & 1120 & 240 & 32 & 2 \end{bmatrix} } \tag {35} $$
Appendix B. Cosine Binomial Calculations
These are the raw calculations for the first few cosine binomial expansions. Each equation is then tested for accuracy by plugging in $0$, $\pi$, and $\pi/2$ to ensure the results are $2^k$, $0$, and $1$ respectively. $$ \begin{aligned} \left[1 + \cos(\theta)\right]^2 &= 1 + 2 \cos(\theta) + \cos^2(\theta) \\ &= 1 + 2 \cos(\theta) + \left[ \frac{e^{i\theta}+e^{-i\theta}}{2}\right]^2 \\ &= 1 + 2 \cos(\theta) + \frac{e^{i2\theta}+2+e^{-i2\theta}}{4} \\ &= \frac{3}{2} + 2 \cos(\theta) + \frac{1}{2} \cos(2\theta) \end{aligned} \tag {36} $$ $$ \begin{aligned} \theta = 0 \quad &\Rightarrow \quad \frac{3}{2} + 2 + \frac{1}{2} = 4 \\ \theta = \pi \quad &\Rightarrow \quad \frac{3}{2} - 2 + \frac{1}{2} = 0 \\ \theta = \frac{\pi}{2} \quad &\Rightarrow \quad \frac{3}{2} - \frac{1}{2} = 1 \end{aligned} \tag {37} $$ $$ \begin{aligned} \left[1 + \cos(\theta)\right]^3 &= 1 + 3 \cos(\theta) + 3 \cos^2(\theta) +\cos^3(\theta) \\ &= 1 + 3 \frac{e^{i\theta}+e^{-i\theta}}{2} + 3 \left[ \frac{e^{i\theta}+e^{-i\theta}}{2}\right]^2 + \left[ \frac{e^{i\theta}+e^{-i\theta}}{2}\right]^3 \\ &= 1 + \frac{3}{2} ( e^{i\theta}+e^{-i\theta} ) + \frac{3}{4} ( e^{i2\theta}+2+e^{-i2\theta} ) \\ &+ \frac{1}{8} ( e^{i3\theta}+3e^{i\theta}+3e^{-i\theta}+e^{-i3\theta} ) \\ &= \frac{5}{2} + \frac{15}{8} ( e^{i\theta}+e^{-i\theta} ) + \frac{3}{4} ( e^{i2\theta}+e^{-i2\theta} ) + \frac{1}{8} ( e^{i3\theta}+e^{-i3\theta} ) \\ &= \frac{5}{2} + \frac{15}{4} \cos(\theta) + \frac{3}{2} \cos(2\theta) + \frac{1}{4} \cos(3\theta) \end{aligned} \tag {38} $$ $$ \begin{aligned} \theta = 0 \quad &\Rightarrow \quad \frac{5}{2} + \frac{15}{4} + \frac{3}{2} + \frac{1}{4} = 8 \\ \theta = \pi \quad &\Rightarrow \quad \frac{5}{2} - \frac{15}{4} + \frac{3}{2} - \frac{1}{4} = 0 \\ \theta = \frac{\pi}{2} \quad &\Rightarrow \quad \frac{5}{2} -\frac{3}{2} = 1 \end{aligned} \tag {39} $$ $$ \begin{aligned} \left[1 + \cos(\theta)\right]^4 &= 1 + 4 \cos(\theta) + 6 \cos^2(\theta) + 4 \cos^3(\theta) + \cos^4(\theta) \\ &= 1 + 4 \frac{e^{i\theta}+e^{-i\theta}}{2} + 6 \left[ \frac{e^{i\theta}+e^{-i\theta}}{2}\right]^2 \\ &+ 4 \left[ \frac{e^{i\theta}+e^{-i\theta}}{2}\right]^3 + \left[ \frac{e^{i\theta}+e^{-i\theta}}{2}\right]^4 \\ &= 1 + 2 ( e^{i\theta}+e^{-i\theta} ) + \frac{3}{2} ( e^{i2\theta}+2+e^{-i2\theta} ) \\ &+ \frac{1}{2} ( e^{i3\theta}+3e^{i\theta}+3e^{-i\theta}+e^{-i3\theta} ) \\ &+ \frac{1}{16} ( e^{i4\theta}+4e^{i2\theta}+6+4e^{-i2\theta}+e^{-i4\theta} ) \\ &= \frac{35}{8} + \frac{7}{2} ( e^{i\theta}+e^{-i\theta} ) + \frac{7}{4} ( e^{i2\theta}+e^{-i2\theta} ) \\ &+ \frac{1}{2} ( e^{i3\theta}+e^{-i3\theta} ) + \frac{1}{16} ( e^{i4\theta}+e^{-i4\theta} ) \\ &= \frac{35}{8} + 7 \cos(\theta) + \frac{7}{2} \cos(2\theta) + \cos(3\theta) + \frac{1}{8} \cos(4\theta) \end{aligned} \tag {40} $$ $$ \begin{aligned} \theta = 0 \quad &\Rightarrow \quad \frac{35}{8} + 7 + \frac{7}{2} + 1 + \frac{1}{8} = 16 \\ \theta = \pi \quad &\Rightarrow \quad \frac{35}{8} - 7 + \frac{7}{2} - 1 + \frac{1}{8} = 0 \\ \theta = \frac{\pi}{2} \quad &\Rightarrow \quad \frac{35}{8} - \frac{7}{2} + \frac{1}{8} = 1 \end{aligned} \tag {41} $$ The general form looks like this: $$ \begin{aligned} \left[1 + \cos(\theta)\right]^n &= \sum_{k=0}^{n} { \binom{n}{k} \cos^k(\theta) } \\ &= \sum_{k=0}^{n} { \binom{n}{k} \left[ \frac{ e^{i\theta} + e^{-i\theta} }{ 2 } \right]^k } \\ &= \sum_{k=0}^{n} { \binom{n}{k} \frac{ 1 }{ 2^k } \sum_{j=0}^{k} { \binom{k}{j} e^{i\theta (k-j)} e^{-i\theta j } } } \\ &= \sum_{k=0}^{n} { \sum_{j=0}^{k} { \frac{ 1 }{ 2^k } \binom{n}{k} \binom{k}{j} e^{i\theta (k-2j)} } } \\ &= \sum_{k=0}^{n} { \sum_{j=0}^{k} { \frac{ 1 }{ 2^k } \frac{n!}{ (n-k)! \cdot j! \cdot (k-j)! } e^{i\theta (k-2j)} } } \end{aligned} \tag {42} $$
Appendix C. Calculation Coefficients
These are the coefficients for the first nine signal calculation. When used in the frequency calculation it is better to use the simplified coefficient values found in Appendix A.
1: | 1 | $\frac{1}{2}$ | ||||||||
---|---|---|---|---|---|---|---|---|---|---|
2: | $\frac{3}{2}$ | 1 | $\frac{1}{4}$ | |||||||
3: | $\frac{5}{2}$ | $\frac{15}{8}$ | $\frac{3}{4}$ | $\frac{1}{8}$ | ||||||
4: | $\frac{35}{8}$ | $\frac{7}{2}$ | $\frac{7}{4}$ | $\frac{1}{2}$ | $\frac{1}{16}$ | |||||
5: | $\frac{63}{8}$ | $\frac{105}{16}$ | $\frac{15}{4}$ | $\frac{45}{32}$ | $\frac{5}{16}$ | $\frac{1}{32}$ | ||||
6: | $\frac{231}{16}$ | $\frac{99}{8}$ | $\frac{495}{64}$ | $\frac{55}{16}$ | $\frac{33}{32}$ | $\frac{3}{16}$ | $\frac{1}{64}$ | |||
7: | $\frac{429}{16}$ | $\frac{3003}{128}$ | $\frac{1001}{64}$ | $\frac{1001}{128}$ | $\frac{91}{32}$ | $\frac{91}{128}$ | $\frac{7}{64}$ | $\frac{1}{128}$ | ||
8: | $\frac{6435}{128}$ | $\frac{715}{16}$ | $\frac{1001}{32}$ | $\frac{273}{16}$ | $\frac{455}{64}$ | $\frac{35}{16}$ | $\frac{15}{32}$ | $\frac{1}{16}$ | $\frac{1}{256}$ | |
9: | $\frac{12155}{128}$ | $\frac{21879}{256}$ | $\frac{1989}{32}$ | $\frac{4641}{128}$ | $\frac{1071}{64}$ | $\frac{765}{128}$ | $\frac{51}{32}$ | $\frac{153}{512}$ | $\frac{9}{256}$ | $\frac{1}{512}$ |
- Comments
- Write a Comment Select to add a comment
Yo Dawg!
I have read Rick Lyons and your series of articles on time domain, single frequency estimation with great interest. I'm trying to estimate the fundamental frequency of a relatively harmonically rich signal. The signal is non-stationary but changes frequency relatively slowly. Although there is noise in the signal, the peaks of the signal are easy to discern by eye. Simple peak finding algorithms seem to function well.
Frequency domain, transform based methods seem to be a poor fit. I'm only interested in the fundamental. The block by block operation of transforms also introduces too much latency in frequency estimation.
Are any of the algorithms that you have so rigorously presented be suitable for my situation? If so , then I will give the most appropriate algorithm a try. Of course, I would be happy to share my results with you.
Thanks so much for the considerable effort that you have expended in developing and documenting this work.
Wayne
Thank you for the kind words.
It was indeed Lyons' article on time domain frequency estimation that took me off my path of DFT analysis on a detour of time domain equations. I have one more article coming up in the time domain and then I'm back to the frequency domain. Unfortunately, the premise of this article and the previous one is that the signal is a single pure tone. Your harmonics are going to render these formulas useless.
The latency issue is interesting. The main benefit of using the formulas I have derived is to minimize latency. I am puzzled that you say the frequency domain methods are a poor fit for you. Particularly if you have a rough idea of the waveform size and can pick your DFT frame accordingly.
If your signal is non-stationary then any zero crossing based methods are going to be useless.
Since the amplitude and phase aren't important, you can apply exponential smoothing to your signal to mitigate the effect of noise without any impact on latency. The smoothing will also reduce your higher harmonics relative to your lower ones.
I believe autocorrelation is the best solution for your situation. Since you have a rough idea of the frequency to begin with, you can make an implementation that is relatively low cost computation wise. If you have smoothed the signal, you will even be able to reduce computation load further by spacing the samples you use in your calculation.
You can send me your email address using the little red envelope under "About Cedron Dawg" on the right side of the page. I'm interested in how it works out for you.
There are some other smoothing techniques, similar to exponential smoothing, that sample by intervals instead of adjacent points to knock out harmonics. Somewhat equivalent to comb filter effects.
Ced
Ced
I will be pleased to share more info with you via your email address.
Whether your methods are applicable or not, I greatly enjoyed how you formulated and attacked the problem. I remembered a lot of nearly forgotten math as well (blush).
I tried auto-correlation but unless I do the correlation via an FFT, the N squared time complexity of direct correlation takes too much time to be of much use on the tiny embedded processor that I am forced to use.
In MATLAB, auto-correlation is "da bomb". I even made a MATLAB animation of auto-correlation estimating the fundamental frequency of actual signals captured from the "machinery".
Stay tuned for more info. I will share some plots of the waveform as well.
Thanks again for all your insights and effort. I believe that I owe you more than a few beers (if you're not a teatoaler)!
Wayne
Yes, I am interested to see your results and I'm very glad I could be helpful. No, I don't drink in the real world, but feel free to dish me out some DspRelated beers. I haven't gotten anything through the message box so I'll give you my email here: cedron ta exede tod net.
Auto-correlation need not be NxN, in your case it can by Nx3 since you roughly know the frequency ahead of time. The auto-correlation is going to be roughly a parabola near the optimal point. All you have to do is calculate three of those points and then do a parabola fit. This will greatly reduce your computation load.
Also, when you are calculating the auto-correlation, you can use every other sample, or every third sample, to bring the calculation load down to (N/2)x3 or (N/3)x3 respectively.
Ced
To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.
Please login (on the right) if you already have an account on this platform.
Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers: