Hauptinhalt

fft

Schnelle Fourier-Transformation

Beschreibung

Y = fft(X) berechnet die diskrete Fourier-Transformation (DFT) von X mithilfe eines schnellen Fourier-Transformations-Algorithmus (FFT). Y hat dieselbe Größe wie X.

  • Wenn X ein Vektor ist, gibt fft(X) die Fourier-Transformation des Vektors zurück.

  • Wenn X eine Matrix ist, behandelt fft(X) die Spalten von X als Vektoren und gibt die Fourier-Transformation jeder Spalte zurück.

  • Wenn X ein mehrdimensionales Array ist, behandelt fft(X) die Werte entlang der ersten Array-Dimension mit Größe ungleich 1 als Vektoren und gibt die Fourier-Transformation jedes Vektors zurück.

Beispiel

Y = fft(X,n) gibt die n-Punkt-DFT zurück.

  • Wenn X ein Vektor ist und die Länge von X weniger als n ist, wird X mit nachstehenden Nullstellen auf die Länge n aufgefüllt.

  • Wenn X ein Vektor ist und die Länge von X mehr als n ist, wird X auf die Länge von n gekürzt.

  • Wenn X eine Matrix ist, wird jede Spalte wie im Vektorfall behandelt.

  • Wenn X ein mehrdimensionales Array ist, wird die erste Array-Dimension, deren Größe ungleich 1 ist, wie im Vektorfall behandelt.

Beispiel

Y = fft(X,n,dim) gibt die Fourier-Transformation entlang der Dimension dim zurück. Wenn X beispielsweise eine Matrix ist, gibt fft(X,n,2) die n-Punkt-Fourier-Transformation jeder Zeile zurück.

Beispiel

Beispiele

alle reduzieren

Finden Sie die Frequenzkomponenten eines verrauschten Signals und finden Sie die Amplituden der Spitzenfrequenzen mithilfe der Fourier-Transformation.

Geben Sie die Parameter eines Signals mit einer Abtastfrequenz von 1 kHz und einer Signaldauer von 1,5 Sekunden an.

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

Bilden Sie ein Signal mit einem DC-Offset mit Amplitude 0,8, einer 50-Hz-Sinuskurve mit Amplitude 0,7 und einer 120-Hz-Sinuskurve mit Amplitude 1.

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

Verrauschen Sie das Signal mit einem mittelwertfreien Zufallsrauschen mit Varianz 4.

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

Plotten Sie das verrauschte Signal in der Zeitdomäne. Die Frequenzkomponenten sind im Diagramm nicht sichtbar.

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

Figure contains an axes object. The axes object with title Signal Corrupted with Zero-Mean Random Noise, xlabel t (milliseconds), ylabel X(t) contains an object of type line.

Berechnen Sie die Fourier-Transformation des Signals.

Y = fft(X);

Da eine Fourier-Transformation mit komplexen Zahlen vorgenommen wird, plotten Sie den komplexen Betrag des fft-Spektrums.

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

Figure contains an axes object. The axes object with title Complex Magnitude of fft Spectrum, xlabel f (Hz), ylabel |fft(X)| contains an object of type line.

Das Diagramm zeigt fünf Frequenzspitzen einschließlich der Spitze des DC-Offsets bei 0 Hz. In diesem Beispiel werden für das Signal drei Frequenzspitzen bei 0 Hz, 50 Hz und 120 Hz erwartet. Hier ist die zweite Hälfte des Diagramms das Spiegelbild der ersten Hälfte ohne die Spitze bei 0 Hz. Der Grund hierfür ist, dass die diskrete Fourier-Transformation eines Signals im Zeitbereich eine periodische Natur aufweist, bei der die erste Hälfte des Spektrums in den positiven Frequenzen und die zweite Hälfte in den negativen Frequenzen liegt, wobei das erste Element für die Nullfrequenz reserviert ist.

Bei reellen Signalen ist das fft-Spektrum ein zweiseitiges Spektrum, wobei das Spektrum in den positiven Frequenzen das komplexe Konjugat des Spektrums in den negativen Frequenzen ist. Um das fft-Spektrum in den positiven und negativen Frequenzen darzustellen, können Sie fftshift verwenden. Für eine gerade Länge von L beginnt die Frequenzdomäne mit der negativen Nyquist-Frequenz -Fs/2 und reicht bis Fs/2-Fs/L, mit einem Abstand (Frequenzauflösung) von 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)|")

Figure contains an axes object. The axes object with title fft Spectrum in the Positive and Negative Frequencies, xlabel f (Hz), ylabel |fft(X)| contains an object of type line.

Um die Amplituden der drei Frequenzspitzen zu finden, konvertieren Sie das fft-Spektrum in Y in das einseitige Amplitudenspektrum. Da die Funktion fft einen Skalierungsfaktor L zwischen dem ursprünglichen und dem transformierten Signal enthält, skalieren Sie Y durch Teilen durch L neu. Nehmen Sie den komplexen Betrag des fft-Spektrums. Das zweiseitige Amplitudenspektrum P2, bei dem das Spektrum in den positiven Frequenzen das komplexe Konjugat des Spektrums in den negativen Frequenzen ist, weist die Hälfte der Spitzenamplituden des Zeitdomänen-Signals auf. Nehmen Sie zum Konvertieren in das einseitige Spektrum die erste Hälfte des zweiseitigen Spektrums P2. Multiplizieren Sie das Spektrum in den positiven Frequenzen mit 2. Sie müssen P1(1) und P1(end) nicht mit 2 multiplizieren, da diese Amplituden den Null- bzw. Nyquist-Frequenzen entsprechen und keine komplexen Konjugatpaare in den negativen Frequenzen aufweisen.

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

Definieren Sie die Frequenzdomäne f für das einseitige Spektrum. Plotten Sie das einseitige Amplitudenspektrum P1. Wie erwartet, liegen die Amplituden nah an 0,8, 0,7 und 1, sind aber aufgrund des zusätzlichen Rauschens nicht genau. In den meisten Fällen erzeugen längere Signale bessere Frequenzannäherungen.

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)|")

Figure contains an axes object. The axes object with title Single-Sided Amplitude Spectrum of X(t), xlabel f (Hz), ylabel |P1(f)| contains an object of type line.

Nehmen Sie nun die Fourier-Transformation des ursprünglichen, nicht verrauschten Signals und ermitteln Sie die präzisen Amplituden bei 0,8, 0,7 und 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)|")

Figure contains an axes object. The axes object with title Single-Sided Amplitude Spectrum of S(t), xlabel f (Hz), ylabel |P1(f)| contains an object of type line.

Konvertieren Sie einen Gaußimpuls von der Zeitdomäne in die Frequenzdomäne.

Geben Sie die Parameter eines Signals mit einer Abtastfrequenz von 44,1 kHz und einer Signaldauer von 1 ms an. Erstellen Sie einen Gaußimpuls mit einer Standardabweichung von 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)));

Plotten Sie den Impuls in der Zeitdomäne.

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

Figure contains an axes object. The axes object with title Gaussian Pulse in Time Domain, xlabel Time (t), ylabel X(t) contains an object of type line.

Die Ausführungszeit von fft hängt von der Länge der Transformation ab. Transformationslängen mit nur kleinen Primfaktoren resultieren in einer erheblich schnelleren Ausführungszeit als Transformationslängen mit großen Primfaktoren.

In diesem Beispiel beträgt die Signallänge L 44.101, eine sehr große Primzahl. Um die Leistung von fft zu verbessern, ermitteln Sie eine Eingabelänge, die der nächsten Zweierpotenz ausgehend von der ursprünglichen Signallänge entspricht. Wenn Sie fft mit dieser Eingabelänge aufrufen, wird der Impuls X mit nachstehenden Nullstellen auf die angegebene Transformationslänge aufgefüllt.

n = 2^nextpow2(L);

Konvertieren Sie den Gaußimpuls in die Frequenzdomäne.

Y = fft(X,n);

Definieren Sie die Frequenzdomäne und plotten Sie die eindeutigen Frequenzen.

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)|")

Figure contains an axes object. The axes object with title Gaussian Pulse in Frequency Domain, xlabel f (Hz), ylabel |P(f)| contains an object of type line.

Vergleichen Sie Kosinuswellen in der Zeitdomäne und der Frequenzdomäne.

Geben Sie die Parameter eines Signals mit einer Abtastfrequenz von 1 kHz und einer Signaldauer von 1 Sekunde an.

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

Erstellen Sie eine Matrix, in der jede Zeile einer Kosinuswelle mit skalierter Frequenz entspricht. Das Ergebnis X ist eine 3-mal-1000-Matrix. Die erste Zeile weist eine Wellenfrequenz von 50 auf, die zweite Zeile eine Wellenfrequenz von 150 und die dritte Zeile eine Wellenfrequenz von 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];

Plotten Sie die ersten 100 Einträge jeder Zeile von X in Reihenfolge auf einer Abbildung und vergleichen Sie die Frequenzen.

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

Figure contains 3 axes objects. Axes object 1 with title Row 1 in the Time Domain contains an object of type line. Axes object 2 with title Row 2 in the Time Domain contains an object of type line. Axes object 3 with title Row 3 in the Time Domain contains an object of type line.

Geben Sie das Argument dim an, um fft entlang den Zeilen von X (also für jedes Signal) zu verwenden.

dim = 2;

Berechnen Sie die Fourier-Transformation der Signale.

Y = fft(X,L,dim);

Berechnen Sie das zweiseitige Spektrum und das einseitige Spektrum jedes Signals.

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

Plotten Sie in der Frequenzdomäne das einseitige Amplitudenspektrum für jede Zeile auf einer Abbildung.

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

Figure contains 3 axes objects. Axes object 1 with title Row 1 in the Frequency Domain contains an object of type line. Axes object 2 with title Row 2 in the Frequency Domain contains an object of type line. Axes object 3 with title Row 3 in the Frequency Domain contains an object of type line.

Erstellen Sie ein Signal, das aus zwei Sinuskurven der Frequenzen 15 Hz und 40 Hz besteht. Die erste Sinuskurve ist eine Kosinuswelle mit der Phase -π/4 und die zweite ist eine Kosinuswelle mit der Phase π/2. Tasten Sie das Signal 1 Sekunde lang bei 100 Hz ab.

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

Berechnen Sie die Fourier-Transformation des Signals. Plotten Sie den Betrag der Transformation als Funktion der Frequenz.

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

Figure contains an axes object. The axes object with title Double-Sided Amplitude Spectrum of x(t), xlabel Frequency (Hz), ylabel |y| contains an object of type stem.

Berechnen Sie die Phase der Transformation und entfernen Sie dabei Transformationswerte mit kleinem Betrag. Plotten Sie die Phase als Funktion der Frequenz.

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

Figure contains an axes object. The axes object with title Phase Spectrum of x(t), xlabel Frequency (Hz), ylabel Phase/ pi contains an object of type stem.

Interpolieren Sie die Fourier-Transformation eines Signals durch Auffüllen mit Nullstellen.

Geben Sie die Parameter eines Signals mit einer Abtastfrequenz von 80 kHz und einer Signaldauer von 0,8 ms an.

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

Erstellen Sie eine Überlagerung eines 2-Hz-Sinuswellensignals und dessen höheren Oberwellen. Das Signal umfasst eine 2-Hz-Kosinuswelle, eine 4-Hz-Kosinuswelle und eine 6-Hz-Sinuswelle.

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

Plotten Sie das Signal in der Zeitdomäne.

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

Figure contains an axes object. The axes object with title Signal superposition in time domain, xlabel t (ms), ylabel X(t) contains an object of type line.

Berechnen Sie die Fourier-Transformation des Signals.

Y = fft(X);

Berechnen Sie das einseitige Amplitudenspektrum des Signals.

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

Plotten Sie in der Frequenzdomäne das einseitige Spektrum. Da die Zeitabtastung des Signals recht kurz ist, ist die Frequenzauflösung der Fourier-Transformation nicht präzise genug, um die Spitzenfrequenz in der Nähe von 4 Hz darzustellen.

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

Figure contains an axes object. The axes object with title Single-Sided Spectrum of Original Signal, xlabel f (Hz), ylabel |P1(f)| contains an object of type line.

Um die Spitzenfrequenzen besser zu beurteilen, können Sie die Länge des Analysefensters steigern, indem Sie das Originalsignal mit Nullstellen auffüllen. Diese Methode interpoliert automatisch die Fourier-Transformation des Signals mit einer präziseren Frequenzauflösung.

Ermitteln Sie eine neue Eingabelänge, die der nächsten Zweierpotenz ausgehend von der ursprünglichen Signallänge entspricht. Füllen Sie das Signal X mit nachstehenden Nullstellen auf, um dessen Länge zu steigern. Berechnen Sie die Fourier-Transformation des mit Nullstellen aufgefüllten Signals.

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

Berechnen Sie das einseitige Amplitudenspektrum des aufgefüllten Signals. Da die Signallänge n von 65 auf 128 gestiegen ist, wird die Frequenzauflösung Fs/n; dies entspricht 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);

Plotten Sie das einseitige Spektrum des aufgefüllten Signals. Dieses neue Spektrum zeigt die Spitzenfrequenzen in Nähe von 2 Hz, 4 Hz und 6 Hz innerhalb der Frequenzauflösung von 0,625 Hz.

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

Figure contains an axes object. The axes object with title Single-Sided Spectrum of Padded Signal, xlabel f (Hz), ylabel |P1(f)| contains an object of type line.

Eingabeargumente

alle reduzieren

Eingabearray, angegeben als Vektor, Matrix oder mehrdimensionales Array.

Wenn X eine leere 0-mal-0-Matrix ist, gibt fft(X) eine leere 0-mal-0-Matrix zurück.

Datentypen: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical
Unterstützung komplexer Zahlen: Ja

Transformationslänge, angegeben als [] oder als nichtnegativer ganzzahliger Skalar. Einen positiven ganzzahligen Skalar für die Transformationslänge anzugeben, kann die Leistung von fft verbessern. Die Länge wird üblicherweise als Potenz von 2 oder Wert angegeben, der in ein Produkt kleiner Primzahlen (mit Primfaktoren kleiner gleich 7) faktoriert werden kann. Wenn n kleiner als die Länge des Signals ist, ignoriert fft die verbleibenden Signalwerte nach dem n. Eintrag und gibt das gekürzte Ergebnis zurück. Wenn n gleich 0 ist, gibt fft eine leere Matrix zurück.

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

Datentypen: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical

Dimension, entlang der die Operation erfolgt, angegeben als positiver ganzzahliger Skalar. Wenn Sie die Dimension nicht angeben, wird standardmäßig die erste Array-Dimension verwendet, deren Größe nicht 1 ist.

  • fft(X,[],1) operiert entlang den Spalten von X und gibt die Fourier-Transformation jeder Spalte zurück.

    fft(X,[],1) column-wise operation

  • fft(X,[],2) operiert entlang den Spalten von X und gibt die Fourier-Transformation jeder Zeile zurück.

    fft(X,[],2) row-wise operation

Wenn dim größer als ndims(X) ist, gibt fft(X,[],dim) X zurück. Ist n angegeben, gleicht fft(X,n,dim) X auf die Länge n entlang Dimension dim an.

Datentypen: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical

Ausgabeargumente

alle reduzieren

Frequenzdomänen-Darstellung, zurückgegeben als Vektor, Matrix oder mehrdimensionales Array.

Wenn X den Typ single aufweist, berechnet fft nativ mit einfacher Genauigkeit, und Y weist ebenfalls den Typ single auf. Andernfalls wird Y als Typ double zurückgegeben.

Die Größe von Y wird folgendermaßen festgelegt:

  • Für Y = fft(X) oder Y = fft(X,[],dim) ist die Größe von Y identisch zur Größe von X.

  • Für Y = fft(X,n,dim) ist der Wert von size(Y,dim) identisch zu n; die Größe der anderen Dimensionen verbleibt wie in X.

Wenn X reell ist, ist Y konjugatsymmetrisch und die Anzahl eindeutiger Punkte in Y ist ceil((n+1)/2).

Datentypen: double | single

Mehr über

alle reduzieren

Tipps

  • Die Ausführungszeit von fft hängt von der Länge der Transformation ab. Transformationslängen mit nur kleinen Primfaktoren (nicht größer als 7) resultieren in einer erheblich schnelleren Ausführungszeit als Transformationslängen, die selbst Primzahlen sind oder große Primfaktoren aufweisen.

  • Für die meisten Werte von n benötigt eine DFT an einer reellen Eingabe etwa halb so viel Rechenzeit wie eine DFT an einer komplexen Eingabe. Wenn n jedoch große Primfaktoren aufweist, zeigt sich kein oder nur ein geringer Geschwindigkeitsunterschied.

  • Unter Umständen können die Geschwindigkeit von fft mithilfe der Hilfsfunktion fftw steigern. Diese Funktion steuert die Optimierung des Algorithmus, der zur Berechnung der FFT bei einer bestimmten Größe und Dimension verwendet wird.

Algorithmen

Die FFT-Funktionen (fft, fft2, fftn, ifft, ifft2, ifftn) basieren auf einer Bibliothek namens FFTW [1] [2].

Referenzen

[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.

Erweiterte Fähigkeiten

alle erweitern

GPU-Codegenerierung
Generieren von CUDA® Code für NVIDIA® Grafikprozessoren mit dem GPU Coder™.

Versionsverlauf

Eingeführt vor R2006a