-1

I have a periodic signal on a given grid, say :

t = 1:30;
omega = 2*pi/18.431;
phi = -pi+2*pi*rand(1);   % a random phase [-pi,pi]
x = sin(omega*t+phi);     % the signal
x = x+0.5*rand(1,length(x)); % add some noise

enter image description here

Now I want to retrieve the phase phi. There are a few approaches to that, for example, to fit this to a sin, but that takes way too long if I need to do it 1e6 times (unless there's a way to parallelize it?). Another is to use fft. The problem is that my grid is not good enough to pick exactly that frequency, hence the phase associated with it (and I cant change that). How can I manage to get that phase by other means? (and how can I estimate the error in that phase retrieval? I expect traces too noisy to be well estimated and I'd like to know the phase error in that case)

Max
  • 85
  • 6
  • Is this assuming you know the frequency in advance? –  Dec 28 '14 at 22:16
  • yes, I know the frequency in advance, but it is not commensurate with the grid, and undersampled in a sense... – Max Dec 28 '14 at 22:17

2 Answers2

3

The FT is just a correlation of your signal with a cosine and sine wave, and the results give you phase and magnitude. Instead of doing a FT of the whole signal, just calculate the coefficient at your frequency of interest.

As for the accuracy - the longer your signal, the more 'distinct' your inspection frequency will be from adjacent frequencies, so I'm not sure you can get much more accuracy with this method without a longer sample length.

EDIT:

Haven't got Matlab to hand, so going to do this in Python (I've tried to make it look similar to Matlab):

from numpy import arange, pi, cos, dot, exp, angle
from numpy.random import rand

N = 1  # Number of cycles
t = arange(30 * N) + 1  # [1, 2 ... 30 * N]
omega = 2 * pi / 18.431
phi = -pi + 2 * pi * rand(1)[0]
x = cos(omega * t + phi)
x = x + 0.5 * rand(len(x))

coeff = dot(x, exp(-1j * omega * t))
phase = angle(coeff)

print "Actual phase:", phi
print "Calculated phase:", phase

This gives fairly good results even for a single cycle. Occasionally the calculated phase is 2*pi out, which doesn't mean it's giving the wrong answer, but is something to bear in mind for subsequent calculations.

Increasing the number of cycles N produces consistently good results, but this will obviously increase the calculation time so you might want to play around with it and find a compromise for speed vs accuracy (depending on your application).

Hope this helps!

  • can you give me an example of how to do that? – Max Dec 28 '14 at 22:31
  • thanks for the edit, I guess I can learn how to calculate that Fourier coefficient, but am not sure that this is a viable technique (for 1e6 traces). Can you please show me how can I calculate that integral in matlab in an efficient manner? (also `angle` in matlab does just what the `phase` you wrote does) – Max Dec 28 '14 at 22:46
  • 2
    Instead of `atan`, I'd use `atan2` so that you get the full zero to 2*pi phase angle. So, you'd do: `phi = atan2(imag(coeff),real(coeff));` – chipaudette Dec 28 '14 at 23:40
  • I haven't used Matlab in a while, so probably best to treat the example as pseudo-code :-) –  Dec 28 '14 at 23:47
  • 1
    thanks for the edit, I've tried it but didn't get an agreement, shouldn't `fs=1/(dt*N)` with `dt=t(2)-t(1)` and N the # of points in the grid?, also `[0:length(t)]` should be `[1:length(t)]`? (because `x` has length of `t` and in matlab you start with `1` as the first index? (maybe its [0:length(t)-1] ?) can you please add an example that shows that this indeed is the answer? – Max Dec 29 '14 at 07:05
0

You can try searching local minimums/maximum to get both phase and frequency.

EDIT

If you know frequency you can normalise signal, multiply and integrate once with cos and then with sin. This get you directly A and B if signal is presented as 

A*sin(x) +B*cos(x)
Luka Rahne
  • 10,336
  • 3
  • 34
  • 56