fft

Transformée de Fourier rapide

Syntaxe

Y = fft(X)

Y = fft(X,n)

Y = fft(X,n,dim)

Description

Y = fft(X) calcule la transformée de Fourier discrète (DFT) de X avec un algorithme de transformée de Fourier rapide (FFT). Y est de la même taille que X.

• Si X est un vecteur, fft(X) renvoie la transformée de Fourier du vecteur.

• Si X est une matrice, fft(X) traite les colonnes de X comme des vecteurs et renvoie la transformée de Fourier de chaque colonne.

• Si X est un tableau multidimensionnel, fft(X) traite les valeurs le long de la première dimension du tableau dont la taille n’est pas égale à 1 comme des vecteurs et renvoie la transformée de Fourier de chaque vecteur.

Y = fft(X,n) renvoie la DFT sur n points.

• Si X est un vecteur et que la longueur de X est inférieure à n, X est rempli avec des zéros à droite pour atteindre la longueur n.

• Si X est un vecteur et que la longueur de X est supérieure à n, X est tronqué pour atteindre la longueur n.

• Si X est une matrice, chaque colonne est traitée comme dans le cas d’un vecteur.

• Si X est un tableau multidimensionnel, la première dimension du tableau dont la taille n’est pas égale à 1 est traitée comme dans le cas d’un vecteur.

Y = fft(X,n,dim) renvoie la transformée de Fourier le long de la dimension dim. Par exemple, si X est une matrice, fft(X,n,2) renvoie la transformée de Fourier sur n points de chaque ligne.

Exemples

Signal bruité

Identifiez les composantes fréquentielles d’un signal noyé dans un bruit et déterminez les amplitudes des fréquences de crête avec une transformée de Fourier.

Spécifiez les paramètres d’un signal avec une fréquence d’échantillonnage de 1 kHz et une durée de signal de 1,5 seconde.

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

Formez un signal contenant un offset DC d’amplitude 0,8, une sinusoïde de 50 Hz d’amplitude 0,7 et une sinusoïde de 120 Hz d’amplitude 1.

S = 0.8 + 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

Corrompez le signal avec un bruit aléatoire de moyenne nulle et de variance 4.

X = S + 2*randn(size(t));

Tracez le signal bruité dans le domaine temporel. Les composantes fréquentielles ne sont pas visibles sur le tracé.

plot(1000*t,X)
title("Signal Corrupted with Zero-Mean Random Noise")
xlabel("t (milliseconds)")
ylabel("X(t)")

Calculez la transformée de Fourier du signal.

Y = fft(X);

Comme les transformées de Fourier impliquent des nombres complexes, tracez l’amplitude complexe du spectre fft.

plot(Fs/L*(0:L-1),abs(Y),"LineWidth",3)
title("Complex Magnitude of fft Spectrum")
xlabel("f (Hz)")
ylabel("|fft(X)|")

Le tracé affiche cinq pics de fréquence, y compris le pic à 0 Hz correspondant à l’offset DC. Dans cet exemple, le signal est censé avoir trois pics de fréquence à 0 Hz, 50 Hz et 120 Hz. Ici, la deuxième moitié du tracé est le reflet miroir de la première moitié sans le pic à 0 Hz. En effet, la transformée de Fourier discrète d’un signal dans le domaine temporel est de nature périodique, avec la première moitié de son spectre dans les fréquences positives et la deuxième moitié dans les fréquences négatives, le premier élément étant réservé pour la fréquence nulle.

Pour les signaux réels, le spectre fft est bilatéral, le spectre dans les fréquences positives étant le conjugué complexe du spectre dans les fréquences négatives. Pour afficher le spectre fft dans les fréquences positives et négatives, vous pouvez utiliser fftshift. Pour une longueur paire L, le domaine fréquentiel commence à la valeur négative de la fréquence de Nyquist -Fs/2 et s’étend jusqu’à Fs/2-Fs/L avec un espacement ou une résolution en fréquence de Fs/L.

plot(Fs/L*(-L/2:L/2-1),abs(fftshift(Y)),"LineWidth",3)
title("fft Spectrum in the Positive and Negative Frequencies")
xlabel("f (Hz)")
ylabel("|fft(X)|")

Pour déterminer les amplitudes des trois pics de fréquence, convertissez le spectre fft en Y en spectre d’amplitude unilatéral. Comme la fonction fft inclut un facteur d’échelle L entre les signaux d’origine et transformés, remettez Y à l’échelle en le divisant par L. Prenez l’amplitude complexe du spectre fft. Le spectre d’amplitude bilatéral P2, où le spectre dans les fréquences positives est le conjugué complexe du spectre dans les fréquences négatives, a des amplitudes de crête égales à la moitié de celles du signal dans le domaine temporel. Pour le convertir en spectre unilatéral, prenez la première moitié du spectre bilatéral P2. Multipliez le spectre dans les fréquences positives par 2. Vous n’avez pas besoin de multiplier P1(1) et P1(end) par 2, car ces amplitudes correspondent respectivement à la fréquence nulle et de Nyquist, et elles n’ont pas de paires conjuguées complexes dans les fréquence négatives.

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

Définissez le domaine fréquentiel f pour le spectre unilatéral. Tracez le spectre d’amplitude unilatéral P1. Comme on peut s’y attendre, les amplitudes sont proches de 0,8, 0,7 et 1 mais pas exactement égales à ces valeurs en raison du bruit ajouté. Dans la plupart des cas, les signaux plus longs produisent de meilleures approximations de fréquence.

f = Fs/L*(0:(L/2));
plot(f,P1,"LineWidth",3) 
title("Single-Sided Amplitude Spectrum of X(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")

À présent, prenez la transformée de Fourier du signal non bruité d’origine et récupérez les amplitudes exactes de 0,8, 0,7 et 1,0.

Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

plot(f,P1,"LineWidth",3) 
title("Single-Sided Amplitude Spectrum of S(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Impulsion gaussienne

Convertissez une impulsion gaussienne du domaine temporel vers le domaine fréquentiel.

Spécifiez les paramètres d’un signal avec une fréquence d’échantillonnage de 44,1 kHz et une durée de signal de 1 ms. Créez une impulsion gaussienne avec un écart-type de 0,1 ms.

Fs = 44100;         % Sampling frequency
T = 1/Fs;           % Sampling period
t = -0.5:T:0.5;     % Time vector
L = length(t);      % Signal length

X = 1/(0.4*sqrt(2*pi))*(exp(-t.^2/(2*(0.1*1e-3)^2)));

Tracez l’impulsion dans le domaine temporel.

plot(t,X)
title("Gaussian Pulse in Time Domain")
xlabel("Time (t)")
ylabel("X(t)")
axis([-1e-3 1e-3 0 1.1]) 

La vitesse d’exécution de fft dépend de la longueur de la transformée. L’exécution est beaucoup plus rapide pour les longueurs de transformée ayant seulement de petits facteurs premiers que pour celles qui ont de grands facteurs premiers.

Dans cet exemple, la longueur de signal L est 44 101, qui est un très grand nombre premier. Pour améliorer les performances de fft, identifiez une longueur en entrée correspondant à la prochaine puissance de 2 par rapport à la longueur de signal d’origine. L’appel à la fonction fft avec cette longueur en entrée remplit l’impulsion X avec des zéros à droite jusqu’à atteindre la longueur de transformée spécifiée.

n = 2^nextpow2(L);

Convertissez l’impulsion gaussienne vers le domaine fréquentiel.

Y = fft(X,n);

Définissez le domaine fréquentiel et tracez les fréquences uniques.

f = Fs*(0:(n/2))/n;
P = abs(Y/sqrt(n)).^2;

plot(f,P(1:n/2+1)) 
title("Gaussian Pulse in Frequency Domain")
xlabel("f (Hz)")
ylabel("|P(f)|")

Ondes cosinusoïdales

Comparez des ondes cosinusoïdales dans le domaine temporel et le domaine fréquentiel.

Spécifiez les paramètres d’un signal avec une fréquence d’échantillonnage de 1 kHz et une durée de signal de 1 seconde.

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector

Créez une matrice dont chaque ligne représente une onde cosinusoïdale avec une fréquence mise à l’échelle. Le résultat X est une matrice de 3 x 1 000. La première ligne a une fréquence d’onde de 50, la deuxième ligne a une fréquence d’onde de 150 et la troisième ligne a une fréquence d’onde de 300.

x1 = cos(2*pi*50*t);          % First row wave
x2 = cos(2*pi*150*t);         % Second row wave
x3 = cos(2*pi*300*t);         % Third row wave

X = [x1; x2; x3];

Tracez les 100 premières entrées de chaque ligne de X dans l’ordre sur une même figure et comparez leurs fréquences.

for i = 1:3
    subplot(3,1,i)
    plot(t(1:100),X(i,1:100))
    title("Row " + num2str(i) + " in the Time Domain")
end

Spécifiez l’argument dim pour utiliser fft le long des lignes de X, c’est-à-dire pour chaque signal.

dim = 2;

Calculez la transformée de Fourier des signaux.

Y = fft(X,L,dim);

Calculez le spectre bilatéral et le spectre unilatéral de chaque signal.

P2 = abs(Y/L);
P1 = P2(:,1:L/2+1);
P1(:,2:end-1) = 2*P1(:,2:end-1);

Dans le domaine fréquentiel, tracez le spectre d’amplitude unilatéral pour chaque ligne sur une même figure.

for i=1:3
    subplot(3,1,i)
    plot(0:(Fs/L):(Fs/2-Fs/L),P1(i,1:L/2))
    title("Row " + num2str(i) + " in the Frequency Domain")
end

Phase de sinusoïdes

Créez un signal constitué de deux sinusoïdes dont les fréquences sont de 15 et 40 Hz. La première sinusoïde est une onde cosinusoïdale de phase $-\pi /4$, et la deuxième est une onde cosinusoïdale de phase $\pi /2$. Échantillonnez le signal à 100 Hz pendant 1 s.

Fs = 100;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*15*t - pi/4) + cos(2*pi*40*t + pi/2);

Calculez la transformée de Fourier du signal. Tracez l’amplitude de la transformée en tant que fonction de la fréquence.

y = fft(x);
z = fftshift(y);

ly = length(y);
f = (-ly/2:ly/2-1)/ly*Fs;
stem(f,abs(z))
title("Double-Sided Amplitude Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("|y|")
grid

Calculez la phase de la transformée en supprimant les valeurs de faible amplitude. Tracez la phase comme une fonction de la fréquence.

tol = 1e-6;
z(abs(z) < tol) = 0;

theta = angle(z);

stem(f,theta/pi)
title("Phase Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("Phase/\pi")
grid

Interpolation de FFT

Interpolez la transformée de Fourier d’un signal en le remplissant avec des zéros.

Spécifiez les paramètres d’un signal avec une fréquence d’échantillonnage de 80 Hz et une durée de signal de 0,8 s.

Fs = 80;
T = 1/Fs;
L = 65;
t = (0:L-1)*T;

Superposez-lui un signal sinusoïdal de 2 Hz et ses harmoniques supérieures. Le signal contient une onde cosinusoïdale de 2 Hz, une onde cosinusoïdale de 4 Hz et une onde sinusoïdale de 6 Hz.

X = 3*cos(2*pi*2*t) + 2*cos(2*pi*4*t) + sin(2*pi*6*t);

Tracez le signal dans le domaine temporel.

plot(t,X)
title("Signal superposition in time domain")
xlabel("t (ms)")
ylabel("X(t)")

Calculez la transformée de Fourier du signal.

Y = fft(X);

Calculez le spectre d’amplitude unilatéral du signal.

f = Fs*(0:(L-1)/2)/L;
P2 = abs(Y/L);
P1 = P2(1:(L+1)/2);
P1(2:end) = 2*P1(2:end);

Dans le domaine fréquentiel, tracez le spectre unilatéral. Comme l’échantillonnage temporel du signal est très court, la résolution en fréquence de la transformée de Fourier n’est pas suffisamment précise pour que la fréquence de crête vers 4 Hz soit visible.

plot(f,P1,"-o") 
title("Single-Sided Spectrum of Original Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Pour mieux évaluer les fréquences de crête, vous pouvez augmenter la longueur de la fenêtre d’analyse en remplissant le signal d’origine avec des zéros. Cette méthode interpole automatiquement la transformée de Fourier du signal avec une résolution en fréquence plus précise.

Identifiez une nouvelle longueur en entrée correspondant à la prochaine puissance de 2 par rapport à la longueur de signal d’origine. Remplissez le signal X avec des zéros à droite pour étendre sa longueur. Calculez la transformée de Fourier du signal rempli avec des zéros.

n = 2^nextpow2(L);
Y = fft(X,n);

Calculez le spectre d’amplitude unilatéral du signal rempli avec des zéros. Comme la longueur de signal n est passée de 65 à 128, la résolution en fréquence passe à Fs/n, c’est-à-dire 0,625 Hz.

f = Fs*(0:(n/2))/n;
P2 = abs(Y/L);
P1 = P2(1:n/2+1);
P1(2:end-1) = 2*P1(2:end-1);

Tracez le spectre unilatéral du signal rempli aves des zéros. Ce nouveau spectre indique les fréquences de crêtes vers 2 Hz, 4 Hz et 6 Hz dans la résolution en fréquence de 0,625 Hz.

plot(f,P1,"-o") 
title("Single-Sided Spectrum of Padded Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

Arguments d'entrée

X — Tableau en entrée

vecteur | matrice | tableau multidimensionnel

Tableau en entrée, spécifié sous forme de vecteur, de matrice ou de tableau multidimensionnel.

Si X est une matrice vide de 0 x 0, fft(X) renvoie une matrice vide de 0 x 0.

Types de données : double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical
Support des nombres complexes : Oui

n — Longueur de la transformée

[] (par défaut) | scalaires de nombres entiers non négatifs

Longueur de la transformée, spécifiée comme [] ou un scalaire entier non négatif. La spécification d’un scalaire entier positif comme longueur de la transformée peut améliorer les performances de fft. La longueur est généralement spécifiée comme une puissance de 2 ou comme une valeur pouvant être factorisée en un produit de petits nombres premiers (avec des facteurs premiers inférieurs ou égaux à 7). Si n est inférieur à la longueur du signal, fft ignore les valeurs de signal restantes, au-delà de la n-ième entrée, et renvoie le résultat tronqué. Si n est égal à 0, fft renvoie une matrice vide.

Exemple : n = 2^nextpow2(size(X,1))

Types de données : double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical

dim — Dimension sur laquelle opérer

scalaire entier positif

Dimension sur laquelle opérer, spécifiée en tant que scalaire entier positif. Si vous ne spécifiez pas la dimension, la valeur par défaut est la première dimension de tableau dont la taille n’est pas égale à 1.

• fft(X,[],1) opère le long des colonnes de X et renvoie la transformée de Fourier de chaque colonne.

  

• fft(X,[],2) opère le long des lignes de X et renvoie la transformée de Fourier de chaque ligne.

  

Si dim est supérieur à ndims(X), fft(X,[],dim) renvoie X. Lorsque n est spécifié, fft(X,n,dim) remplit ou tronque X pour atteindre la longueur n le long de la dimension dim.

Types de données : double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical

Arguments de sortie

Y — Représentation dans le domaine fréquentiel

vecteur | matrice | tableau multidimensionnel

Représentation dans le domaine fréquentiel, renvoyée sous forme de vecteur, de matrice ou de tableau multidimensionnel.

Si X est de type single, fft effectue les calculs de manière native en simple précision et Y est également de type single. Sinon, Y est renvoyé avec le type double.

La taille de Y est la suivante :

• Pour Y = fft(X) ou Y = fft(X,[],dim), la taille de Y est égale à celle de X.

• Pour Y = fft(X,n,dim), la valeur de size(Y,dim) est égale à n tandis que la taille de toutes les autres dimensions restent identiques à celles de X.

Si X est réel, Y est symétrique conjugué et le nombre de points uniques dans Y est ceil((n+1)/2).

Types de données : double | single

En savoir plus

Transformée de Fourier discrète d’un vecteur

Y = fft(X) et X = ifft(Y) implémentent respectivement la transformée de Fourier et la transformée de Fourier inverse. Pour X et Y de longueur n, ces transformées sont définies comme suit :

$$\begin{array}{l}Y(k)={\displaystyle \sum _{j=1}^{n}X}(j)W_n^{(j-1)\text{}(k-1)}\\ X(j)=\frac{1}{n}{\displaystyle \sum _{k=1}^{n}Y}(k)W_n{}^{-(j-1)\text{}(k-1)},\end{array}$$

où

$$W_n=e^{(-2\pi i)/n}$$

est l’une des n racines de l’unité.

Conseils

• La vitesse d’exécution de fft dépend de la longueur de la transformée. L’exécution est beaucoup plus rapide pour les longueurs de transformée ayant seulement de petits facteurs premiers (inférieurs ou égaux à 7) que pour celles qui sont premières ou ont de grands facteurs premiers.

• Pour la plupart des valeurs de n, les DFT à entrées réelles nécessitent environ deux fois moins de temps de calcul que celles à entrées complexes. Toutefois, quand n a de grands facteurs premiers, il n’y a que peu ou pas de différence de vitesse.

• Vous pouvez éventuellement augmenter la vitesse de fft avec la fonction utilitaire fftw. Cette fonction contrôle l’optimisation de l’algorithme utilisé pour calculer une FFT ayant une taille et une dimension particulières.

Algorithmes

Les fonctions FFT (fft, fft2, fftn, ifft, ifft2, ifftn) sont basées sur une bibliothèque nommée FFTW [1] [2].

Références

[1] FFTW (https://www.fftw.org)

[2] Frigo, M., and S. G. Johnson. “FFTW: An Adaptive Software Architecture for the FFT.” Proceedings of the International Conference on Acoustics, Speech, and Signal Processing. Vol. 3, 1998, pp. 1381-1384.

Capacités étendues

Génération de code C/C++

Générez du code C et C++ avec MATLAB® Coder™.

Notes d’usage et limitations :

• Pour l’argument en entrée X :

  • Si vous ne spécifiez pas de dimension, le générateur de code opère le long de la première dimension du tableau en entrée dont la taille est variable ou n’est pas égale à 1. Si cette dimension est de taille variable au moment de la génération de code et égale à 1 au run-time, une erreur run-time peut survenir. Pour éviter cette erreur, spécifiez la dimension.

• Pour une sortie MEX, MATLAB^® Coder™ utilise la bibliothèque dont se sert MATLAB pour les algorithmes FFT. Pour le code C/C++ autonome, le générateur de code produit par défaut du code pour les algorithmes FFT au lieu d’appels de bibliothèque FFT. Pour générer des appels à une bibliothèque FFTW installée spécifique, utilisez une classe de rappel de bibliothèque FFT. Pour plus d’informations sur la classe de rappel de bibliothèque FFT, consultez coder.fftw.StandaloneFFTW3Interface (MATLAB Coder).

• Pour la simulation d’un bloc MATLAB Function, le logiciel de simulation utilise la bibliothèque dont se sert MATLAB pour les algorithmes FFT. Pour la génération de code C/C++, le générateur de code produit par défaut du code pour les algorithmes FFT au lieu d’appels de bibliothèque FFT. Pour générer des appels à une bibliothèque FFTW installée spécifique, utilisez une classe de rappel de bibliothèque FFT. Pour plus d’informations sur la classe de rappel de bibliothèque FFT, consultez coder.fftw.StandaloneFFTW3Interface (MATLAB Coder).

• En utilisant la bibliothèque de remplacement de code (CRL), vous pouvez générer du code optimisé exécutable sur les processeurs ARM^® Cortex^®-A avec l’extension Neon. Pour générer ce code optimisé, vous devez installer Embedded Coder^® Support Package for ARM Cortex-A Processors (Embedded Coder). Le code généré pour ARM Cortex-A utilise la bibliothèque Ne10. Pour plus d’informations, consultez Ne10 Conditions for MATLAB Functions to Support ARM Cortex-A Processors (Embedded Coder) (Conditions Ne10 pour que les fonctions MATLAB supportent les processeurs ARM Cortex-A).

• En utilisant la bibliothèque de remplacement de code (CRL), vous pouvez générer du code optimisé exécutable sur les processeurs ARM Cortex-M. Pour générer ce code optimisé, vous devez installer Embedded Coder Support Package for ARM Cortex-M Processors (Embedded Coder). Le code généré pour ARM Cortex-M utilise la bibliothèque CMSIS. Pour plus d’informations, consultez CMSIS Conditions for MATLAB Functions to Support ARM Cortex-M Processors (Embedded Coder) (Conditions CMSIS pour que les fonctions MATLAB supportent les processeurs ARM Cortex-M).

Génération de code GPU

Générez du code CUDA® pour les GPU NVIDIA® avec GPU Coder™.

Consultez les notes d’usage et limitations à la section « Génération de code C/C++ ». Les mêmes notes d’usage et limitations s’appliquent à la génération de code GPU.

Environnement basé sur les threads

Exécutez du code en arrière-plan avec MATLAB® backgroundPool ou accélérez le code avec Parallel Computing Toolbox™ ThreadPool.

La fonction fft supporte entièrement les environnements basés sur des threads. Pour plus d’informations, consultez Run MATLAB Functions in Thread-Based Environment.

GPU Arrays

Accélérez le code en exécutant les calculs sur une unité de traitement graphique (GPU) avec Parallel Computing Toolbox™.

La fonction fft supporte les GPU arrays avec les notes d’usage et limitations suivantes :

• La sortie Y est toujours complexe, même si toutes les parties imaginaires sont nulles.

Pour plus d’informations, consultez Exécuter les fonctions MATLAB sur un GPU (Parallel Computing Toolbox).

Distributed arrays

Divisez les tableaux volumineux dans la mémoire combinée de votre cluster avec Parallel Computing Toolbox™.

La fonction fft supporte les distributed arrays avec les notes d’usage et limitations suivantes :

• Pour les tableaux distribués, au lieu d’utiliser un algorithme FFT parallèle, fft rassemble les vecteurs sur un seul worker pour exécuter des FFT de longueur première. Pour les FFT de vecteurs ayant une grande longueur première, des erreurs de mémoire insuffisante peuvent survenir.

Pour plus d’informations, consultez Exécuter les fonctions MATLAB avec des tableaux distribués (Parallel Computing Toolbox).

Historique des versions

Introduit avant R2006a

Voir aussi

Fonctions

• fft2 | fftn | fftw | fftshift | ifft | interpft

Rubriques

• Transformées de Fourier

Ressources pédagogiques

• Analyse de Fourier (Ressources pédagogiques MathWorks)
• La convolution en traitement numérique du signal (Ressources pédagogiques MathWorks)