How can I correctly plot phase spectrum of fourier series with matlab?

Question

I've made a rectangular pulse x, with fundamental period 2seconds, to fourier transform

t=-2:0.01:2; 
dt=t(2)-t(1); %increment of time
fs=1/dt % sampling rate 
n=length(t) %number of samples
X=fftshift(fft(x,n))/n %fourier transform of x%

f=linspace(-fs/2,fs/2, n) %making frequency axis

X_angle=angle(X); %the phase of X

I was expecting the phase spectrum alternates -pi/2 and pi/2,

but the graph(too bad that I cannot post it due to lack of my reputation)

shows me that X_angle gradually increases as the frequency increases, ranges from -pi to pi

it just worked fine when I plotted magnitude spectrum, with

plot(f,X_mag), X_mag=abs(X)

I think I have to make some magnitude adjustment with angle(X)

but isn't angle irrelevant with magnitude of X?

I don't know why the absolute value of X_angle increases as the absolute frequency increases.

Answer 1

You made a couple of mistakes in your code.

1. You did not print out the values for x in your code. I assume you want to analyse a symmetric rectangular pulse. However your number of samples is an odd value. This means, your signal is not symmetric, which leads to small differences compared to the theoretical text book values.

   2. A rectangular pulse has an infinite bandwidth. However, for a FFT the sampling rate has to be at least twice the highest frequency of the signal. This means your sampling rate has to be at least 2*infinity. Bad luck, you can't do that. Nobody can do that. As a result you will get aliasing, which means that your results contain errors. The good news is, that in the case of a rectangular function this effect can be compensated with the help of the sinc function.

   3. If you made anything right, you would get the correct textbook coefficients. If your input signal has the form like (N=10) 11111-1-1-1-1-1, the function is an odd function. This means f(-t) = -f(t). In this case the rectangle can be constructed by a series of sine functions. The theoretical function is:

   f(t) = 4/pi( sin(wt) + 1/3 sin(3wt) + 1/5 sin(5wt) + 1/7 sin(7wt) ...)

   There are no frequencies with even numbers like sin(2wt) or sin(4wt). This means, every second frequency in the spectrum has the value zero. Caused by numerical noise these values are not exactly zero, but close to zero. Calculating a phase from these values produces meaningless values. The other frequencies are the Fourier transformed values of sine functions. The FFT of a sine function is pure imaginary. Because all elements of the sum have the same sign, all angles have the same value, which is pi/2.

Below you find the modified code:

close all;
clear all;
clc;

t=-2:0.01:(2-0.01); 
dt=t(2)-t(1); %increment of time
fs=1/dt; % sampling rate 
n=length(t); %number of samples
x = [ones(1,n/2), -ones(1,n/2)];
X=fft(x)/n; %fourier transform of x%

f=0:n-1; %making frequency axis
sincComp = @(N) (exp(-i*pi*(f/N)).*sinc(f/N)).';
%  Perform the compensation
comp = transpose(sincComp(n));
Xcorrected = X.*comp;

X_angle=angle(Xcorrected(2:2:n)); %the phase of X

figure;
subplot(3,1,1);
hold on;
stem(f(1:10), pi*abs(X(1:10))/2, 'ob', 'LineWidth', 3);
plot(f(1:10), pi*abs(Xcorrected(1:10))/2, 'or', 'LineWidth', 3);
hold off;
grid on;
title('Absolute value of FFT result', 'FontSize', 18);
xlabel('frequency', 'FontSize', 18);
ylabel('abs', 'FontSize', 18);
legend(['FFT'], ['FFT compensated']);
subplot(3,1,2);
hold on;
stem(f(1:10), pi*real(X(1:10))/2, 'ob', 'LineWidth', 3);
plot(f(1:10), pi*real(Xcorrected(1:10))/2, 'or', 'LineWidth', 3);
hold off;
grid on;
title('Real value of FFT result', 'FontSize', 18);
xlabel('frequency', 'FontSize', 18);
ylabel('real', 'FontSize', 18);
subplot(3,1,3);
hold on;
stem(f(1:10), pi*imag(X(1:10))/2, 'ob', 'LineWidth', 3);
plot(f(1:10), pi*imag(Xcorrected(1:10))/2, 'or', 'LineWidth', 3);
hold off;
grid on;
title('Imaginary value of FFT result', 'FontSize', 18);
xlabel('frequency', 'FontSize', 18);
ylabel('imag', 'FontSize', 18);

Result of FFT and zinc-compensated FFT

The phase values of the first ten non-zero frequencies:

>> X_angle(1:10)*2/pi
ans =
-1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000   -1.0000