Bonjour à toutes et à tous,

Je viens ici pour essayer de comprendre la différence en transformé de Fourier discrète et l'algorithme fft. J'ai réalisé un script qui me calcule la transformé de Fourier d'une fonction en réalisant une simple intégration. Par soucis de rapidité j'aimerais utiliser l'algorithme fft mais je n'arrive pas à comprendre comment il fonctionne et surtout je n'arrive pas à obtenir les résultats que j'avais avec mon ancien script (résultats validés).

Par exemple, si je veux réalise la fft d'une fonction porte, je n'arrive pas à retrouvé un sinus cardinal. Aussi quand je réalise une transformée de Fourier, je suis sensée obtenir une fonction dans l'espace des nombres d'onde qui dépend de s, je ne comprends pas comment retrouver ce s avec les fft.

Merci par avance

-
Edité par Lucie2356 7 novembre 2020 à 11:45:45

Bonjour Lucie,

l'algorithme FFT permet de calculée la transformée de Fourier discrète en \(\mathcal O(N\cdot\log(N))\) opérations contre \(\mathcal O(N^2)\) pour une implémentation naïve. Tu peux voir la transformée discrète comme un calcul de la transformée continue par la méthode des rectangles. En quel langage est écrit ton script; tu pourrais peut-être nous le montrer ?

Bonjour Nozio,

Merci de ta réponse. Si je comprends bien l'algorithme FFT est une sorte de réarrangement des sommes de la transformée discretes?

alors j'ai tenté de calculer la transformée fourrier de ma fonction de trois manières différentes (j'ai pris une fonction porte d'amplitude 0.5) :

Transformé continue :

def fourier(a, b, dx):
    if a > b:
        a,b = b,a
    n = int((b - a) / dx)
    t = np.linspace(-20,20,n//2)
    f_i = np.empty_like(t,dtype=np.complex64)
    for k in range(0,t.shape[0]):
        s = 0.0
        x = a
        for i in range(n):
            f_i[k] = 0.5*np.exp(-complex(0,1)*x*t[k])
            s += f_i[k]
            x += dx
        f_i[k] = s * dx
    
    return f_i,t

Transformé discrète :

def fourier_discret(N):
    n = np.linspace(0,20,N)
    k = np.linspace(-20,20,200)
    X = np.empty_like(k,dtype=np.complex64)
    for i in range(0,k.size):
        s = 0.0
        for j in range(0,N-1):
            s += 0.5*np.exp(-complex(0,1)*k[i]*n[j]/N)
        X[i] = s
    
    return X,k

Avec les fft :

y = 0.5*np.ones_like(x)
mask_neg = x<-1
mask_pos= x>1
y[mask_neg]=0
y[mask_pos]=0

yf = np.fft.fft(y)
xf = np.fft.fftfreq(x.size,x[1]-x[0])

J'ai du mal avec les deux dernières formulations, je ne comprends pas pourquoi je n'ai pas les mêmes courbes finales entre la méthode qui consiste à calculé une intégrale et celle dite discrète. Quand j'augmente ma valeur de N (algo 2), l'amplitude de mon signal augmente alors que ça ne devrait pas avoir lieu...

Salut,

je vois plusieurs choses :

• pour ton algo 1 tu calcules la FFT sur \((-\infty,\infty)\) ? Il faut bien voir que appliquer une FFT suppose que ton signal est périodique. Par exemple, si tu as un créneau sur \([-1/2,1/2]\) et que tu calcules sa FFT sur \([-1,1]\), tu supposes qu'il est de période 2. Si tu veux calculer la version "continue", tu dois calculer la version continue de la transformée de Fourier périodique (je me trompe peut-être mais je pense qu'il te manque une division par \(b-a\) dans ton algo 1.
• Ensuite, point relatif à l'algo 1 et au 2 par rapport au 3 : tu as bien utilisé la définition standard de la transformée de Fourier (celle en \(\omega \Leftrightarrow k\)) sauf qu'un coup d'oeil ici : https://docs.scipy.org/doc/scipy/reference/generated/scipy.fft.fft.html indique que le calcul est fait en fréquence donc pour \(f = \dfrac{\omega}{2\pi}\) comme tu l'as bien vu puisque tu utilises fftfreq.
• Dans ton algo 2 (comme pour le 3), il te manque une normalisation par \(1/N\). Pour t'en convaincre, tu peux réécrire ce que ça signifie d'intégrer l'intégrale de Fourier périodique par méthode des rectangles. Tu n'as pas ce problème dans l'algo 1 parce que cette division est cachée dans le \(dx\).

Pour ce qui est de l'algo FFT original, il s'agit de décomposer la somme en termes pairs et impairs (en terme d'indice de la somme), factoriser une exponentielle complexe pour obtenir deux termes pairs (\(E_k,O_k\) sur cette page), puis procéder récursivement sur les termes obtenus jusqu'à ce que ce ne soit plus possible (d'où la nécessité d'avoir un nombre de termes en puissances de 2).

Bonjour Nozio,

Merci de ta réponse complète, ça m'aide beaucoup mais j'ai encore deux trois petites questions.

• Je pensais qu'une transformée de Fourier devait toujours être calculé sur $(-\infty,\infty)$, je ne comprends donc pas très bien par version continue de la transformée de Fourier périodique ? (Oui, je pense qu'il manque la division par la période, je n'y avais pas pensé ;))
• Si je comprends bien le second point, cela veut dire que je suis obligée de travailler en fréquentiel avec les FFT. Dans ce cas, comment faire pour calculé la transformée de Fourier d'un simple sin(x) et récupérer l'abscisse xf de mon signal transformé sans passer par les fréquences ?
• Merci pour le troisième point !

T = 2*np.pi # periode du signal
N = 20      #nbperiode
x = np.linspace(0,T,N)
y = np.exp(-x**2/100)*np.sin(x)

yf = np.fft.fftshift(y)
xf = np.fft.fftshift(np.fft.fftfreq(x.size, T)) #fréquences correspondante au signal yf

plt.plot(xf,1/N*yf)

-
Edité par Lucie2356 10 novembre 2020 à 9:54:16

La transformée de Fourier périodique (décomposition en série de Fourier) est simplement (voir ici)

\(\dfrac{1}{T}\int_{-T/2}^{T/2}{f(t) \,e^{-i \xi t} dt}\)

Je ne comprends pas trop ton exponentielle devant ton sinus, mais si tu veux simplement regarder pour le sinus de fréquence 1 sur \([0,1]\), il suffit de faire

import numpy as np
import matplotlib.pyplot as plt

T  = 1
N  = 20
dt = T/N

x = np.arange(0,T,dt)
y = np.sin(2*np.pi*x)

yf = np.fft.fftshift(np.fft.fft(y))/N
xf = np.fft.fftshift(np.fft.fftfreq(N,dt))

plt.figure()
plt.plot(xf,2*np.abs(yf))
plt.title("Spectre en amplitude du sinus")
plt.show()

Evite le linspace car tu dois exclure le point  \(x = 1\) (rappel: on suppose que le signal est périodique, or \(x=1\) est le premier point de la seconde période ! Attention aussi à l'utilisation de fftfreq (c'est dt et pas T). Enfin, dernier piège : il faut multiplier le spectre par 2. La raison est que la transformée discrète produit des fréquences "négatives" qui correspondent en fait aux fréquences opposées positives (en gros le spectre est symétrique) et l'amplitude est répartie à part en égale entre la fréquence positive et la fréquence négative correspondante. Donc si tu as un sinus d'amplitude 1, tu obtiendras 1/2 et 1/2 en les points -1 et 1. Pour l'affichage, j'ai tout remultiplié par 2.

-
Edité par Nozio 10 novembre 2020 à 13:07:20

D'accord, je crois que je comprends. Je viens d'essayer d'utiliser mon algo pour essayer de voir si j'obtenais la même chose qu'avec ton code mais quelque chose cloche... Je ne vois pas quoi.

def fourier_discret():
    T = 2*np.pi
    N = 20
    dt = T/N
    
    n = np.arange(0,T,dt)
    k = np.arange(-200,200,dt)
    
    X = np.empty_like(k,dtype=np.complex64)
    for i in range(0,k.size):
        s = 0.0
        for j in range(0,N-1):
            s += np.sin(n[j])*np.exp(-complex(0,1)*k[i]*n[j]/N)
            
        X[i] = s/N
    
    return X,k

yf,xf = fourier_discret()
plt.plot(xf,(yf.real)*2)

Je me aussi s'il est possible d'augmenter le nombre de point de la fft sans jouer sur les fréquences. Le signal me paraît imprécis. Je me pose aussi une autre question sur ce N. Si j'ai bien compris il dois être une puissance de deux pour que la fft fonctionne correctement mais lorsque je met un N = 19 par exemple, cela semble quand même fonctionner. Comment est possible si le signal doit être décomposé en somme paire ?

Dernière question : si je veux calculer le TF de \( F(-s) = \int_{-\infty}^{\infty} f(t)e^{jst}\)

J'avais pensé passer par les transformées de Fourier inverse ce qui me ferais sortie un facteur N si je suis ce lien. Est-ce viable comme méthode ?

 Edit : je m'excuse pour la formule, je ne sais pas comment utiliser Latex sur ce forum




-
Edité par Lucie2356 11 novembre 2020 à 9:45:16

Pour Latex, c'est

\( ... \)

Ensuite, la contrainte \(N = 2^n\) permet juste d'être optimal pour la calcul de la FFT, mais c'est uniquement un aspect algorithmique qui t'es caché par l'interface que tu appelles. La transformée discrète est définie quelque soit N. Tu peux augmenter le nombre de fréquences en augmentant ta fréquence d'échantillonnage (en diminuant \(dt\)).

Pour ce qui est de ton code, si tu veux répliquer le comportement de la FFT, il doit bien sûr calculer ce que l'interface Python calcule ! D'après la doc, et pour \(k \in [\![0,N-1]\!]\), on récupère

y[k] = np.sum(x * np.exp(-2j * np.pi * k * np.arange(n)/n))

or toi tu donnes explicitement les fréquences. Dans un premier temps, contente-toi juste de prendre des k entiers entre 0 et N-1 (k[i] = i). Tu devrais alors retrouver le résultat voulu.

Merci ! Les deux méthodes fonctionnent et me donnent les mêmes résultats ! Concernant les FFT d'ordre n, il s'agit bien d'appliquer la fft selon les deux axes ?

Normalement oui

Je crois que ça fonctionne mais j'ai un dernier petit point dont j'aimerais être sûre. Si je fais l'analogie avec le 1D, je dois normaliser mon signal par 1/MN. Au niveau de son amplitude je serais tenté de tout multiplier par 4 mais j'en suis moins certaine.Ce qui me donnerai quelque chose comme ça :

T = 2*np.pi
N = 20
M = 20
dx = T/N
dy = T/M
x = np.arange(0,T,dx)
y = np.arange(0,T,dy)

y1 = np.sin(x)
y2 = np.sin(y)


Y = y1[:,None]*y2[None,:]

yfft = np.fft.fftshift(np.fft.fft2(Y)/(M*N))

xf = np.fft.fftshift(np.fft.fftfreq(M,dx))
zf = np.fft.fftshift(np.fft.fftfreq(N,dy)) 

%% Post-traitement
X, Y = np.meshgrid(xf, zf)
fig = plt.figure()
my_col = cm.jet(yfft.real/np.amax(yfft.real))
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y,4*np.abs(yfft), rstride=1, cstride=1, facecolors = my_col,
        linewidth=0, antialiased=False)

Au niveau de la période de mon signal, celle ci sera bien toujours égale à 2 pi vu que j'ai un sinus ?

-
Edité par Lucie2356 11 novembre 2020 à 10:11:06

J'ai envie de dire oui, sinon ton spectre n'a pas la "bonne" amplitude (\(0.25\) au lieu de \(1\) quand on multiplie par 4). Ta période est bien \(2\pi\), ta fréquence est \(1/(2\pi)\).

C'est encore moi !! J'ai toujours des soucis avec mes transformées de Fourier. J'ai tourné le problème dans tous les sens sans résultats.

def nextpow2(A):     
    p = 0
    while 2**p < np.abs(A):  
        p += 1 
    return p

et mon code, il tourne mais j'essaie d'obtenir des résultats issu d'un article mais je n'y arrive pas. Je suis presque sûre que cela viens de la fft. Le restant des calculs étant plutôt classique et je pense pas qu'il y est des soucis là dessus.

f = 1000
L = 500

nfft = 2**nextpow2(L)                      #nombres de points de la fft

freq = f/2 * np.linspace(0,1,nfft//2+1)

# longueur des signaux
M = 50
N = 70

# nombres de points des fft
nfftx = 2**nextpow2(M)
nffty = 2**nextpow2(N)

# période spatiale
Dx = Lx/(M-1);
Dy = Ly/(N-1);

# vecteurs espace
x = np.arange(0,M)*Dx
y = np.arange(0,N)*Dy

kx = 2*np.pi/(2*Dx)*np.arange(1,nfftx//2+1)/(nfftx//2+1)
ky = 2*np.pi/(2*Dy)*np.arange(1,nffty//2+1)/(nffty//2+1)

dkx = kx[1]-kx[0]
dky = ky[1]-ky[0]
# %% ------------- FFT -------------
Phi_mn = np.sin(x)[:,None]*np.sin(y)[None,:]
PHI_mn2 = np.fft.fft2(Phi_mn)/(M*N)

PHI_mn = (PHI_mn2[0:nfftx//2, 0:nffty//2])

Il faut que tu nous dises ce que ton code doit faire, commente tes différentes étapes ligne par ligne et donne nous le résultat attendu.

Bonjour Nozio,

Je t'ai répondu en message privé mais je peux le faire ici si c'est préférable.