# Fourier Analysis - Magnitude spectrum and phase shift

124 views (last 30 days)
Stefan on 5 Apr 2024
Edited: Stefan on 9 Apr 2024
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;

Paul on 5 Apr 2024
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, /
Paul on 8 Apr 2024
Edited: Paul on 9 Apr 2024
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 on 9 Apr 2024
Edited: Stefan on 9 Apr 2024
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