interpreting frequency vs phase plot in fft

19 views (last 30 days)
joonhak lee
joonhak lee on 9 Nov 2021
Edited: David Goodmanson on 1 Dec 2021
I have problems understanding the second plot above. Why do the phases look discontinuous at the frequencies with peaks in magnitude? How can I determine the desirable initial phases from the plot?
  4 Comments

Sign in to comment.

Answers (2)

Chunru
Chunru on 10 Nov 2021
Edited: Chunru on 10 Nov 2021
The output of the FFT has zero values in most frequency bins(with some leaking effect due to the fact that frequencies of sinosoidal signal are not falling at the integer bins). The phase for a zero value is meaningless. The strange phase plot also has something to do with the unwrap. If we apply a tolerance threshold to make those small numbers (due to leaking and computational error) to be zeros, the phase plot will look "nicer".
t = 0:1/100:10-1/100; % Time vector
x = sin(2*pi*15*t) + sin(2*pi*40*t); % Signal
x = sin(2*pi*25*t) + sin(2*pi*12.5*t); % integral freq bin
n = 512;
y = fft(x,n);
tol = 1e-6;
y(abs(y) < tol) = 0;
m = abs(y);
%p = unwrap(angle(y));
p = (angle(y));
f = (0:length(y)-1)*100/length(y);
subplot(2,1,1)
plot(f,m)
title('Magnitude')
ax = gca;
ax.XTick = [15 40 60 85];
subplot(2,1,2)
plot(f,p*180/pi)
title('Phase')
ax = gca;
ax.XTick = [15 40 60 85];

David Goodmanson
David Goodmanson on 12 Nov 2021
Edited: David Goodmanson on 12 Nov 2021
Hello JL,
Since
sin(2*pi*f*t) = ( e^(2*pi*i*f*t)-e^(-2*pi*i*f*t) )/(2i), (1)
the sine signal has both positive and negative frequency components. The frequency plots have f=0 as the first point in the array, so the positive frequencies f=15 and f=40 appear in the lower half of the plot, and by the properties of the fft, the negative frequencies f =-15 and f = -40 appear in the upper half as +85 and +60 respectively.
oo The first part of the code you cited with n = 1000 instead of 512, has an exact number of cycles in the time array. In the frequency domain you get spikes at +-15, +-40 and zeroes everywhere else. The amplitudes at the splkes are pure imaginary because of the factor of i in the denominator of (1), so the phase is +-90 degrees there. Pretty straightforward. Everywhere else the amplitude is 0, and Matlab uses 0 for the angle, although the answer should really be NaN. Wrong of Matlab not to use NaN in my opinion, but there you are.
oo The second part of the code cuts the points to 512. Now there are not an exact number of cycles in the window, so the peaks broaden, and there are four 180 degree drops. That means as a function of frequency, a factor of -1 appears when crossing those locations.
First of all, those drops have nothing to do with phase wrap. If you edit out the unwrap function you will se that there is a true 2*pi phase discontinuity at the center of the plot (removed by unwrap), but the same four drops are still there. Taking a positive signal frequency f0 as an example, and approximating the fft with an integral for simplicity, the fft is
g(f) = Int{0,a} e^(2*pi*i*(f0-f)*t) dt
where 'a' is some value appropriate to 512 points. The result is
g(f) = ( e^(2*pi*i*(f0-f)*a) -1 )/(2*pi*i*(f0-f))
= e^(2*pi*i*(f0-f)*a/2)*sin(2*pi*(f0-f)*a/2) / (pi*(f0-f))
The first factor has slowly varying phase. The sine factor is real and does not affect the phase. But the denominator changes sign as you cross from f<f0 to f>f0. Hence the 180 degree drops.
  2 Comments
David Goodmanson
David Goodmanson on 12 Nov 2021
Hello Chunru,
Whoops, I meant to address my answer to the original questioner, JL, and have modified my answer accordingly.
In the 512 point case, though, I don't believe that computational error is involved or that the results might be unpredictable. In the 512 point case the magnitude in the frequency domain has a maximum of 239 and a minimum of .288, only 830 times less and not near zero. As you say, the phase can change quickly around the peaks, but here the phase is quite smooth in the domain between the drops. I believe that while the phase curve may not be obvious, the phases are correct at all frequencies.

Sign in to comment.

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!