Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
how to interpret phase spectrum of of FFT of sinusoid?

Subject: how to interpret phase spectrum of of FFT of sinusoid?

From: shinchan

Date: 17 Nov, 2009 23:32:07

Message: 1 of 3

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

Subject: how to interpret phase spectrum of of FFT of sinusoid?

From: Claire Mulcock

Date: 18 Nov, 2009 00:01:54

Message: 2 of 3

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.

Subject: how to interpret phase spectrum of of FFT of sinusoid?

From: shinchan

Date: 18 Nov, 2009 01:35:04

Message: 3 of 3

Claire Mulcock <c.mulcock@gmail.com> wrote in message <b6623b2c-81b4-41a9-bd48-0f8c6ffda0cf@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?

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us