0
$\begingroup$

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

$\endgroup$
5
  • $\begingroup$ Have you attempted to take a much smaller step such as $dt=1/10000$ ? $\endgroup$
    – Jean Marie
    Oct 13, 2021 at 9:01
  • $\begingroup$ Yes I did but it did not solve the issue. $\endgroup$
    – bfgt
    Oct 13, 2021 at 9:05
  • $\begingroup$ $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. $\endgroup$
    – bfgt
    Oct 13, 2021 at 10:30
  • $\begingroup$ Sorry, I understand now: it is $t$ which is a list... $\endgroup$
    – Jean Marie
    Oct 13, 2021 at 10:33
  • $\begingroup$ You are simply off by a constant factor. ($2$ ?) $\endgroup$
    – user974557
    Oct 13, 2021 at 17:32

1 Answer 1

0
$\begingroup$

The approach I propose as a reference is to take

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

But your approach with

$$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
$\endgroup$
12
  • $\begingroup$ 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... $\endgroup$
    – bfgt
    Oct 13, 2021 at 13:42
  • $\begingroup$ 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. $\endgroup$
    – Jean Marie
    Oct 13, 2021 at 17:26
  • $\begingroup$ Oh yes you are perfectly right, the problem was in the creation of the time vector. Thank you for spotting this out! $\endgroup$
    – bfgt
    Oct 13, 2021 at 18:50
  • $\begingroup$ I could have said it otherwise by taking the same $dt$ but with $(-\frac12+dt:dt:\frac12)$. $\endgroup$
    – Jean Marie
    Oct 13, 2021 at 19:03
  • $\begingroup$ 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. $\endgroup$
    – bfgt
    Oct 14, 2021 at 6:09

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .