Phase angle from FFT using atan2 - weird behaviour. Phase shift offset? Unwrapping?

Question

I'm testing and performing simple FFT's and I'm interested in phase shift. I geneate simple array of 256 samples with sinusoid with 10 cycles.

I perform an FFT of those samples and receiving complex data (2x128). Than I calculate magnitude of those data and FFT looks like expected:

Then I want to calculate phase shift from fft complex output. I'm using atan2. Combined output fft_magnitude (blue) + fft+phase (red) looks like this: This is pretty much what I expect with a "small" problem. I know this is wrapping but if I imagine unwrapping it, the phase shift in the magnitude peak is reading 36 degrees and I think it should be 0 because my input sinusiod was not shifted at all.

If I shift this -36 deg (blue is in-phase, red is shifted, blue is printed only for reference) the sinusiod looks like this:

And than if I perform an FFT of this red data the magnitude + phase output looks like this:

So it is easy to imagine that unwrapped phase will be close to 0 at the magniture peak. So there is 36 deg offset. But what happenes if I prepare data with 20 cycles per 256 samples and 0 phase shift

If I then perform an FFT, this is an output (magnitude + phase): And I can tell you if will cross the peak point at 72 degrees. So there is now 72 degrees offset.

Can anyone give me a hint why is that happening? Is it right that atan2() phase output is frequency dependent with offset of 2pi/cycles (360 deg/cycles) ? How to unwrap it and get correct results (I couldn't find working C library to unwrap).

This is running on ARM Cortex-M7 processor (embedded).

#define phaseShift 0
#define cycles 20

#include <arm_math.h>
#include <arm_const_structs.h>

float32_t phi = phaseShift * PI / 180; //phase shift in radians
float32_t data[256];  //input data for fft
float32_t output_buffer[256]; //output buffer from fft
float32_t phase_data[128];  //will contain atan2 values of output from fft (output values are complex)
float32_t magnitude[128];  //will contain absolute values of output from fft (output values are complex)
float32_t incrFactorRadians = cycles * 2 * PI / 255;

arm_rfft_fast_instance_f32 RealFFT_Instance;

void setup()
{
  Serial.begin(115200);
  delay(2000);
  arm_rfft_fast_init_f32(&RealFFT_Instance, 256);   //initializing fft to be ready for 256 samples

  for (int i = 0; i < 256; i++)  //print sinusoids
  {
    data[i] = arm_sin_f32(incrFactorRadians * i + phi);
    Serial.print(arm_sin_f32(incrFactorRadians * i), 8); Serial.print(","); Serial.print(data[i], 8); Serial.print("\n"); //print reference in-phase sinusoid and shifted sinusoid (data for fft)
  }
  Serial.print("\n\n");
  delay(10000);

  arm_rfft_fast_f32(&RealFFT_Instance, data, output_buffer, 0); //perform fft

  for (int i = 0; i < 128; i++) //calculate absolute values of an fft output (fft output is complex), and phase shift
  {
    magnitude[i] = output_buffer[i * 2] * output_buffer[i * 2] + output_buffer[(i * 2) + 1] * output_buffer[(i * 2) + 1];
    __ASM("VSQRT.F32 %0,%1" : "=t"(magnitude[i]) : "t"(magnitude[i])); //fast square root ARM DSP function

    phase_data[i] = atan2(output_buffer[i * 2], output_buffer[i * 2 +1]) * 180 / PI;
  }

}

void loop() //print magnitude of fft and phase output every 10 seconds 
{
  for (int i = 0; i < 128; i++)
  {
    Serial.print(magnitude[i], 8); Serial.print(","); Serial.print(phase_data[i], 8); Serial.print("\n");
  }
  
  Serial.print("\n\n");
  delay(10000);
}

Answer 1

An bare FFT plus an atan2() only correctly measures the starting phase of an input sinusoid if that sinusoid is exactly integer periodic in the FFT's aperture width.

If the signal is not exactly integer periodic (some other frequency), then you have to recenter the data by doing an FFTshift (rotate the data by N/2) before the FFT. The FFT will then correctly measure the phase at the center of the original data, and away from the circular discontinuity produced by the FFT's finite length rectangular window on non-periodic-in-aperture signals.

If you want the phase at some point in the data other than the center, you can use the estimate of the frequency and phase at the center to recalculate the phase at other positions.

There are other window functions (Blackman-Nutall, et.al.) that might produce a better phase estimate than a rectangular window, but usually not as good an estimate as using an FFTShift.

Answer 2

To break down the excellent answer by hotpaw2. (Their answers are always so loaded with golden nuggets of information that I spend days learning enough to comprehend the brilliance.)

When an engineer says "integer periodic" they mean your samples that you are feeding into the FFT (the aperture) needs to sample in a way the ensures you capture one full wave of the frequency sin wave.

Think of the sin wave starting at zero and cresting at one then falling below zero into the trough at negative one and then coming back up to zero.

This is one "full cycle". Now if your wave has a period of 10 cycles per second and you sample at 100 samples per second you will have 10 samples per wave.

So now you put 13 samples into an FFT and your phase is off. Why? Well the phase is looking for the wave to smoothly continue forever. You just started a zero for the first sample and dropped off as .25 on the 13th sample. Now the phase calculation tries to connect the two ends and has this jump in the wave. This causes the phase to come out wrong.

What you need to do is select a number of samples to feed into your FFT that you know will contain full waves only.

(NOTE) You are only concerned with the phase of one freq at a time.

AND your sample aperture must not start and end at the sin waves same point. IF you start at zero and end at zero the calculation pasting the two ends together in a forever circle will get two zeros at the transition. So you have to stop one sample short of the repeat point.

Code demonstrating this can be found: Scipy FFT - how to get phase angle