15
$\begingroup$

If you have a function $f(t)=A \cdot \sin(\omega t+\phi)$, and reference sin wave $\sin(\omega x)$ what would be a fast algorithm to compute $\phi$?

I was looking at Goertzel algorithm, but it doesn't seem to deal with phase?

$\endgroup$
5
$\begingroup$

Use a DFT at the specific frequency. Then compute amplitude and phase from the real/imag parts. It gives you the phase referenced to the start of the sampling time.

In a 'normal' FFT (or a DFT computed for all N harmonics), you typically compute frequency with f = k*(sample_rate)/N, where k is an integer. Although it may seem sacrilegious (especially to members of the Church of the Wholly Integer), you can actually use non-integer values of k when doing a single DFT.

For instance, suppose you've generated (or obtained) N = 256 points of a sine wave of 27 Hz. (let's say, sample_rate = 200). Your 'normal' frequencies for a 256 point FFT (or N point DFT) would correspond to: f = k*(sample_rate)/N = k*(200)/256, where k is an integer. But a non-integer 'k' of 34.56 would correspond to a frequency of 27 Hz., using the parameters listed above. It's like creating a DFT 'bin' that is exactly centered at the frequency of interest (27 Hz.). Some C++ code (DevC++ compiler) might look as follows:

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;

// arguments in main needed for Dev-C++ I/O
int main (int nNumberofArgs, char* pszArgs[ ] ) {
const long N = 256 ;
double sample_rate = 200., amp, phase, t, C, S, twopi = 6.2831853071795865; 
double  r[N] = {0.}, i[N] = {0.}, R = 0., I = 0. ;
long n ;

// k need not be integer
double k = 34.56;

// generate real points
for (n = 0; n < N; n++) {
    t =  n/sample_rate;
    r[n] = 10.*cos(twopi*27.*t - twopi/4.);
}  // end for

// compute one DFT
for (n = 0; n < N; n++) {
    C = cos(twopi*n*k/N); S = sin(twopi*n*k/N);
    R = R + r[n]*C + i[n]*S;
    I = I + i[n]*C - r[n]*S;
} // end for

cout<<"\n\ndft results for N = " << N << "\n";
cout<<"\nindex k     real          imaginary       amplitude         phase\n";

amp = 2*sqrt( (R/N)*(R/N) + (I/N)*(I/N) ) ;
phase = atan2( I, R ) ;
// printed R and I are scaled
printf("%4.2f\t%11.8f\t%11.8f\t%11.8f\t%11.8f\n",k,R/N,I/N,amp,phase);

cout << "\n\n";
system ("PAUSE");
return 0;
} // end main

//**** end program

(PS: I hope the above above translates well to stackoverflow – some of it might wrap around)

The result of the above is a phase of -twopi/4, as shown in the generated real points (and amp is doubled to reflect the pos/neg frequency).

A few things to note – I use cosine to generate the test waveform and interpret results – you have to be careful about that – phase is referenced to time = 0, which is when you started sampling (ie: when you collected r[0]), and cosine is the correct interpretation).

The above code is neither elegant nor efficient (eg: use a look-up tables for the sin/cos values, etc.).

Your results will get more accurate as you use larger N, and there's a little bit of error due to the fact that the sample rate and N above are not multiples of each other.

Of course, if you want to change your sample rate, N, or f, you'd have to change the code and the value of k. You can plunk down a DFT bin anywhere on the continuous frequency line – just make sure that you're using a value of k that corresponds to the frequency of interest.

$\endgroup$
  • $\begingroup$ This approach can be improved by adjusting N to make k closer to a whole. I posted a separate answer that improres the accuracy of this algorithm. $\endgroup$ – mojuba Nov 15 at 2:48
10
$\begingroup$

The problem can be formulated as (nonlinear) least-squares problem:

$$F(\phi) = \frac{1}{2}\sum^{n}_{i=1}\left[ A \cdot \sin(\omega i+\phi) - f_{i}(\omega) \right]^{2}$$

where $F(\phi)$ is the objective function to minimize with respect to $\phi$.

The derivative is very simple:

$$F'(\phi)=\sum^{n}_{i=1} A\cdot \cos(\omega i+\phi) \left[A \cdot \sin(\omega i+\phi) - f_{i}(\omega)\right]$$

The above objective function can be minimized iteratively using Gradient descent method (first order approximation), Newton method, Gauss-Newton method or Levenberg-Marquardt method (second order approximation - $F''(\phi)$ need to be provided in these).

Obviously, the above objective function has multiple minima because of periodicity, hence some penalty term can be added to discriminate other minima (for example, adding $\phi^{2}$ to the model equation). But I think the optimization will just converge to the nearest minima and you can update the result subtracting $2\pi k, k\in N $.

$\endgroup$
  • $\begingroup$ I dont think you need to penalize because of periodicity no? You can just take whatever minima in phase space it converges to and do a modulu $2\pi$, no? $\endgroup$ – Spacey Sep 6 '12 at 22:48
  • $\begingroup$ @Mohammad Yes, but some optimization techniques may use multiple starting points which should converge to same value or assume a convex function with single global minimizer which can be approximated well with a quadratic. The other benefit is that we end with same result for any starting point $\phi_{0}$. $\endgroup$ – Libor Sep 6 '12 at 23:30
  • $\begingroup$ Interesting. Might I invite you to also take a crack at this related question? :-) $\endgroup$ – Spacey Sep 7 '12 at 14:16
  • $\begingroup$ @Mohammad OK, I have contributed a little there :) $\endgroup$ – Libor Sep 7 '12 at 15:30
  • $\begingroup$ Where does the function fi(w) go? fi(w) is not a constant so when you take a derivative of a non constant how does it become zero? $\endgroup$ – SamFisher83 Sep 10 '12 at 16:13
5
$\begingroup$

There are several different formulations of the Goertzel algorithm. The ones that provide 2 state variables (orthogonal or close to), or a complex state variable, as possible outputs often can be used to calculate or estimate phase with reference to some point in the Goertzel window, such as the middle. The ones that provide a single scalar output alone usually cannot.

You will also need to know where your Goertzel window is in relation to your time axis.

If your signal is not exactly integer periodic in your Goertzel window, the phase estimate around a reference point in middle of the window may be more accurate then referencing phase to the beginning or end.

A full FFT is overkill if you know the frequency of your signal. Plus a Goertzel can be tuned to a frequency not periodic in the FFT length, whereas an FFT will need additional interpolation or zero padding for non-periodic-in-window frequencies.

A complex Goertzel is equivalent to 1 bin of a DFT that uses a recurrence for the cosine and sine basis vectors or FFT twiddle factors.

$\endgroup$
  • $\begingroup$ Isn't the phase estimate anywhere within the window of exactly the same accuracy, because you'd just add $\omega k$ to the phase estimate at the beginning of the window to calculate the phase estimate at sample $k$ within the window ($k=0$ being the beginning of the window)? $\endgroup$ – Olli Niemitalo Jan 8 at 8:30
  • $\begingroup$ No, because adding wk results in a different phase at the end of the window than at the beginning for a non-integer-periodic-in-aperture sinusoid. But a 1-bin DFT computes a single circular phase at that same point. Thus the 3 values will all be different. But the center phase is always related to the ratio of odd/even function, no matter what f0. $\endgroup$ – hotpaw2 Jan 8 at 13:19
  • $\begingroup$ Trying, but I don't get that. $\endgroup$ – Olli Niemitalo Jan 8 at 16:11
  • $\begingroup$ Use a cosine (phase of zero at k=0), tweak the frequency slightly (by a tiny irrational number, but without changing the phase at k=0). A DFT reports the phase has changed! Try the same with a cosine exactly centered at k=N/2. No change at k=N/2 for any df. Same for sin or any mix. Centering the phase reference point show less changes in measured phase with changes in f0. e.g. frequency error does not contribute to increased phase measurement errors. $\endgroup$ – hotpaw2 Jan 8 at 16:56
  • 1
    $\begingroup$ Yes, the phase estimate error being less at the center of the window makes sense if the sinusoid and the Goertzel filter are at different frequencies. In that case, the phase estimate say at the end of the window is biased by a constant which is the product of the distance between the center and the end of the window and the difference between the sinusoid and Goertzel filter frequencies. Subtracting this bias gives the same size error as for the center estimate, but it requires knowing the frequency of the sinusoid. $\endgroup$ – Olli Niemitalo Jan 9 at 8:25
4
$\begingroup$

If your signals are noise-free, you can identify zero crossings in both and determine frequency and relative phase.

$\endgroup$
3
$\begingroup$

That depends on what your definition of "fast" is, how accurate you want your estimate, whether you want $\phi$ or the phase relative to your samplings, and how much noise there is on your function and reference sine wave.

One way to do this is to just take the FFT of $f(t)$ and just look at the bin closest to $\omega$. However, this will depend on $\omega$ being close to the bin center frequency.

So:

  • What do you mean by "fast"?
  • How accurate do you need the estimate?
  • Do you want $\phi$ (phase relative to reference) or phase relative to the start of sampling? Does it matter?
  • What is the noise level on each signal?

PS: I'm assuming you meant $f(t)=A\sin(ωt+ϕ)$, rather than $f({\Huge t})=A\sin(ω{\Huge x}+ϕ)$.

$\endgroup$
2
$\begingroup$

Start point:
1) multiply your signal and reference sin wave:
$F(t) $= A⋅sin(ωt+ϕ)⋅sin(ωt) = 0.5⋅A⋅(cos(ϕ) - cos(2⋅ωt+ϕ))
2) find integral on period $T= \pi /\omega$:
$ I(\phi) = \int_0^T F(t)dt\ = 0.5⋅A⋅cos(ϕ) \cdot T$
3) you can calculate $\phi$:
$cos(\phi) = I(t)/(0.5 \cdot A \cdot T)$

Think about:
how to measure A?
how to determine $\phi$ in $0..(2 \cdot \pi)$ interval? (think about " reference cos wave")

For discrete signal change the integral to sum and carefully choose T!

$\endgroup$
0
$\begingroup$

This is an improvement on @Kevin McGee's suggestion to use a single frequency DFT with a fractional bin index. Kevin's algorithm doesn't yield great results: while at half bins and whole bins it's very precise, also close to the wholes and halves it's also pretty good, but otherwise the error can be within 5%, which is probably not acceptable for most tasks.

I suggest to improve Kevin's algorithm by adjusting $N$, i.e. the length of the DFT window so that $k$ gets as close to a whole as possible. This works since unlike FFT, DFT doesn't require $N$ to be a power of 2.

The code below is in Swift, but should be intuitively clear:

let f = 27.0 // frequency of the sinusoid we are going to generate
let S = 200.0 // sampling rate
let Nmax = 512 // max DFT window length
let twopi = 2 * Double.pi

// First, calculate k for Nmax, and then round it
var k = round(f * Double(Nmax) / S)

// The magic part: recalculate N to make k as close to whole as possible
// We also need to recalculate k once again due to rounding of N. This is important.
let N = Int(k * S / f)
k = f * Double(N) / S

// Generate the sinusoid
var r: [Double] = []
for i in 0..<N {
    let t = Double(i) / S
    r.append(sin(twopi * f * t))
}

// Compute single-frequency DFT
var R = 0.0, I = 0.0
let twopikn = twopi * k / Double(N)
for i in 0..<N {
    let x = Double(i) * twopikn
    R += r[i] * cos(x)
    I += r[i] * sin(x)
}
R /= Double(N)
I /= Double(N)

let amp = 2 * sqrt(R * R + I * I)
let phase = atan2(I, R) / twopi

print(String(format: "k = %.2f    R = %.8f    I = %.8f    A = %.8f    φ/2π = %.8f", k, R, I, amp, phase))
$\endgroup$
  • $\begingroup$ The FFT is simply a way to calculate a DFT efficiently. With modern libraries, the power of two restriction is no longer there. If you only need one or two bin values, then it is better to calculate them directly like you did. For a single pure tone (real or complex), only two bin values are needed to calculate the frequency, phase, and amplitude exactly. See dsprelated.com/showarticle/1284.php. The math is quite sophisticated, but there are links to the articles where the derivations are explained. Linear Algebra is a prerequisite for a true understanding. $\endgroup$ – Cedron Dawg Nov 20 at 4:40
0
$\begingroup$

You could also do this (in numpy notation):

np.arctan( (signal*cos).sum() / (signal*sin).sum() ))

where signal is your phase-shifted signal, cos and sin are the reference signals, and you generate an approximation of an integral over a certain time via summing over the two products.

$\endgroup$

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.