[Math] Computing Fast Fourier Transform in R

fourier analysis

This is a math question with programming application. So I am trying to find 3 things given a certain function in the x domain when transformed into the spectral domain.
1) the Amplitude
2) The Frequency
3) The Phase
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 have two quick questions:

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

Second) How do I convert the amplitude which is abs(fty) into the constant in my consine function above, which is 7? abs(fty)[12]=34351.41 , how do I convert this number to 7?. I appreciate your help!

Best Answer

There are several reasons why you don't get what you expect to see:

  • Your definition of the phase is not what you get from the FFT. The phase you get from the FFT is the phase of the first point of the time domain signal, which in your case is $$\frac{-2 \pi\cdot 50}{9}+32=-2.91\text{ rad}$$ The reason why you only get an approximate value is answered in the next point.

  • The frequency you chose does not lie exactly on one of the FFT bins. This is why you get spectral leakage and neither the magnitude nor the phase will be exact estimates of your original signal. Another way to see this is that within the time window there is not an integer number of periods of your signal.

  • The magnitude of the output of the FFT for a sinusoidal input signal is scaled by $N/2$ (where $N$ is the length of the FFT). So your magnitude estimate is $34351\cdot 2/10001=6.87$ which is a reasonable value given the spectral leakage mentioned above.

If you want to generate a sinusoidal signal which lies exactly on an FFT bin you can use the following formula:

$$y(n)=A\cos\left(2\pi n\frac{i}{N}+\phi\right),\quad n=0,1,\ldots,N-1;\quad i=0,1,\ldots N/2$$

where $n$ is the discrete time index, $N$ is the FFT length (assumed even), $i$ is the desired bin index (note that the index starts at $0$), and $A$ and $\phi$ are some arbitrary constants. In this case the FFT will give you exact values for $A$ and $\phi$. Note that for $\phi$ you get of course the principal value $\phi\in [-\pi,\pi)$.

Related Question