The fourier series coefficient phase

I am trying to get the x[n] from fourier coefficients, then I used the fft to find the coefficients of my x[n] and compare it with the prompt. So Here is the code

N = 8; % the period
figure;
k = -4:3; % the index for fourier transform
a_k = [0,0,1,1,1,1,1,0]; % fourier series
subplot(2,2,1);
stem(k,abs(a_k)); % plot the magnitude of the fourier coefficient
xlabel('k'); ylabel('|a_k|'); title('Magnitude of the Fourier Series Coefficient Directly From Prompt'); %xlabel tilte and ylabel
subplot(2,2,3);
stem(k,angle(a_k)); % plot the phase of the fourier coefficient
xlabel('k'); ylabel('phase of a_k'); title('Phase of the Fourier Series Coefficient Directly Prompt');%xlabel title and ylabel

x = zeros(1,N);
n = -4:3; % the n matrix
w = 2*pi/N; % omega
A = [1,2,2,0,0,0,0,0]; % the cos wave amplitude
phi = [0,0,0,0,0,0,0,0]; % the phase
for i = 1:N
x = x+A(i)*cos((i-1)*w*n+phi(i));
end

ak_re = fft(fftshift(x)); % the coefficients  
ak_shifted = fftshift(ak_re); % shifted version
subplot(2,2,2);
stem(k,abs(ak_shifted)/N);% stem plot
xlabel('k'); ylabel('|a_k|'); title('Magnitude of the Fourier Series Coefficient by FFT'); %xlabel tilte and ylabel
subplot(2,2,4);
stem(k,angle(ak_shifted));% stem plot
xlabel('k'); ylabel('phase of a_k'); title('Phase of the Fourier Series Coefficient by FFT');%xlabel title and ylabel

x_re = N*fftshift(ifft(ifftshift(a_k))); % Rebuilt x[n]
figure;
subplot(1,2,1);
stem(n,x_re); % plot the original graph
xlabel('n'); ylabel('x[n]'); title('Plot of x[n] Rebuilt from Coefficients'); %xlabel tilte and ylabel
subplot(1,2,2);
stem(n,x); % stem plot
xlabel('n'); % xlabel of the x[n]
ylabel('x[n]'); % ylabel
title('Plot of x[n] Directly from Expression'); % title

I find out that x[n] plots are same, but a_k plots are a little bit off on phase plots. Apparently the Phase of the Fourier Series Coefficient by FFT is not equal to the original phase of the coefficients by prompt. They are like half periods off. So I tried to use fftshift, but it did not help. So please let me know what is wrong here.

Answer 1

[David Goodmanson]

Hi Zeyuan

>> [ak_shifted; angle(ak_shifted)]

0.0000   -0.0000   8.0000   8.0000   8.0000   8.0000   8.0000   -0.0000      % value
0         3.1416        0        0        0        0        0    3.1416      % angle 

>> ak_shifted([2 8])
ans =  1.0e-15 *
       -0.8882   -0.8882

For the 'bad' phases you are finding the angle of quantities that are only different from zero because of numerical imprecision. They so happen to be negative, so you get a misleading result. Down in this region of size (precisionwise), angles don't have any meaning.

>> ak_shifted(1)
ans =   8.8818e-16

The first element of ak_shifted comes out tiny but positive, so you get an answer that looks great but in reality is just as unreliable as the other two.

[Walter Roberson]

By the way,

w = 2*pi/N; % omega

A = [1,2,2,0,0,0,0,0]; % the cos wave amplitude

phi = [0,0,0,0,0,0,0,0]; % the phase

for i = 1:N

x = x+A(i)*cos((i-1)*w*n+phi(i));

end

you are using cos() with something that is a multiple of numeric pi. cos() with numeric pi can be inaccurate:

cos(pi/2)

ans = 6.1232e-17

whereas

Pi = sym(pi)

Pi = 

π

cos(Pi/2)

ans = 

0

You can in theory compensate for the inaccuracy by using

W = 2/N; % omega

A = [1,2,2,0,0,0,0,0]; % the cos wave amplitude

phi = [0,0,0,0,0,0,0,0]; % the phase

for i = 1:N

x = x+A(i)*cospi((i-1)*W*n+phi(i)/pi);

end

but then the phi(i)/pi introduces theoretical inaccuracies. So you can in theory replace cos(a+b) by cos(a)*cos(b) - sin(a)*sin(b) and then in this context that would become cospi((i-1)*W*n).*cos(phi(i)) + sinpi((i-1)*W*n)*sin(phi(i)) . This pushes the inaccuracies to the representation of phi(i)