Traitement du Signal
Séries de Fourier

Séries de Fourier

Constat

Tout signal Périodique est la somme de sinusoïdes de différentes fréquences.

Ex: Signal carré

à mesure que j’ajoute des sinusoïdes, le signal carré se forme.
Je note l’amplitude de chacune ces sinusoides dans le plan des fréquences.

sf_carre.svg

Définition

y(t) Fonction T-périodique

y(t)=a0+∑n=1+∞[ancos⁡nωt+bnsin⁡nωt] \boxed{ y(t)=a_0+\sum_{n=1}^{+\infty}[a_n\cos n\omega t + b_n \sin n\omega t] }

a0 valeur moyenne:a0=1T∫0Ty(t)dt a_0\ valeur\ moyenne : a_0=\frac{1}{T}\int_{0}^{T}y(t)dt

Toute fonction T-périodique résulte d’une somme infinie de cosinus et sinus.

{an=2T∫0Ty(t)cos⁡(nωt)dtbn=2T∫0Ty(t)sin⁡(nωt)dt \left\{\begin{matrix} a_n=\frac{2}{T}\int_{0}^{T}y(t)\cos( n\omega t) dt \\ b_n=\frac{2}{T}\int_{0}^{T}y(t)\sin( n\omega t) dt \end{matrix}\right.

sf3d.svg

Démonstration (facultatif)

y(t).cos⁡(kωt)=a02.cos(kωt)+cos⁡(kωt)∑n=1N(ancos⁡(nωt)+bnsin⁡(nt)) y(t).\cos(k\omega t)=\frac{a_0}{2}.cos(k\omega t)+\cos(k\omega t)\sum_{n=1}^{N}(a_n \cos(n\omega t)+b_n \sin(nt))
y(t).cos⁡(kωt)=a02.cos(kωt)+∑n=1N(ancos⁡(nωt)cos⁡(kωt)+bnsin⁡(nωt)cos⁡(kωt)) y(t).\cos(k\omega t)=\frac{a_0}{2}.cos(k \omega t)+\sum_{n=1}^{N}(a_n \cos(n \omega t) \cos(k \omega t)+b_n \sin(n \omega t)\cos(k \omega t))
∫−T+Ty(t)cos(kωt)dt=∫−T+Ta02cos⁡(kωt)dt+∫−T+T∑n=1N(ancos⁡(nωt)cos⁡(kωt)+bnsin⁡(nωt)cos⁡(kωt))dt \int_{-T}^{+T} y(t)cos(k\omega t)dt = \int_{-T}^{+T} \frac{a_0}{2} \cos(k \omega t)dt + \int_{-T}^{+T} \sum_{n=1}^{N}(a_n \cos(n \omega t) \cos(k \omega t)+b_n \sin(n \omega t)\cos(k \omega t))dt
∫−T+Ty(t)cos(kωt)dt=a02∫−T+Tcos⁡(kωt)dt+∑n=1Nan∫−TTcos⁡(nωt)cos⁡(kωt)dt+∑n=1Nbn∫−TTsin⁡(nωt)cos⁡(kωt)dt \int_{-T}^{+T} y(t)cos(k\omega t)dt = \frac{a_0}{2} \int_{-T}^{+T} \cos(k \omega t)dt + \sum_{n=1}^{N} a_n \int_{-T}^{T} \cos(n \omega t) \cos(k \omega t) dt + \sum_{n=1}^{N} b_n \int_{-T}^{T} \sin(n \omega t) \cos(k \omega t)dt
parmi les n, il existe un et un seul n égal à k tel que :

∑n=1Nan∫−TTcos⁡(nωt)cos⁡(kωt)dt=∑n=1,n≠kNan∫−TTcos⁡(nωt)cos⁡(kωt)dt+ak∫−TTcos⁡2(kωt)dt=akπ \sum_{n=1}^{N} a_n \int_{-T}^{T} \cos(n \omega t) \cos(k \omega t) dt = \sum_{n=1, n\neq k}^{N} a_n \int_{-T}^{T} \cos(n \omega t) \cos(k \omega t) dt + a_k \int_{-T}^{T} \cos^2(k\omega t) dt = a_k \pi
ak=1π∫−TTy(t)cos⁡(kωt)dt a_k = \frac{1}{\pi}\int_{-T}^{T}y(t)\cos(k \omega t) dt

Forme Complexe

y(t)=a0+∑n=1+∞[ancos⁡nωt+bnsin⁡nωt] y(t)=a_0+\sum_{n=1}^{+\infty} \left [a_n\cos n\omega t + b_n \sin n\omega t \right ]
y(t)=a0+∑n=1+∞[anejnωt+e−jnωt2+bnejnωt−e−jnωt2j] y(t)=a_0+\sum_{n=1}^{+\infty} \left [a_n \frac{e^{jn\omega t}+e^{-jn\omega t}}{2} + b_n \frac{e^{jn\omega t}-e^{-jn\omega t}}{2j} \right ]
y(t)=a0+∑n=1+∞[an−jbn2ejnωt+an+jbn2e−jnωt] y(t)=a_0+\sum_{n=1}^{+\infty} \left [\frac{a_n-jb_n}{2}e^{jn\omega t} + \frac{a_n+jb_n}{2}e^{-jn\omega t} \right ]

On pose : cn=an−jbn2 c_n= \frac{a_n-jb_n}{2} , c−n=an+jbn2 c_{-n}= \frac{a_n+jb_n}{2} , c0=a0 c_0=a_0

y(t)=a0+∑n=1+∞cnejnωt+∑−n=1+∞c−ne−jnωt y(t)=a_0+\sum_{n=1}^{+\infty} c_n e^{jn\omega t} + \sum_{-n=1}^{+\infty} c_{-n} e^{-jn\omega t}
y(t)=a0+∑n=1+∞cnejnωt+∑n=−1−∞cnejnωt y(t)=a_0+\sum_{n=1}^{+\infty} c_n e^{jn\omega t} + \sum_{n=-1}^{-\infty} c_{n} e^{jn\omega t}

y(t)=∑n=−∞+∞cnejnωt \boxed{ y(t)=\sum_{n=-\infty}^{+\infty}c_n e^{jn\omega t} }

cn=1T∫−T2T2y(t)e−jnωtdt \boxed{ c_n=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}y(t)e^{-jn\omega t} dt }

La forme complexe permet d’estimer en une seule expression la contribution sinusoïdale et cosinusoïdale.

La contrepartie de cette simplification mathématique est désormais la nécessité de considérer des ‘fréquences négatives’ pour le calcul.

Exemples

Signal carré y(t) d’amplitude U

y(t)=4Uπ∑n=0+∞sin⁡[(2n+1)ωt]2n+1 y(t)=\frac{4U}{\pi}\sum_{n=0}^{+\infty}\frac{\sin[(2n+1)\omega t]}{2n+1}
y(t)=4Uπ[sin⁡(ωt)+13sin⁡(3ωt)+15sin(5ωt)+...+sin[(2n+1)ωt]2n+1] y(t)=\frac{4U}{\pi}[\sin(\omega t)+\frac{1}{3}\sin(3\omega t)+\frac{1}{5}sin(5\omega t)+...+\frac{sin[(2n+1)\omega t]}{2n+1}]

Signal triangulaire y(t) d’amplitude U

y(t)=−8Uπ2∑n=0+∞cos⁡[(2n+1)ωt](2n+1)2 y(t)=-\frac{8U}{\pi^2}\sum_{n=0}^{+\infty}\frac{\cos[(2n+1)\omega t]}{(2n+1)^2}
y(t)=−8Uπ2[cos⁡ωt+19cos⁡(3ωt)+125cos⁡(5ωt)+...+cos[(2n+1)ωt](2n+1)2] y(t)=-\frac{8U}{\pi^2}[\cos\omega t + \frac{1}{9}\cos(3\omega t)+\frac{1}{25}\cos(5\omega t) +...+ \frac{cos[(2n+1)\omega t]}{(2n+1)^2} ]

sf_triangle.svg

Interprétation Mathématique

Espace Euclidien et Produit scalaire

vecteur.svg

Je considère une base géométrique orthogonale (orthogonale # indépendant, il me faut ici 2 coordonnées pour exprimer un point)

V→=2.X→+5.Y→ \overrightarrow{V} = 2.\overrightarrow{X} + 5.\overrightarrow{Y}

Si je veux connaître la coordonnée en X, j’applique le produit scalaire :

V→[25] \overrightarrow{V} \begin{bmatrix} 2\\ 5 \end{bmatrix}

X→[10] \overrightarrow{X} \begin{bmatrix} 1\\ 0 \end{bmatrix}

⟨V,X⟩=V→.X→=[20] \langle V,X \rangle = \overrightarrow{V}.\overrightarrow{X} = \begin{bmatrix} 2\\ 0 \end{bmatrix}

Ainsi, mon vecteur existe dans une base donnée.
Je peux le décrire à travers une somme faisant apparaitre les composantes de cette base ( X et Y).
Pour calculer une composante en particulier, j’utilise le produit scalaire.

Retour aux séries de Fourier : Produit scalaire hermitien

Ce concept de base peut s’étendre au domaine spectral, l’objectif étant alors d’exciter chacune des composantes sinusoïdales de notre signal.

( REM : on utilise le concept de produit scalaire dans une base hermitienne : <f,g>=∫f(t)g(t)ˉdt <f,g>= \int f(t) \bar{g(t)} dt )

Considérons un signal y(t)=cos⁡(2πFt) y(t)= \cos(2 \pi Ft) .
Afin de représenter ce signal dans un espace fréquentiel, je vais comparer y(t) avec des signaux sinusoidaux.
J’utilise la forme exponentielle pour détecter la contribution sinusoidale et cosinusoidale.

Supposons dans un premier temps que mon signal e−2πFat e^{-2 \pi F_at } ait une fréquence Fa=FF_a = F
En décomposant y(t) en une somme d’exponentielles, et en appliquant le produit scalaire ∫−T/2T/2y(t)e−2πFatdt \int_{-T/2}^{T/2}y(t)e^{-2\pi F_at} dt
J’obtiens un résultat non nul.
Cette intégration, ou produit scalaire, permet d’estimer le degré de ressemblance ( ou corrélation ) entre mon signal y(t) et un signal sinusoidal.
S’agissant de 2 sinusoïdes de même fréquence, ce résultat est non nul.

serie_fourier_tune.svg

ANALYSE SPECTRALE :

Je peux “ranger” ce résultat dans une base fréquentielle :

spectre_cos.svg

Si je fais le même calcul pour une fréquence Fa≠FF _a \neq F ( ou -F), le résultat de mon intégrale sera nul.

serie_fourier_tune_2.svg

ANALYSE SPECTRALE :

spectre_cos_2.svg

Script Scilab et Matlab/Octave pour générer les figures ci dessus

SCILAB
serie_fourier.sce 
//==============================================================================
//              SQUARE SIGNAL
//==============================================================================
N = 1000;       
t = [0:N-1]
F=10            
Fe=1000 
Te=1/Fe
F_norm=F/Fe

Nval = [1, 3, 5, 51];
figure(1);

for in = 1:length(Nval),
    v = 0;
    subplot(4,2,(2*in-1));               // TIME
            for n = 1:2:Nval(in),
                addto = 4*sin(2*n*%pi*F_norm*t)/(n*%pi);
                v     = v + addto;
                plot(t/(N),addto,'b--','LineWidth',1);
            end
            plot(t/(N),v,'r-','LineWidth',2);    
            a=gca()
            t_min=0
            t_max=2/F
            y_min=-1.5
            y_max=1.5
            a.data_bounds=[t_min,t_max,y_min,y_max];  
            xlabel('t (s)','FontSize',1); ylabel('FS Approx to y(t)','FontSize',1);
            title(['Fourier Series Approx with',string(Nval(in)),' Terms'],'FontSize',2);
            
 
    subplot(4,2,(2*in));                // SPECTRUM

            X = abs(fft(v));
            X = fftshift(X);
            vect_Freq = [-N/2:N/2-1]/N;
            plot2d2(Fe*vect_Freq,X/(N/2),rect=[-200,0,200,1.5])       // [xmin,ymin,xmax,ymax]  
            p = get("hdl");
            p.children.thickness = 2;
            p.children.foreground = 3; 
            xlabel('f (Hz)','FontSize',1);   
            title(['Spectrum Analysis with',string(Nval(in)),' Terms'],'FontSize',2);            
end

//==============================================================================
//              TRIANGLE SIGNAL
//==============================================================================
N = 1000;       
t = [0:N-1]
F=10            
Fe=1000 
Te=1/Fe
F_norm=F/Fe

Nval = [1, 3, 5, 51]; 
figure(2);

for in = 1:length(Nval),
    v = 0;
    subplot(4,2,(2*in-1));               // TIME
            for n = 1:2:Nval(in),
                addto = 8*cos(2*n*%pi*F_norm*t)/(n*n*%pi*%pi);
                v     = v + addto;
                plot(t/(N),addto,'b--','LineWidth',1);
            end
            plot(t/(N),v,'r-','LineWidth',2);    
            a=gca()
            t_min=0
            t_max=2/F
            y_min=-1
            y_max=1
            a.data_bounds=[t_min,t_max,y_min,y_max];  
            xlabel('t (s)','FontSize',1); ylabel('FS Approx to y(t)','FontSize',1);
            title(['Fourier Series Approx with',string(Nval(in)),' Terms'],'FontSize',2);
            
 
    subplot(4,2,(2*in));                // SPECTRUM
            X = abs(fft(v));
            X = fftshift(X);
            vect_Freq = [-N/2:N/2-1]/N;
            plot2d2(Fe*vect_Freq,X/(N/2),rect=[-200,0,200,1])       // [xmin,ymin,xmax,ymax]  
            p = get("hdl");
            p.children.thickness = 2;
            p.children.foreground = 3; 
            xlabel('f (Hz)','FontSize',1);   
            title(['Spectrum Analysis with',string(Nval(in)),' Terms'],'FontSize',1);            
end
MATLAB / OCTAVE
serie_fourier_tune_freq.m 
Delta_t=1e-3          % Pas de Calcul
N=1/Delta_t           % Nombre de points
A=1                   % Amplitude
F=10                  % Frequence
T=1/F                 % Periode
Nb_per=3              % Nombre de Periodes
t=[-N*Nb_per*T/2:(N-1)*Nb_per*T/2]/N;

F_analyze = 3*F
%===============================================================================
% 1 %                    SAME FREQUENCY ANALYSIS
%===============================================================================
re=cos(-2*pi*F_analyze*t);
im=sin(-2*pi*F_analyze*t);
y=cos(2*pi*F*t);

%...............................................................................
figure()
subplot(2,2,1)

plot3(t,re,im,"b", "linewidth", 2)
hold on
plot3(t,y,"r","linewidth", 5)
hold on

xlabel("t(s)")
ylabel('real : cos\omega t')
zlabel('imag : sin\omega t')

str = {strcat('y(t)=cos(2\piFt) with F = ' , num2str(F),'Hz, analyzed with e^{-j2\pi F_at} with F_a=', num2str(F_analyze),'Hz')} 
title(str)
% str = {strcat('e^{-j2\pi Fa t} with Fa=' , num2str(F_analyze),'Hz')} 
str = {strcat('e^{-j2\pi Fa t} with Fa=' , num2str(F_analyze),'Hz')} 
legend(str, 'y(t)=cos(2\piFt)',"location", "northeast")
grid
%===============================================================================
% 2 %                COS BECOMES A SUM OF EXP
%===============================================================================
re=cos(-2*pi*F_analyze*t);
im=sin(-2*pi*F_analyze*t);
y_1_re=0.5*cos(2*pi*F*t);
y_1_im=0.5*sin(2*pi*F*t);
y_2_re=0.5*cos(-2*pi*F*t);
y_2_im=0.5*sin(-2*pi*F*t);

%...............................................................................
subplot(2,2,2)

plot3(t,re,im,"b", "linewidth", 2)
hold on
plot3(t,y_1_re,y_1_im,"m","linewidth", 5)
hold on
plot3(t,y_2_re,y_2_im,"g","linewidth", 5)
hold on

xlabel("t(s)")
ylabel('real : cos\omega t')
zlabel('imag : sin\omega t')
str = {strcat('y(t)=(e^{j\omega t}+e^{-j\omega t})/2 with F = ' , num2str(F),'Hz')} 
title(str)
str = {strcat('e^{-j2\pi Fa t} with Fa=' , num2str(F_analyze),'Hz')} 
legend (str,'y(t)=e^{j2\pi Ft}','y(t)=e^{-j2\pi Ft}',"location", "northeast");
grid
%===============================================================================
% 3 %                PRODUCT FOR EACH t : exp(+jwt)
%===============================================================================
re=cos(-2*pi*F_analyze*t);
im=sin(-2*pi*F_analyze*t);
y_1_re=0.5*cos(2*pi*F*t);
y_1_im=0.5*sin(2*pi*F*t);

product=0.5*exp(j*2*pi*F*t).*exp(-j*2*pi*F_analyze*t);

nd=100
product_area = product.*[zeros(1,(length(t)-nd)/2), ones(1, nd), zeros(1, (length(t)-nd)/2)] ; 

%...............................................................................
subplot(2,2,3)

plot3(t,re,im,"b", "linewidth", 2)
hold on
plot3(t,y_1_re,y_1_im,"m","linewidth", 5)
hold on
plot3(t,product,"y","linewidth", 5)
area(t,product_area)

xlabel("t(s)")
ylabel('real : cos\omega t')
zlabel('imag : sin\omega t')
str = {strcat('y_1(t)=e^{j\omega t}/2 with F = ' , num2str(F),'Hz')} 
title(str)
str = {strcat('e^{-j2\pi Fa t} with Fa=' , num2str(F_analyze),'Hz')} 
legend (str,'y_1(t)=0.5.e^{j2\pi Ft}','y_1(t).e^{-j2\pi Ft}','\int_{-T/2}^{T/2}y_1(t).e^{-j2\pi Ft}dt',"location", "northeast");
grid
%===============================================================================
% 4 %                PRODUCT FOR EACH t : exp(+jwt)
%===============================================================================
re=cos(-2*pi*F_analyze*t);
im=sin(-2*pi*F_analyze*t);
y_2_re=0.5*cos(-2*pi*F*t);
y_2_im=0.5*sin(-2*pi*F*t);

product=0.5*exp(-j*2*pi*F*t).*exp(-j*2*pi*F_analyze*t);

nd=100
product_area = product.*[zeros(1,(length(t)-nd)/2), ones(1, nd), zeros(1, (length(t)-nd)/2)] ; 

%...............................................................................
subplot(2,2,4)

plot3(t,re,im,"b", "linewidth", 2)
hold on
plot3(t,y_2_re,y_2_im,"g","linewidth", 5)
hold on
plot3(t,product,"y","linewidth", 5)
area(t,product_area)

xlabel("t(s)")
ylabel('real : cos\omega t')
zlabel('imag : sin\omega t')
str = {strcat('y_2(t)=e^{-j\omega t}/2 with F = ' , num2str(F),'Hz')} 
title(str)
str = {strcat('e^{-j2\pi Fa t} with Fa=' , num2str(F_analyze),'Hz')} 
legend (str,'y_2(t)=0.5.e^{-j2\pi Ft}','y_2(t).e^{-j2\pi Ft}','\int_{-T/2}^{T/2}y_2(t).e^{-j2\pi Ft}dt',"location", "northeast");
grid

%===============================================================================
%                 FREQUENCY ANALYSIS
%===============================================================================
figure()
freq = [0:1:100];
amplitudes=[zeros(1,length(freq))];

Delta_t=1e-3          % Pas de Calcul
N=1/Delta_t           % Nombre de points
A=1                   % Amplitude
F=10                  % Frequence
T=1/F                 % Periode
Nb_per=1              % Nombre de Periodes
t=[-N*Nb_per*T/2:(N-1)*Nb_per*T/2]/N;

y=cos(2.*pi.*F.*t);
TF=exp(-i.*2.*pi.*F_analyze.*t).*y;

acc_TF=0;
for n=1:length(t)
        acc_TF = acc_TF+TF(n).*Delta_t;
end
acc_TF=acc_TF/T;

amplitudes(F_analyze)=abs(acc_TF);
%-------------------------------------------------------------------------
stem(freq,amplitudes)
axis([0 max(freq) 0 1])
hold on
%===============================================================================
Produit de ConvolutionTransformée de Fourier