# Matlab: trigonometric form of Fourier Series

I have tried to implement a matlab function that computes a Fourier series of a discrete periodic signal using its trigonometric form. The Fourier series can be obtained with:

\begin{align} s_{\scriptscriptstyle N}(x) = \frac{a_0}{2} + \sum_{n=1}^N \left(a_n \cos\left(\tfrac{2\pi}{T} nx \right) + b_n \sin\left(\tfrac{2\pi}{T} nx \right) \right). \end{align} \begin{align} a_0 &= \frac{1}{T}\int_{T} f(x) dx\\ b_0 &= 0\\ \end{align} for $$n>0$$: \begin{align} a_n &= \frac{2}{T}\int_{T} f(x)\cdot \cos\left( \tfrac{2\pi}{T} nx \right)\, dx\\ b_n &= \frac{2}{T}\int_{T} f(x)\cdot \sin\left( \tfrac{2\pi}{T} nx \right)\, dx \end{align}

I have therefor implemented the following Matlab function:

function fFS = get_fourier_trigo(f, T, t, dt, N, k_plot)
figure()
plot(t, f,'-k', 'LineWidth',1.5, 'DisplayName', 'Original')
hold on

% Compute Fourier series
A0 = (1/T)*sum(f)*dt % This corresponds to the mean value
fFS = A0/2;
for n=1:N
A(n) = (2/T)*sum(f.*cos(n*t*2*pi/T))*dt;
B(n) = (2/T)*sum(f.*sin(n*t*2*pi/T))*dt;
fFS = fFS + A(n)*cos(n*t*2*pi/T) + B(n)*sin(n*t*2*pi/T);
if mod(n, k_plot) == 0
plot(t, fFS, '-', 'LineWidth', 0.5, 'DisplayName', int2str(n))
end
end
title('Fourier Series')
legend show
end


But unfortunately my reconstructed Fourier series does not match my original signal:

dt = 1/100;
T = 1; % signal period
t = (-1+dt:dt:1)*T/2;
f = 2*sin(2*pi*t) + 0.5*sin(2*pi*5*t)
N = 5;
y_s = get_fourier_trigo(f, T, t, dt, N, 1);


Results of Fourier series approximation

I found out that if I modify my function by removing the $$\frac{1}{T}, \frac{2}{T}$$ for the estimation of $$a_0$$, $$a_n$$ and $$b_n$$ respectively then my obtained Fourier series converges to my original signal:

dt = 1/100;
T = 1; % signal period
t = (-1+dt:dt:1)*T/2;
f = 2*sin(2*pi*t) + 0.5*sin(2*pi*5*t)
N = 5;
y_s = get_fourier_trigo_corrected(f, T, t, dt, N, 1);

function fFS = get_fourier_trigo_corrected(f, T, t, dt, N, k_plot)
figure()
plot(t, f,'-k', 'LineWidth',1.5, 'DisplayName', 'Original')
hold on

% Compute Fourier series
A0 = sum(f)*dt % This corresponds to the mean value
fFS = A0/2;
for n=1:N
A(n) = sum(f.*cos(n*t*2*pi/T))*dt;
B(n) = sum(f.*sin(n*t*2*pi/T))*dt;
fFS = fFS + A(n)*cos(n*t*2*pi/T) + B(n)*sin(n*t*2*pi/T);
if mod(n, k_plot) == 0
plot(t, fFS, '-', 'LineWidth', 0.5, 'DisplayName', int2str(n))
end
end
title('Fourier Series')
legend show
end


Results of the modified Fourier series approximation

Can anyone explain me where is my mistake / misunderstanding?

Update: I have implemented the changes suggested by @jean marie and it works well for sum of sines waves but I now face an issue with obtaining the Fourier series of a pulse signal. I have also replaced the sums by a numerical integration using trapz:

dt = 1/1000; % sampling time
T = pi; % periodicity
t = (0:dt:1-dt)*T;

% split interval in quarters
n = length(t);
nquart = floor(n/4);

% Pulse function
f = 0*t;
f(nquart:3*nquart) = 1;

% Fourier approximation
N = 50;
kplot = 10;

fs = get_fourier_trigo(f, T, t, dt, N, kplot);
%% Functions
function fFS = get_fourier_trigo(f, T, t, dt, N, k_plot)
figure()
plot(t, f,'-k','LineWidth',1.5, 'DisplayName', 'Original')
hold on

% Compute Fourier series
A0 = (1/T)*trapz(t,f) % This corresponds to the mean value
fFS = A0/2;
wt = t*2*pi/T;
for n=1:N
A(n) = (2/T)*trapz(t ,f.*cos(n*wt));
B(n) = (2/T)*trapz(t, f.*sin(n*wt));
fFS = fFS + A(n)*cos(n*wt) + B(n)*sin(n*wt);
if mod(n, k_plot) == 0
plot(t, fFS, '-', 'LineWidth', 0.5, 'DisplayName', int2str(n))
end
end
title('Fourier Series')
legend show

end


Fourier series for a pulse signal

• Have you attempted to take a much smaller step such as $dt=1/10000$ ? Oct 13, 2021 at 9:01
• Yes I did but it did not solve the issue.
– bfgt
Oct 13, 2021 at 9:05
• $n$ is always an integer in the function: $n$ successively takes the integer values from 1 to N but is never a list per se.
– bfgt
Oct 13, 2021 at 10:30
• Sorry, I understand now: it is $t$ which is a list... Oct 13, 2021 at 10:33
• You are simply off by a constant factor. ($2$ ?)
– user974557
Oct 13, 2021 at 17:32

The approach I propose as a reference is to take

$$t=(0:dt:1)*T\tag{1}$$

$$t = (-1+dt:dt:1)*T/2;\tag{2}$$

is perfectly valid under the condition that the "differential element" $$dt$$ (= $$\Delta t$$) is divided by 2 (see below the line with comments %%%%%%%%%).

Why that ? Because in case (2), one uses twice as many values of $$f$$ than in case (1), therefore it has to be compensated by taking a differential element divided by 2 !

Please note that factors $$(2/T)$$ are to be kept in front of formulas for $$A(n)$$ and $$B(n)$$.

function m;
clear all;close all;hold on;
dt = 1/1000;
T = 1; % signal period
% Approach (1):
% t = (0:dt:1)*T;
% f=2*sin(2*pi*t)+0.5*cos(2*pi*5*t)+3*cos(2*pi*3*t);
% plot(t,f,'c','linewidth',3);
% "Your" approach (2)... perfectly valid:
t=(-1+dt:dt:1)*T/2;
dt=dt/2; % <- %%%%%%%%%%%%... under this condition
% (that will be applied later on)
f =2*sin(2*pi*t)+0.5*cos(2*pi*5*t)+3*cos(2*pi*3*t);
plot(t,f,'c','linewidth',3)
plot(t,f,'c','linewidth',3)
N = 5;
get_fourier_trigo(f, T, t, dt, N, 1);
%
function get_fourier_trigo(f, T, t, dt, N, k_plot)
s='rgbcmyrgbcmy';
A0 = (2/T)*sum(f)*dt;
fFS = A0/2;
u=t*2*pi/T;
for n=1:N
A(n) = (2/T)*sum(f.*cos(n*u))*dt
B(n) = (2/T)*sum(f.*sin(n*u))*dt
fFS = fFS + A(n)*cos(n*u) + B(n)*sin(n*u);
c=s(n); % plotting color
plot(t, fFS, '-','color',c, 'LineWidth', 0.5, 'DisplayName', int2str(n));
end
title('Fourier Series')
legend show

• So if I understand correctly, you need to integrate on the interval $[0; T]$ and not $[-T/2; T/2]$ ? I thought I could integrate on any interval of length T...
– bfgt
Oct 13, 2021 at 13:42
• Rectification: indeed, you are perfectly right when saying that we should obtain the same result using any interval of length $T$, be it $[0,T]$ or $[T/2,T/2]$ ; it has taken me some time before realizing that the issue is in the second case to take a "dt" divided by two because in this case we take twice as many values of $f$. See my updated answer. Oct 13, 2021 at 17:26
• Oh yes you are perfectly right, the problem was in the creation of the time vector. Thank you for spotting this out!
– bfgt
Oct 13, 2021 at 18:50
• I could have said it otherwise by taking the same $dt$ but with $(-\frac12+dt:dt:\frac12)$. Oct 13, 2021 at 19:03
• Unfortunately I still have some issues. I now try to obtain the Fourier series of a pulse. % Define function dt = 1/1000; % sampling time T = 1; % periodicity t = (0:dt:1-dt)*T; % split interval in quarters n = length(t); nquart = floor(n/4); % Pulse function f = 0*t; f(nquart:3*nquart) = 1; I have alse replaced the sums by numerical integration using trapz but it does not help.
– bfgt
Oct 14, 2021 at 6:09