[Math] Computing the phase in a Fourier Transform

fourier analysisharmonic-analysis

This is a math question with programming application. So I am trying to find the phase from the following signal $y$.

$$y=A\cos\left(2\pi \omega +\phi\right)$$

I am specifically trying to estimate $\phi$ from the signal.

In R (statistical software) I have coded the following function:

y=7*cos(2*pi*(seq(-50,50,by=.01)*(1/9))+32)
fty=fft(y,inverse=F)
angle=atan2(Im(fty), Re(fty))
x=which(abs(fty)[1:(length(fty)/2)]==max(abs(fty)[1:(length(fty)/2)]))
par(mfcol=c(2,1))
plot(seq(-50,50,by=.01),y,type="l",ylab = "Cosine Function")
plot(abs(fty),xlim=c(x-30,x+30),type="l",ylab="Spectral Density in hz")

I know I can compute the frequency by taking the bin value and diving it by the size of the interval(total time of the domain). Since the bins started at 1, when it should be zero, it would thus be
frequency=(BinValue-1)/MaxTime, which does get me the 1/9'th I have in the function above.

I am having trouble computing the phase. the density function peaks at 12 (see bottom graph), shouldn't the value then the value of the phase be 2*pi+angle[12]
but I am getting a value of

angle[12]
[1] -2.558724

which puts the phase at 2*pi+angle[12]=3.724462. But that's wrong the phase should be 32 radians. What am I doing wrong?

Best Answer

Two (maybe three) things are going on.

First, atan2() will only ever return angles in the interval $(-\pi, \pi]$. It can't possibly know how many times you've wrapped around by $2\pi$ radians. You should notice that

$$-2.558724 \space \mathrm{rad} + 2\pi(5) \space \mathrm{rad} + \pi \space \mathrm{rad}= 32 \space \mathrm{rad}$$

which shows $5$ full rotations that atan2() couldn't possibly know about. But you also seem to have an error of $\pi$ radians, so where is that coming from?

Seconds, let's look at the Fourier Transform of

$$y = A \cos(2\pi f_ct+\phi_0)$$ $$\begin{align*}\mathscr{F}\{y\} &= \int_{-\infty}^{\infty} A \cos(2\pi f_ct+\phi_0) e^{-2\pi i s t} \mathrm{d}t \\ &= \dfrac{A}{2}\left[\int_{-\infty}^{\infty} e^{i(2\pi f_ct+\phi_0)} e^{-2\pi i s t} \mathrm{d}t + \int_{-\infty}^{\infty} e^{-i(2\pi f_ct+\phi_0)} e^{-2\pi i s t} \mathrm{d}t\right]\\ &= \dfrac{A}{2}\left[e^{i\phi_0}\int_{-\infty}^{\infty} e^{-2\pi i (s-f_c) t} \mathrm{d}t + e^{-i\phi_0}\int_{-\infty}^{\infty} e^{-2\pi i (s+f_c) t} \mathrm{d}t\right]\\ &= \dfrac{A}{2}\left[e^{i\phi_0}\delta(s-f_c) + e^{-i\phi_0}\delta(s+f_c)\right]\\ \end{align*}$$

Note how the phase term of the delta function sitting at the negative frequency, $\delta(s+f_c)$, has a negative phase, $e^{i(-\phi_0)}$. That could be a problem if you've selected that peak, but I doubt it in this case.

Third, two nuances of the DFT that most people don't immediately realize is that the phase is relative to the first element of the time series data ($t=0$ is the time of the first time series data point), and that the time series data is cyclic as far as the DFT math is concerned (the second half of your time series data can be used for the negative time values with the time index of the last element of the time series being just less than $t=0$).

Your cosine wave is shifted to the right in time, as far as the DFT is concerned. This will result in a modulation by a frequency dependent complex exponential in the frequency domain, and will impart a frequency dependent phase rotation on the peak you are examining.

This property of the DFT is likely giving you the additional phase shift of $\pi$ in your code. Fix your cosine time series data so that the first element corresponds to $t=0$, to get rid of this problem.

Related Question