Fourier Analysis - Magnitude spectrum and phase shift

Hi All,

After numerous attempts to plot the correct values for Magnitude and phase shift of triangular wave signal. The approximation has been calculated using Complex fourier series and the result is plugged in the code. the wave form looks pretty good depending the on n=... I have attempted to plot analysis of coeficients but, its not near what i am getting on the paper. The Amplitude of first harmonic seems alright, but multiiples should decay, and according to my knowlage and calculations on paper, there should be no phase shift. I have worked on the code applying for loops, but no luck, I wonder if someone with more MATLAB coding experiance could throw an eye over it and snagg it up. Any advise is much appriciated

Best

Stefan

%% Triangle wave Fourier approximation + Phase & Magnitude Spectra

%% Setting Maximum Harmonic to be used (important to be odd number)

h = 11;

%% Defining period of the function

p = 2;

%% Creating an array for x from -11 to 11

n_Pos = 1:1:h;

n_Neg = -h:1:-1;

n_f = [n_Neg,0,n_Pos];

%% Define the range of x values to evaluate the function

S = linspace(-h, h, 1000);

%% Defining the triangle wave function

f = @(S) sawtooth(pi*S+3/2,1/2);

%% Evaluate the function

y = f(S);

%% processing amplitude spectra

mag_Cn_0 = 0 ; % at 0

mag_Cn_Pos = (((4*1i*(-1).^n_Pos)/(pi.^2*n_Pos.^2))*sin((pi*n_Pos)/2)); % for >0

mag_Cn_Neg = (((4*1i*(-1).^n_Neg)/(pi.^2*n_Neg.^2))*sin((pi*n_Neg)/2)); % for <0

mag_Cn = [mag_Cn_Neg mag_Cn_0 mag_Cn_Pos]; % combined

figure;

name='Fourier series of Triangle waveform';

subplot(311);

set(gcf, 'Name', name);

stem(n_f, abs(mag_Cn));

title("Amplitude Spectra");

xlabel("Fourier Coefficients");

ylabel("Magnitude");

grid on;

xlim([0 11]);

%% defining variables for triangle wave reconstruction

w0 = 2*pi/p; % wo definition using the function's Period

q = -2:0.01:11; % Reconstruction data for processing

f_x = zeros(1, length(q)); % initialising output array

c = @(nf) ((4*1i*(-1)^nf)/(pi^2*nf^2))*sin((pi*nf)/2); % c(n) for processing

%% generating fourier Approximation to h'th harmonic

for n = 1:1:3

f_x = f_x+(c(n)*exp(1i*w0*n.*q))-(c(n)*exp(-1i*w0*n.*q));

end

f_x = mag_Cn_0 + f_x;

%% processing phase spectra

mag_Cn(abs(imag(mag_Cn))<0.0001) = 0; % replace real + imag 0 as real 0 only

subplot(312);

stem(n_f, angle(mag_Cn)); % plot phase Spectrum using angle function

title("Phase Spectra");

xlabel("Fourier Coefficients");

ylabel("Phase (pi)");

xlim([0 11]);

ylim([-2 2]);

grid on;

subplot(313);

hold on;

plot(q, f_x,'linestyle','-.',"LineWidth",2.5);

plot(S, y,'linestyle','--',"LineWidth",2);

title("Triangle wave & Fourier Reconstruction");

xlabel("x");

ylabel("f(x)");

xlim([-0.1 11]);

ylim([-1.2 1.2]);

legend("Fourier Transform","Rect Pulse Train")

grid on;

Answer 1

[Paul]

Hi Stefan,

%% Triangle wave Fourier approximation + Phase & Magnitude Spectra
%% Setting Maximum Harmonic to be used (important to be odd number) 
h = 11;
%% Defining period of the function
p = 2;
%% Creating an array for x from -11 to 11
n_Pos = 1:1:h;
n_Neg = -h:1:-1;
n_f = [n_Neg,0,n_Pos];
%% Define the range of x values to evaluate the function
S = linspace(-h, h, 1000);
%% Defining the triangle wave function
f = @(S) sawtooth(pi*S+3/2,1/2);
%% Evaluate the function
y = f(S);

Plotting one period of the function shows that f isn't being defined quite correctly, though I think it's only used for plotting.

figure
plot(S,y),xlim([0 2])

Here are the equations for the Fourier series coefficients.

%% processing amplitude spectra 
mag_Cn_0 = 0 ;                                                          % at 0 
mag_Cn_Pos = (((4*1i*(-1).^n_Pos)/(pi.^2*n_Pos.^2))*sin((pi*n_Pos)/2));    % for >0 
mag_Cn_Neg = (((4*1i*(-1).^n_Neg)/(pi.^2*n_Neg.^2))*sin((pi*n_Neg)/2));    % for <0
mag_Cn = [mag_Cn_Neg mag_Cn_0 mag_Cn_Pos]; 

n_Pos and n_Neg are vectors. I suggest reviewing the equations above, particularly accounting for the difference between element wise multiplication and division as implemented by times, .* and rdivide, ./ respectively, versus mtimes, * and mrdivide, /

[Stefan]

Hi Paul,

Thanks for highlighting the issues. I have gone through them and I made some improvements.

see the triangle wave plot, much better than previous definition of "f" function

I have improved the code as follow:

%% Triangle wave Fourier approximation + Phase & Magnitude Spectra

%% Setting Maximum Harmonic to be used (important to be odd number)

h = 11;

%% Defining period of the function

p = 2;

%% Creating an array for x from -11 to 11

n_Pos = 1:1:h;

n_Neg = -h:1:-1;

n_f = [n_Neg,0,n_Pos];

%% Define the range of x values to evaluate the function

S = linspace(-h, h, 1000);

%% Defining the triangle wave function

f = @(S) triangleWave(S,1,p,0);

%% Evaluate the function

y = f(S);

%% processing amplitude spectra

mag_Cn_0 = 0; % at 0

mag_Cn_Pos = sqrt(((4*(-1).^n_Pos).^2)+1).*sin((pi*n_Pos)/2)./...

((pi^2)*(n_Pos.^2)); % for >0

mag_Cn_Neg = sqrt(((4*(-1).^n_Neg).^2)+1).*sin((pi*n_Neg)/2)./...

((pi^2)*(n_Neg.^2)); % for <0

mag_Cn = [mag_Cn_Neg mag_Cn_0 mag_Cn_Pos]; % combined

figure;

name='Fourier series of Triangle waveform';

subplot(311);

set(gcf, 'Name', name);

stem(n_f, abs(mag_Cn));

title("Amplitude Spectra");

xlabel("Fourier Coefficients");

ylabel("Magnitude");

grid on;

xlim([0 11]);

%% defining variables for triangle wave reconstruction

w0 = 2*pi/p; % wo definition using the function's Period

q = -11:0.01:11; % Reconstruction data for processing

f_x = zeros(1, length(q)); % initialising output array

c = @(n_f) ((4*1i*(-1)^n_f)/(pi^2*n_f^2))*sin((pi*n_f)/2); % c(n) for processing

%% generating fourier Approximation to h'th harmonic

for n = 1:1:3

f_x = f_x+(c(n)*exp(1i*w0*n.*q))-(c(n)*exp(-1i*w0*n.*q));

end

f_x = mag_Cn_0 + f_x;

%% processing phase spectra

mag_Cn(abs(imag(mag_Cn))<0.0001) = 0; % replace real + imag 0 as real 0 only

subplot(312);

stem(n_f, angle(mag_Cn)); % plot phase Spectrum using angle function

title("Phase Spectra");

xlabel("Fourier Coefficients");

ylabel("Phase (pi)");

xlim([0 11]);

ylim([-2 2]);

grid on;

subplot(313);

hold on;

plot(q, f_x,'linestyle','-.',"LineWidth",2.5);

plot(S, y,'linestyle','--',"LineWidth",2);

title("Triangle wave & Fourier Reconstruction");

xlabel("x");

ylabel("f(x)");

xlim([-0.1 11]);

ylim([-1.2 1.2]);

legend("Fourier Transform","Triangle wave")

grid on;

%% Triangle_wave function

function [ y ] = triangleWave( time, amplitude, waveL, phaseShift)

y = 4*amplitude/waveL *( abs( mod((time-waveL/4+phaseShift),waveL) - waveL/2) - waveL/4);

end

resulting plot shows exact ampliotudes as calculated on paper, i have rearanged the eqaution to calculate them, with sqrt()

here is the manual calculation for comper:

Many thanks for help and support

Best

Stefan

[Paul]

The code still doesn't look right

%% Triangle wave Fourier approximation + Phase & Magnitude Spectra
%% Setting Maximum Harmonic to be used (important to be odd number) 
h = 11;
%% Defining period of the function
p = 2;
%% Creating an array for x from -11 to 11
n_Pos = 1:1:h;
n_Neg = -h:1:-1;
n_f = [n_Neg,0,n_Pos];
%% Define the range of x values to evaluate the function
S = linspace(-h, h, 1000);
%% Defining the triangle wave function
f = @(S) triangleWave(S,1,p,0);   
%% Evaluate the function
y = f(S);

Looks like these equations changed from the original code

%% processing amplitude spectra 
mag_Cn_0 = 0;                                                       % at 0 
mag_Cn_Pos = sqrt(((4*(-1).^n_Pos).^2)+1).*sin((pi*n_Pos)/2)./...
                ((pi^2)*(n_Pos.^2));                                % for >0
mag_Cn_Neg = sqrt(((4*(-1).^n_Neg).^2)+1).*sin((pi*n_Neg)/2)./...
                ((pi^2)*(n_Neg.^2));                                % for <0
mag_Cn = [mag_Cn_Neg mag_Cn_0 mag_Cn_Pos];                          % combined

Look at the first few values of mag_Cn_pos

mag_Cn_Pos(1:5)

ans = 1x5

0.4178 0.0000 -0.0464 -0.0000 0.0167

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure;
name='Fourier series of Triangle waveform';
subplot(311);
set(gcf, 'Name', name);
stem(n_f, abs(mag_Cn));
title("Amplitude Spectra");
xlabel("Fourier Coefficients");
ylabel("Magnitude");
grid on;
xlim([0 11]);
%% defining variables for triangle wave reconstruction
w0 = 2*pi/p;                    % wo definition using the function's Period 
q = -11:0.01:11;                % Reconstruction data for processing  
f_x = zeros(1, length(q));      % initialising output array 
c = @(n_f) ((4*1i*(-1)^n_f)/(pi^2*n_f^2))*sin((pi*n_f)/2);   % c(n) for processing
%% generating fourier Approximation to h'th harmonic
for n = 1:1:3
    f_x = f_x+(c(n)*exp(1i*w0*n.*q))-(c(n)*exp(-1i*w0*n.*q));
end
f_x = mag_Cn_0 + f_x;

Check the first couple of values c(n)

[c(1) c(2) c(3) c(4) c(5)]

ans = 
   0.0000 - 0.4053i   0.0000 + 0.0000i   0.0000 + 0.0450i   0.0000 - 0.0000i   0.0000 - 0.0162i

These don't match the values of mag_Cn_Pos, even if the latter is supposed to be interpreted as a magnitude (based on the variable's name), even though mag_Cn_Pos can take on a negative value.

%% processing phase spectra

This line sets all of mag_Cn to 0, though it looks like the intent is only to get rid of small rounding errors in the imaginary part.

mag_Cn(abs(imag(mag_Cn))<0.0001) = 0;    % replace real + imag 0 as real 0 only 
subplot(312);
stem(n_f, angle(mag_Cn));                 % plot phase Spectrum using angle function
title("Phase Spectra");
xlabel("Fourier Coefficients");
ylabel("Phase (pi)");
xlim([0 11]);
ylim([-2 2]);
grid on;
subplot(313);
hold on;
plot(q, f_x,'linestyle','-.',"LineWidth",2.5);
plot(S, y,'linestyle','--',"LineWidth",2);
title("Triangle wave & Fourier Reconstruction");
xlabel("x");
ylabel("f(x)");
xlim([-0.1 11]);
ylim([-1.2 1.2]);
legend("Fourier Transform","Triangle wave")
grid on;

%% Triangle_wave function
function [ y ] = triangleWave( time, amplitude, waveL, phaseShift)
y = 4*amplitude/waveL *( abs( mod((time-waveL/4+phaseShift),waveL) - waveL/2) - waveL/4);
end

[Stefan]

Hi Paul,

After revising the paper calculation and the code, I have found few mistakes. these has been addressed. Pls see the plot. in fact the previos magnitude of fundamental frequency was only approaching 0.5. Now its at 0.8, where if we add all remining frequences the approximation is reaching peak at 1 on y axis.

regarding the phase shift, here is my calculation , I based my code on. In fact the frequencies contributing in to function approximation they do obtain phase shift as shown on the plot.

If you think, the task/code could be improved, I'd be happy to follow the lead from experianced one. I am still just learning how to work with Matlab and Fourier series subject.

Thanks

Stefan

%% Triangle wave Fourier Complex series approximation + Phase & Magnitude Spectra

%% Setting Maximum Harmonic to be used (important to be odd number)

h = 11;

%% Defining period of the function

p = 2;

%% Creating an array for x from -11 to 11

n_Pos = 1:1:h;

n_Neg = -h:1:-1;

n_f = [n_Neg 0 n_Pos];

%% Define the range of x values to evaluate the function

S = linspace(-h, h, 1000);

%% Defining the triangle wave function

% triangle atributes "time, amplitude, waveL, phaseShift"

f = @(S) triangleWave(S,1,p,0);

%% Evaluate the function

y = f(S);

%% processing amplitude spectra

mag_Cn_0 = 0; % at 0

mag_Cn_Pos = sqrt((8.^2)+(1i)^2).*sin((pi*n_Pos)/2)./...

((pi^2)*(n_Pos.^2)); % for > 0

mag_Cn_Neg = sqrt((8.^2)+(1i)^2).*sin((pi*n_Neg)/2)./...

((pi^2)*(n_Neg.^2)); % for < 0

mag_Cn = [mag_Cn_Neg mag_Cn_0 mag_Cn_Pos]; % combined

figure;

name='Fourier series of Triangle waveform';

subplot(311);

set(gcf, 'Name', name);

stem(n_f, abs(mag_Cn)); % plot magnitudes by aplying absolute number function

title("Amplitude Spectra");

xlabel("Fourier Coefficients");

ylabel("Magnitude");

grid on;

xlim([0 11]);

%% defining variables for triangle wave reconstruction

w0 = 2*pi/p; % wo definition using the function's Period

q = -11:0.01:11; % Reconstruction data for processing

f_x = zeros(1, length(q)); % initialising output array

c = @(n_f) ((8*1i*(-1)^n_f)./((pi^2)*(n_f.^2)))*sin((pi.*n_f)./2); % c(n) for processing

%% generating fourier Approximation to h'th harmonic

for n = 1:1:10

f_x = f_x+(c(n)*exp(1i*w0*n.*q));

end

f_x = mag_Cn_0 + f_x;

%% processing phase spectra

phase_Pos = 1i./(((8*(-1).^n_Pos))./(pi^2.*n_Pos.^2)).*sin((pi.*n_Pos)/2);

phase_Neg = 1i./(((8*(-1).^n_Neg))./(pi^2.*n_Neg.^2)).*sin((pi.*n_Neg)/2);

phase_Cn = [phase_Neg mag_Cn_0 phase_Pos] ;

phase_Cn(abs(imag(phase_Cn))<0.001) = 0; % Filter vey small phase shift calculations near 0

phase_Cn = angle(phase_Cn); % aply conversion to phase-shift in radians

subplot(312);

stem(n_f, phase_Cn); % plot phase Spectrum using angle function a

title("Phase Spectra");

xlabel("Fourier Coefficients");

ylabel("Phase (pi)");

xlim([0 11]);

ylim([-2 2]);

grid on;

subplot(313);

hold on;

plot(q, f_x,'linestyle','-.',"LineWidth",2.5);

plot(S, y,'linestyle','--',"LineWidth",2);

title("Triangle wave & Fourier Reconstruction");

xlabel("x");

ylabel("f(x)");

xlim([-0.1 11]);

ylim([-1.2 1.2]);

legend("Fourier Transform","Triangle wave")

grid on;

%% Triangle_wave function

function [ k ] = triangleWave( time, amplitude, waveL, phaseShift)

k = 4*amplitude/waveL *( abs( mod((time-waveL/4+phaseShift),waveL) - waveL/2) - waveL/4);

end

[Paul]

Hi Stefan,

mag_Cn_pos still does not match the c(n_f).

I see now you've multiplied c(n_f) by a factor of two, and then only used the positive values of n to construct fx. I don't know if that approach is correct or not, but it's not they way I would do it, at least not w/o checking it thoroughly. Also, it looks like you've done the same scaling on mag_Cn*, but still computing mag_Cn* for negative values of n, which seems unnecssary if you're trying to get a series that only needs n >= 0. Also, there's no need to compute mag_Cn for positive and negative n separately, at least not as far as I can tell.

My suggestion is to not do the multiply by 2. Just compute the c(n_f) for all n, reconstruct f using positive and negative values of n, and compute the mag and phase from c(n_f).

It looks like you're solving for c(n_f) by hand and then from that computing by hand the magnitude and phase of c(n_f). Unless the closed-form expressions for the mag and phase is a requirement of the problem, I'd just stick with evaluating c(n_f) at values of n and then let Matlab compute the mag and phase of the results using abs and angle.

%% Triangle wave Fourier Complex series approximation + Phase & Magnitude Spectra
%% Setting Maximum Harmonic to be used (important to be odd number) 
h = 11;
%% Defining period of the function
p = 2;
%% Creating an array for x from -11 to 11
n_Pos = 1:1:h;
n_Neg = -h:1:-1;
n_f = [n_Neg 0 n_Pos];
%% Define the range of x values to evaluate the function
S = linspace(-h, h, 1000);
%% Defining the triangle wave function
%  triangle atributes "time, amplitude, waveL, phaseShift"
f = @(S) triangleWave(S,1,p,0);   
%% Evaluate the function
y = f(S);
%% processing amplitude spectra 
mag_Cn_0 = 0;                                          % at 0 
mag_Cn_Pos = sqrt((8.^2)+(1i)^2).*sin((pi*n_Pos)/2)./...
                ((pi^2)*(n_Pos.^2));                   % for > 0   
mag_Cn_Neg = sqrt((8.^2)+(1i)^2).*sin((pi*n_Neg)/2)./...
                ((pi^2)*(n_Neg.^2));                   % for < 0
mag_Cn = [mag_Cn_Neg mag_Cn_0 mag_Cn_Pos];             % combined   

mag_Cn_Pos(1:5)

ans = 1x5

0.8042 0.0000 -0.0894 -0.0000 0.0322

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%{
figure;
name='Fourier series of Triangle waveform';
subplot(311);
set(gcf, 'Name', name);
stem(n_f, abs(mag_Cn));     % plot magnitudes by aplying absolute number function
title("Amplitude Spectra");
xlabel("Fourier Coefficients");
ylabel("Magnitude");
grid on;
xlim([0 11]);
%}

%% defining variables for triangle wave reconstruction
w0 = 2*pi/p;                    % wo definition using the function's Period 
q = -11:0.01:11;                % Reconstruction data for processing  
f_x = zeros(1, length(q));      % initialising output array 
c = @(n_f) ((8*1i*(-1)^n_f)./((pi^2)*(n_f.^2)))*sin((pi.*n_f)./2);   % c(n) for processing

[c(1) c(2) c(3) c(4) c(5)]

ans = 
   0.0000 - 0.8106i   0.0000 + 0.0000i   0.0000 + 0.0901i   0.0000 - 0.0000i   0.0000 - 0.0324i

%% generating fourier Approximation to h'th harmonic
for n = 1:1:10
    f_x = f_x+(c(n)*exp(1i*w0*n.*q));
   
end
f_x = mag_Cn_0 + f_x;
%% processing phase spectra 
phase_Pos = 1i./(((8*(-1).^n_Pos))./(pi^2.*n_Pos.^2)).*sin((pi.*n_Pos)/2);
phase_Neg = 1i./(((8*(-1).^n_Neg))./(pi^2.*n_Neg.^2)).*sin((pi.*n_Neg)/2);
phase_Cn = [phase_Neg mag_Cn_0 phase_Pos]  ;
phase_Cn(abs(imag(phase_Cn))<0.001) = 0;   % Filter vey small phase shift calculations near 0 
phase_Cn = angle(phase_Cn);             % aply conversion to phase-shift in radians

%{
subplot(312);
stem(n_f, phase_Cn);                 % plot phase Spectrum using angle function a
title("Phase Spectra");
xlabel("Fourier Coefficients");
ylabel("Phase (pi)");
xlim([0 11]);
ylim([-2 2]);
grid on;
subplot(313);
hold on;
plot(q, f_x,'linestyle','-.',"LineWidth",2.5);
plot(S, y,'linestyle','--',"LineWidth",2);
title("Triangle wave & Fourier Reconstruction");
xlabel("x");
ylabel("f(x)");
xlim([-0.1 11]);
ylim([-1.2 1.2]);
legend("Fourier Transform","Triangle wave")
grid on;
%}

%% Triangle_wave function
function [ k ] = triangleWave( time, amplitude, waveL, phaseShift)
k = 4*amplitude/waveL *( abs( mod((time-waveL/4+phaseShift),waveL) - waveL/2) - waveL/4);
end

[Stefan]

Hi Paul,

You are right, there is unnecessary stuff, which has been removed. I have used one off formula calculating c(n_f) and then used MATLAB functions to extract Magnitude and phase shift from each point. I hope this keeps the sanity over the code. Also, it corected the phase shift plot. This was great learning curve in terms of subject and coding. Many thanks for taking part as The Devil's Advocate.

%% Triangle wave Fourier Complex series approximation + Phase & Magnitude Spectra

%% Setting Maximum Harmonic to be used (important to be odd number)

h = 11;

%% Defining period of the function

p = 2;

%% Creating an array for x from -11 to 11

n_Pos = 1:1:h;

n_Neg = -h:1:-1;

n_f = [n_Neg 0 n_Pos];

%% Define the range of x values to evaluate the function

S = linspace(-h, h, 1000);

%% Defining the triangle wave function

% triangle atributes "time, amplitude, waveL, phaseShift"

f = @(S) triangle_Wave(S,1,p,0);

%% Evaluate the function

y = f(S);

%% defining variables for triangle wave reconstruction

Cn_0 = 0;

w0 = 2*pi/p; % wo definition using the function's Period

q = -11:0.01:11; % Reconstruction data for processing

f_x = zeros(1, length(q)); % initialising output array

c = @(n_f) ((4*1i*(((-1).^n_f)-1))./((pi^2)*(n_f.^2))).*sin((pi.*n_f)/2); % c(n) for processing

%% generating fourier Approximation to h'th harmonic

for n = 1:1:10

f_x = f_x+(c(n)*exp(1i*w0*n.*q));

end

f_x = Cn_0 + f_x;

%% processing amplitude spectra

mag_Cn = hypot(abs(c(n_f)),imag(c(n_f)));

figure;

name='Fourier series of Triangle waveform';

subplot(311);

set(gcf, 'Name', name);

stem(n_f, mag_Cn); % plot magnitudes by aplying absolute number function

title("Amplitude Spectra");

xlabel("Fourier Coefficients");

ylabel("Magnitude");

grid on;

xlim([0 11]);

%% processing phase spectra

phase_Cn = imag(c(n_f))./real(c(n_f));

phase_Cn(phase_Cn < 0.0001) = 0; % Filter vey small phase shift calculations near 0

phase_Cn = atan(phase_Cn); % aply conversion to phase-shift in radians

%% Ploting

subplot(312);

stem(n_f, phase_Cn); % plot phase Spectrum using angle function a

title("Phase Spectra");

xlabel("Fourier Coefficients");

ylabel("Phase (pi)");

xlim([0 11]);

ylim([-2 2]);

grid on;

subplot(313);

hold on;

plot(q, f_x,'linestyle','-.',"LineWidth",2.5);

plot(S, y,'linestyle','--',"LineWidth",2);

title("Triangle wave & Fourier Reconstruction");

xlabel("x");

ylabel("f(x)");

xlim([-0.1 11]);

ylim([-1.2 1.2]);

legend("Fourier Transform","Triangle wave")

grid on;

%% Triangle_wave function for referance

function [ k ] = triangle_Wave( time, amplitude, waveL, phaseShift)

k = 4*amplitude/waveL *( abs( mod((time-waveL/4+phaseShift),waveL) - waveL/2) - waveL/4);

end

Regards

Stefan

[Paul]

I'm afraid there is still a problem with this updated code.

h = 11;
%% Defining period of the function
p = 2;
%% Creating an array for x from -11 to 11
n_Pos = 1:1:h;
n_Neg = -h:1:-1;
n_f = [n_Neg 0 n_Pos];
%% Define the range of x values to evaluate the function
S = linspace(-h, h, 1000);
%% Defining the triangle wave function
%  triangle atributes "time, amplitude, waveL, phaseShift"
f = @(S) triangle_Wave(S,1,p,0);   
%% Evaluate the function
y = f(S);
%% defining variables for triangle wave reconstruction
Cn_0 = 0;  
w0 = 2*pi/p;                    % wo definition using the function's Period 
q = -11:0.01:11;                % Reconstruction data for processing  
f_x = zeros(1, length(q));      % initialising output array 
c = @(n_f) ((4*1i*(((-1).^n_f)-1))./((pi^2)*(n_f.^2))).*sin((pi.*n_f)/2);   % c(n) for processing
%% generating fourier Approximation to h'th harmonic
for n = 1:1:10
    f_x = f_x+(c(n)*exp(1i*w0*n.*q));   
end
f_x = Cn_0 + f_x;

We see that the real part of f_x is at least reasonable, but the imaginary part is not zero, when it should be

figure
plot(q,real(f_x),q,imag(f_x)),grid,legend('real part','imaginary part')

Unless the code is trying to take a shortcut that I can't see, it should be summing from -N to N, like so

f_x = zeros(1, length(q));      % initialising output array 
%% generating fourier Approximation to h'th harmonic
for n = [-10:1:-1 1:1:10]
    f_x = f_x+(c(n)*exp(1i*w0*n.*q));   
end
f_x = Cn_0 + f_x;

figure
plot(q,real(f_x),q,imag(f_x)),grid,legend('real part','imaginary part')

max(abs(imag(f_x)))

ans = 6.3317e-17

Now the imaginary part is small, but the real part is too large.

Let's go back to this comment and use (changing the operators for elementwise operations)

c = @(n_f) ((4*1i*(-1).^n_f)./(pi^2*n_f.^2)).*sin((pi*n_f)/2);   % c(n) for processing

Then we get

f_x = zeros(1, length(q));      % initialising output array 
%% generating fourier Approximation to h'th harmonic
for n = [-10:1:-1 1:1:10]
    f_x = f_x+(c(n)*exp(1i*w0*n.*q));   
end
f_x = Cn_0 + f_x;

figure
plot(q,real(f_x),q,imag(f_x)),grid,legend('real part','imaginary part')

max(abs(imag(f_x)))

ans = 3.0983e-17

Because the imaginary part is rounding errors, get rid of it

f_x = real(f_x)

f_x = 1x2201

-0.0000 -0.0212 -0.0424 -0.0632 -0.0838 -0.1040 -0.1238 -0.1432 -0.1622 -0.1811 -0.1997 -0.2184 -0.2373 -0.2564 -0.2758 -0.2956 -0.3159 -0.3366 -0.3577 -0.3791 -0.4006 -0.4222 -0.4435 -0.4646 -0.4852 -0.5053 -0.5249 -0.5439 -0.5624 -0.5805

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Why not use abs(c(n_f)) here? And if using hypot, shouldn't the first argument be real(c(n_f))?

%% processing amplitude spectra 
%mag_Cn = hypot(abs(c(n_f)),imag(c(n_f))); 
mag_Cn = abs(c(n_f)); 

figure;
name='Fourier series of Triangle waveform';
subplot(311);
set(gcf, 'Name', name);
stem(n_f, mag_Cn);     % plot magnitudes by aplying absolute number function
title("Amplitude Spectra");
xlabel("Fourier Coefficients");
ylabel("Magnitude");
grid on;
xlim([0 11]);

Why not use angle to compute the phase?

%% processing phase spectra 
%phase_Cn = imag(c(n_f))./real(c(n_f));
%phase_Cn(phase_Cn < 0.0001) = 0;   % Filter vey small phase shift calculations near 0 
%phase_Cn = atan(phase_Cn);             % aply conversion to phase-shift in radians
phase_Cn = angle(c(n_f));
%% Ploting
subplot(312);
stem(n_f, phase_Cn);                 % plot phase Spectrum using angle function a
title("Phase Spectra");
xlabel("Fourier Coefficients");
ylabel("Phase (pi)");
xlim([0 11]);
ylim([-2 2]);
grid on;
subplot(313);
hold on;
plot(q, f_x,'linestyle','-.',"LineWidth",2.5);
plot(S, y,'linestyle','--',"LineWidth",2);
title("Triangle wave & Fourier Reconstruction");
xlabel("x");
ylabel("f(x)");
xlim([-0.1 11]);
ylim([-1.2 1.2]);
legend("Fourier Transform","Triangle wave")
grid on;

If you're interested, we can show how to use Matlab to verify our current expression for c(n_f).

%% Triangle_wave function for referance
function [ k ] = triangle_Wave( time, amplitude, waveL, phaseShift)
k = 4*amplitude/waveL *( abs( mod((time-waveL/4+phaseShift),waveL) - waveL/2) - waveL/4);
end

[Stefan]

Hi Paul,

Thanks for the feedback, and I think you got it right, As we know Fourier series is the sum of the equations from minus infinity to infinity, This duality property actually cancels out the imaginary part.

I have spent hours revising my very first handwritten calculation of c_n coficient, to match the result in Matlab. I thought I have mad mistake and still I have proved it right many times. but it did not match the reference triangle wave.

so, the second statement came later, but by plotting only the positive part of the sum of equations, therefore I had to compromise the handwritten calculation to satisfy the plot result.

Now, it all makes sense.

Regarding Matlab functions, I gave them a good look and, they are a gift. These functions are just a perfect fit for the Fourier series. On a side note, I aimed to show the examiner that I was using the right formulas for calculating magnitude and phase shift.

I have optimized the code, and it works a charm. Thank you

Regards

Stefan

PS: Great to see how you troubleshoot the code, I learned a lot, for sure. I am always up to improvement, I am in the middle of doing Matlab courses