Is there an algorithm to compute the phase for a single frequecy?

Question

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?

Answer 1

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.

Answer 2

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

Answer 3

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.

Answer 4

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

Answer 5

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}+ϕ)$.

Answer 6

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!

Answer 7

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.

Answer 8

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))

Answer 9

Phase, frequency, and amplitude can all be computed exactly using the inverse solution to the sine DFT closed form solution. Here's the worst case, $N=3$, noiseless:

Here's noisy:

It's accomplished using A Two Bin Solution, by Cedron Dawg. It's basically just $O(1)$ on top of FFT, so extremely fast. On frequency, see a competitive comparison, where it possibly achieves state of the art in compute-limited settings.

Reproducing code and fuller discussion at Extract Sine Phase and Amplitude - accurate and robust method.