Fourier Analysis - Magnitude spectrum and phase shift
Show older comments
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;
Accepted Answer
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, /
7 Comments
Stefan
on 6 Apr 2024
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
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)
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)]
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
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
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)
%{
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)]
%% 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
on 8 Apr 2024
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
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)))
Now the imaginary part is small, but the real part is too large.
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)))
Because the imaginary part is rounding errors, get rid of it
f_x = real(f_x)
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]);
%% 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
More Answers (0)
Categories
Find more on Spectral Measurements in Help Center and File Exchange
