4
$\begingroup$

Let's say I have two signals. The first is a cosine wave and the second is a sine wave. Each oscillates at 0.01 Hz. The sample rate is 1 Hz and the length of time series is 1000 seconds. Each has an amplitude of 1.

My understanding is that an FFT of these two signals should recover a spike in the amplitude spectrum of 0.5 at 0.01 Hz (and similarly at the corresponding negative frequency). This all makes sense and I can get this to work as expected.

But the phase of the FFT is a bit perplexing. My expectation is that the cosine wave will give a phase of 0° at 0.01 Hz, and the sine wave will give a phase of 90° at 0.01 Hz. However, the result I get gives a phase of 3.6° for cosine, and 86.4° for sine. (The negative frequencies are complex conjugates).

Why can the FFT recover the precise amplitudes, but can't do so with the phases? Is there some reason for this? Is it just some sort of numerical or indexing issue or is there some deeper reason? Is my Fourier frequency list incorrect and I'm not actually sampling the spectrum at exactly 0.01 Hz?

MATLAB code to replicate is below.

Any help is appreciated.

f = 0.01; %signal frequency
fs = 1; %sample rate
dt = 1./fs; %length of sample
t = (dt:dt:1000)'; %time vector

df = fs/length(t); %frequency spacing
fAxis = (0:df:(fs-df)) - (fs-mod(length(t),2)*df)/2; %frequency axis with negative freqs

b1 = cos(2*pi*f*t); %first time series signal
b2 = sin(2*pi*f*t); %second time series signal

%FFT each signal and scale
B1 = fftshift((fft(b1)./length(fAxis)));
B2 = fftshift((fft(b2)./length(fAxis)));

%Find index of 0.01 frequency
indf = find(abs(fAxis-0.01)<10^-9);

%Magnitudes look okay. Each returns 0.5
abs(B1(indf)) 
abs(B2(indf))

%But phases???
angle(B1(indf))*180/pi
angle(B2(indf))*180/pi
$\endgroup$
1
  • $\begingroup$ 3.6° = 360°/100 $\endgroup$ Sep 29, 2023 at 1:31

1 Answer 1

5
$\begingroup$

Your time vector is slightly off, it should start at t=0, and end at t=999 (this latter modification keeps the length of the FFT at 1000 to keep a resolution $d_f = 1/1000$, allowing for the frequency component at $f = 0.01\texttt{Hz}$ to fall in one bin):

t = (0:dt:1000-1/fs)'; %time vector

Result:

ans =

   8.8369e-14

ans =

  -90.0000
$\endgroup$
6
  • 6
    $\begingroup$ Again, MATLAB's stupid hard-wired 1-based indexing strikes again. 30 years ago I was trying to get them to make the index base a parameter of an array that can be changed to any integer. So even negative indices, as well as 0, could be accommodated. Even talked with Cleve Moler about it once (back in the '90s). Advocacy fell on deaf ears. $\endgroup$ Sep 29, 2023 at 1:03
  • 9
    $\begingroup$ they should have compromised on a base of 0.5 $\endgroup$
    – benxyzzy
    Sep 29, 2023 at 6:28
  • 4
    $\begingroup$ Although in this case I’m not sure the indexing is the problem, just the way the time vector was constructed. Could’ve been constructed badly with 0 indexing as well! $\endgroup$
    – Jdip
    Sep 29, 2023 at 7:03
  • 2
    $\begingroup$ @robertbristow-johnson older PERL has this feature where you can set the indexing base for arrays by assigning to $[, i.e. $[ = 1; makes your PERL behave like matlab. But: that's global, so it of course breaks kind of everything, and I don't think it has any practical use beyond well showing off that PERL could in theory do that. $\endgroup$ Sep 29, 2023 at 9:40
  • 3
    $\begingroup$ @robertbristow-johnson I'll happily join any tirade against Matlab and 1-based indexing, but for this problem the indexing convention is really besides the point. In Python you also get fencepost errors etc.. The real culprit is the tradition of using raw arrays at all. The proper fix is to have abstractions representing "function sampled on an interval", and FFT methods implemented in a suitable way without having to worry about where the exact sample points are. ODL goes in that direction, though IMO such abstractions work best in static languages. $\endgroup$ Sep 29, 2023 at 13:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.