how to interpret phase spectrum of of FFT of sinusoid?

[shinchan]

Hi,
I would appreciate your input enormously.
I am trying to understand phase spectrum so I wrote a little code here, but I am puzzled by the output. My time domain signal is a combination of 50Hz and 120Hz sinusoids, with the phi term of 0.25pi and 0.75pi, respectively. So I expect my phase spectrum to show 0.25 and 0.75 rad at these frequencies respectively. However the code I have below doesn't do so, and I don't know why? I thought angle() gives me the phase of the cosine at that frequency. But I don't know how to interpret the phase spectrum.

2. What if I replace cosine with sine in my x? how do I interpret the phase spectrum then?

Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*cos(2*pi*50*t+0.25*pi) + cos(2*pi*120*t+0.75*pi);
y = x;
plot(Fs*t(1:50),y(1:50))
title('Signal')
xlabel('time (milliseconds)')

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
figure(2); subplot(211); plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

figure(2); subplot(212);
phi=angle(Y);
figure(2); plot(f, phi(1:NFFT/2+1)/pi); %plot in radian.
title('Single-Sided Phase Spectrum of y(t)');
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');
Y(52) %Magnitude at 50Hz.
Y(124) %Magnitude at 120Hz.
angle(Y(52))/pi %Phase angle (rad) at 50Hz
angle(Y(124))/pi %Phase angle (rad) at120Hz

[Claire Mulcock]

Yes, but what are f(52) and f(152).
They are not exactly 50 and 120 Hz, are they?
So, why do you expect the phases to be what you set them for 50 and
120 Hz?
Your error is setting NFFT to be a dyadic number.
Try this:
NFFT=length(y);
Modern FFTs no longer require the data length to be a dyadic number.
In your case, you have introduced leakage by doing this.

[shinchan]

Claire Mulcock <c.mu...@gmail.com> wrote in message <b6623b2c-81b4-41a9...@t11g2000prh.googlegroups.com>...

> On Nov 18, 12:32?pm, "shinchan " <shinchan75...@gmail.com> wrote:
> > Hi,
> > I would appreciate your input enormously.

> > I am trying to understand phase spectrum so I wrote a little code here, but I am puzzled by the output. My time domain signal is a combination of 50Hz and 120Hz sinusoids, with the phi term of 0.25pi and 0.75pi, respectively. So I expect my phase spectrum to show 0.25 and 0.75 rad at these frequencies respectively. However the code I have below doesn't do so, and I don't know why? I thought angle() ?gives me the phase of the cosine at that frequency. But I don't know how to interpret the phase spectrum.
> >
> > 2. What if ?I replace cosine with sine in my x? how do I interpret the phase spectrum then?
> >
> > Fs = 1000; ? ? ? ? ? ? ? ? ? ?% Sampling frequency
> > T = 1/Fs; ? ? ? ? ? ? ? ? ? ? % Sample time
> > L = 1000; ? ? ? ? ? ? ? ? ? ? % Length of signal
> > t = (0:L-1)*T; ? ? ? ? ? ? ? ?% Time vector

> > % Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
> > x = 0.7*cos(2*pi*50*t+0.25*pi) + cos(2*pi*120*t+0.75*pi);
> > y = x;
> > plot(Fs*t(1:50),y(1:50))
> > title('Signal')
> > xlabel('time (milliseconds)')
> >
> > NFFT = 2^nextpow2(L); % Next power of 2 from length of y
> > Y = fft(y,NFFT)/L;
> > f = Fs/2*linspace(0,1,NFFT/2+1);
> >
> > % Plot single-sided amplitude spectrum.
> > figure(2); subplot(211); plot(f,2*abs(Y(1:NFFT/2+1)))
> > title('Single-Sided Amplitude Spectrum of y(t)')
> > xlabel('Frequency (Hz)')
> > ylabel('|Y(f)|')
> >
> > figure(2); subplot(212);
> > phi=angle(Y);
> > figure(2); plot(f, phi(1:NFFT/2+1)/pi); %plot in radian.
> > title('Single-Sided Phase Spectrum of y(t)');
> > xlabel('Frequency (Hz)');
> > ylabel('Phase (rad)');
> > Y(52) %Magnitude at 50Hz.
> > Y(124) %Magnitude at 120Hz.

> > angle(Y(52))/pi ?%Phase angle (rad) at 50Hz

> > angle(Y(124))/pi %Phase angle (rad) at120Hz
>
> Yes, but what are f(52) and f(152).
> They are not exactly 50 and 120 Hz, are they?
> So, why do you expect the phases to be what you set them for 50 and
> 120 Hz?
> Your error is setting NFFT to be a dyadic number.
> Try this:
> NFFT=length(y);
> Modern FFTs no longer require the data length to be a dyadic number.
> In your case, you have introduced leakage by doing this.

Hi Claire,
Thank you very much. I can see now that the 51th and 121th elements of the phase vector indeed contains the correct phase, which is 0.25 and 0.75 respectively. But how about the rest of the phase spectrum, they look very noisy, is it because the magnitude terms are too small and renders phase meaningless? If so, how do I determine a fair threshold so I don't see those meaningless phase values?