Zeichnen einer schnellen Fourier-Transformation in Python

Ich habe Zugriff auf NumPy und SciPy und möchte eine einfache FFT eines Datensatzes erstellen. Ich habe zwei Listen, eine mit yWerten und eine mit Zeitstempeln für diese yWerte.

Was ist der einfachste Weg, diese Listen in eine SciPy- oder NumPy-Methode einzuspeisen und die resultierende FFT zu zeichnen?

Ich habe Beispiele nachgeschlagen, aber alle basieren darauf, einen Satz gefälschter Daten mit einer bestimmten Anzahl von Datenpunkten, einer bestimmten Häufigkeit usw. zu erstellen, und zeigen nicht wirklich, wie dies mit nur einem Satz von Daten und den entsprechenden Zeitstempeln zu tun ist .

Ich habe folgendes Beispiel ausprobiert:

from scipy.fftpack import fft

# Number of samplepoints
N = 600

# Sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
import matplotlib.pyplot as plt
plt.plot(xf, 2.0/N * np.abs(yf[0:N/2]))
plt.grid()
plt.show()

Wenn ich jedoch das Argument von fftin meinen Datensatz ändere und es zeichne, erhalte ich äußerst merkwürdige Ergebnisse, und es scheint, dass die Skalierung für die Frequenz möglicherweise nicht stimmt. Ich bin mir nicht sicher.

Hier ist ein Pastebin der Daten, die ich zu FFT versuche

http://pastebin.com/0WhjjMkb http://pastebin.com/ksM4FvZS

Wenn ich fft()das Ganze benutze, hat es nur eine riesige Spitze bei Null und sonst nichts.

Hier ist mein Code:

## Perform FFT with SciPy
signalFFT = fft(yInterp)

## Get power spectral density
signalPSD = np.abs(signalFFT) ** 2

## Get frequencies corresponding to signal PSD
fftFreq = fftfreq(len(signalPSD), spacing)

## Get positive half of frequencies
i = fftfreq>0

##
plt.figurefigsize = (8, 4));
plt.plot(fftFreq[i], 10*np.log10(signalPSD[i]));
#plt.xlim(0, 100);
plt.xlabel('Frequency [Hz]');
plt.ylabel('PSD [dB]')

Abstand ist gerade gleich xInterp[1]-xInterp[0].

Also führe ich eine funktional äquivalente Form Ihres Codes in einem IPython-Notizbuch aus:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

# Number of samplepoints
N = 600
# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

fig, ax = plt.subplots()
ax.plot(xf, 2.0/N * np.abs(yf[:N//2]))
plt.show()

Ich bekomme das, was ich für eine sehr vernünftige Ausgabe halte.

Es ist länger her, als ich zugeben möchte, seit ich in der Ingenieurschule über Signalverarbeitung nachgedacht habe, aber Spitzen bei 50 und 80 sind genau das, was ich erwarten würde. Also, was ist das Problem?

Als Antwort auf die Rohdaten und Kommentare, die veröffentlicht werden

Das Problem hierbei ist, dass Sie keine periodischen Daten haben. Sie sollten immer die Daten überprüfen, die Sie in einen Algorithmus einspeisen, um sicherzustellen, dass sie angemessen sind.

import pandas
import matplotlib.pyplot as plt
#import seaborn
%matplotlib inline

# the OP's data
x = pandas.read_csv('http://pastebin.com/raw.php?i=ksM4FvZS', skiprows=2, header=None).values
y = pandas.read_csv('http://pastebin.com/raw.php?i=0WhjjMkb', skiprows=2, header=None).values
fig, ax = plt.subplots()
ax.plot(x, y)

Das Wichtige an fft ist, dass es nur auf Daten angewendet werden kann, bei denen der Zeitstempel einheitlich ist ( dh eine einheitliche zeitliche Abtastung, wie oben gezeigt).

Bei ungleichmäßiger Probenahme verwenden Sie bitte eine Funktion zum Anpassen der Daten. Es stehen verschiedene Tutorials und Funktionen zur Auswahl:

https://github.com/tiagopereira/python_tips/wiki/Scipy%3A-curve-fitting http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html

Wenn eine Anpassung nicht möglich ist, können Sie direkt eine Form der Interpolation verwenden, um Daten zu einer einheitlichen Stichprobe zu interpolieren:

https://docs.scipy.org/doc/scipy-0.14.0/reference/tutorial/interpolate.html

Wenn Sie einheitliche Stichproben haben, müssen Sie sich nur um das Zeitdelta ( t[1] - t[0]) Ihrer Stichproben kümmern. In diesem Fall können Sie die fft-Funktionen direkt verwenden

Y    = numpy.fft.fft(y)
freq = numpy.fft.fftfreq(len(y), t[1] - t[0])

pylab.figure()
pylab.plot( freq, numpy.abs(Y) )
pylab.figure()
pylab.plot(freq, numpy.angle(Y) )
pylab.show()

Dies sollte Ihr Problem lösen.

Die hohe Spitze, die Sie haben, ist auf den DC-Anteil (nicht variierend, dh freq = 0) Ihres Signals zurückzuführen. Es ist eine Frage der Größenordnung. Wenn Sie zur Visualisierung Inhalte ohne Gleichstromfrequenz anzeigen möchten, müssen Sie möglicherweise vom Offset 1 und nicht vom Offset 0 der FFT des Signals zeichnen.

Ändern des oben angegebenen Beispiels von @PaulH

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

# Number of samplepoints
N = 600
# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = 10 + np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

plt.subplot(2, 1, 1)
plt.plot(xf, 2.0/N * np.abs(yf[0:N/2]))
plt.subplot(2, 1, 2)
plt.plot(xf[1:], 2.0/N * np.abs(yf[0:N/2])[1:])

Die Ausgabediagramme:

Eine andere Möglichkeit besteht darin, die Daten im Protokollmaßstab zu visualisieren:

Verwenden von:

plt.semilogy(xf, 2.0/N * np.abs(yf[0:N/2]))

Wird zeigen:

Als Ergänzung zu den bereits gegebenen Antworten möchte ich darauf hinweisen, dass es oft wichtig ist, mit der Größe der Behälter für die FFT zu spielen. Es wäre sinnvoll, eine Reihe von Werten zu testen und den Wert auszuwählen, der für Ihre Anwendung sinnvoller ist. Oft ist es in der gleichen Größenordnung wie die Anzahl der Proben. Dies wurde von den meisten Antworten angenommen und führt zu großartigen und vernünftigen Ergebnissen. Für den Fall, dass man das untersuchen möchte, hier ist meine Codeversion:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

fig = plt.figure(figsize=[14,4])
N = 600           # Number of samplepoints
Fs = 800.0
T = 1.0 / Fs      # N_samps*T (#samples x sample period) is the sample spacing.
N_fft = 80        # Number of bins (chooses granularity)
x = np.linspace(0, N*T, N)     # the interval
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)   # the signal

# removing the mean of the signal
mean_removed = np.ones_like(y)*np.mean(y)
y = y - mean_removed

# Compute the fft.
yf = scipy.fftpack.fft(y,n=N_fft)
xf = np.arange(0,Fs,Fs/N_fft)

##### Plot the fft #####
ax = plt.subplot(121)
pt, = ax.plot(xf,np.abs(yf), lw=2.0, c='b')
p = plt.Rectangle((Fs/2, 0), Fs/2, ax.get_ylim()[1], facecolor="grey", fill=True, alpha=0.75, hatch="/", zorder=3)
ax.add_patch(p)
ax.set_xlim((ax.get_xlim()[0],Fs))
ax.set_title('FFT', fontsize= 16, fontweight="bold")
ax.set_ylabel('FFT magnitude (power)')
ax.set_xlabel('Frequency (Hz)')
plt.legend((p,), ('mirrowed',))
ax.grid()

##### Close up on the graph of fft#######
# This is the same histogram above, but truncated at the max frequence + an offset. 
offset = 1    # just to help the visualization. Nothing important.
ax2 = fig.add_subplot(122)
ax2.plot(xf,np.abs(yf), lw=2.0, c='b')
ax2.set_xticks(xf)
ax2.set_xlim(-1,int(Fs/6)+offset)
ax2.set_title('FFT close-up', fontsize= 16, fontweight="bold")
ax2.set_ylabel('FFT magnitude (power) - log')
ax2.set_xlabel('Frequency (Hz)')
ax2.hold(True)
ax2.grid()

plt.yscale('log')

die Ausgabediagramme:

Ich habe eine Funktion erstellt, die sich mit dem Zeichnen der FFT realer Signale befasst. Der zusätzliche Bonus in meiner Funktion in Bezug auf die obigen Nachrichten ist, dass Sie die IST-Amplitude des Signals erhalten. Aufgrund der Annahme eines realen Signals ist die FFT auch symmetrisch, sodass wir nur die positive Seite der x-Achse darstellen können:

import matplotlib.pyplot as plt
import numpy as np
import warnings


def fftPlot(sig, dt=None, plot=True):
    # here it's assumes analytic signal (real signal...)- so only half of the axis is required

    if dt is None:
        dt = 1
        t = np.arange(0, sig.shape[-1])
        xLabel = 'samples'
    else:
        t = np.arange(0, sig.shape[-1]) * dt
        xLabel = 'freq [Hz]'

    if sig.shape[0] % 2 != 0:
        warnings.warn("signal prefered to be even in size, autoFixing it...")
        t = t[0:-1]
        sig = sig[0:-1]

    sigFFT = np.fft.fft(sig) / t.shape[0]  # divided by size t for coherent magnitude

    freq = np.fft.fftfreq(t.shape[0], d=dt)

    # plot analytic signal - right half of freq axis needed only...
    firstNegInd = np.argmax(freq < 0)
    freqAxisPos = freq[0:firstNegInd]
    sigFFTPos = 2 * sigFFT[0:firstNegInd]  # *2 because of magnitude of analytic signal

    if plot:
        plt.figure()
        plt.plot(freqAxisPos, np.abs(sigFFTPos))
        plt.xlabel(xLabel)
        plt.ylabel('mag')
        plt.title('Analytic FFT plot')
        plt.show()

    return sigFFTPos, freqAxisPos


if __name__ == "__main__":
    dt = 1 / 1000

    # build a signal within nyquist - the result will be the positive FFT with actual magnitude
    f0 = 200  # [Hz]
    t = np.arange(0, 1 + dt, dt)
    sig = 1 * np.sin(2 * np.pi * f0 * t) + \
        10 * np.sin(2 * np.pi * f0 / 2 * t) + \
        3 * np.sin(2 * np.pi * f0 / 4 * t) +\
        7.5 * np.sin(2 * np.pi * f0 / 5 * t)

    # res in freqs
    fftPlot(sig, dt=dt)
    # res in samples (if freqs axis is unknown)
    fftPlot(sig)

Es gibt bereits großartige Lösungen auf dieser Seite, aber alle haben angenommen, dass der Datensatz gleichmäßig / gleichmäßig abgetastet / verteilt ist. Ich werde versuchen, ein allgemeineres Beispiel für zufällig ausgewählte Daten bereitzustellen. Ich werde auch dieses MATLAB-Tutorial verwenden als Beispiel verwenden:

Hinzufügen der erforderlichen Module:

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack
import scipy.signal

Beispieldaten generieren:

N = 600 # number of samples
t = np.random.uniform(0.0, 1.0, N) # assuming the time start is 0.0 and time end is 1.0
S = 1.0 * np.sin(50.0 * 2 * np.pi * t) + 0.5 * np.sin(80.0 * 2 * np.pi * t) 
X = S + 0.01 * np.random.randn(N) # adding noise

Sortieren des Datensatzes:

order = np.argsort(t)
ts = np.array(t)[order]
Xs = np.array(X)[order]

Resampling:

T = (t.max() - t.min()) / N # average period 
Fs = 1 / T # average sample rate frequency
f = Fs * np.arange(0, N // 2 + 1) / N; # resampled frequency vector
X_new, t_new = scipy.signal.resample(Xs, N, ts)

Zeichnen der Daten und erneutes Abtasten der Daten:

plt.xlim(0, 0.1)
plt.plot(t_new, X_new, label="resampled")
plt.plot(ts, Xs, label="org")
plt.legend()
plt.ylabel("X")
plt.xlabel("t")

Berechnen Sie jetzt die fft:

Y = scipy.fftpack.fft(X_new)
P2 = np.abs(Y / N)
P1 = P2[0 : N // 2 + 1]
P1[1 : -2] = 2 * P1[1 : -2]

plt.ylabel("Y")
plt.xlabel("f")
plt.plot(f, P1)

PS Ich hatte endlich Zeit, einen kanonischeren Algorithmus zu implementieren, um eine Fourier-Transformation ungleichmäßig verteilter Daten zu erhalten. Sie können den Code, die Beschreibung und das Beispiel eines Jupyter-Notizbuchs hier sehen .

Ich schreibe diese zusätzliche Antwort, um die Ursprünge der Diffusion der Spikes bei der Verwendung von fft zu erklären und insbesondere das scipy.fftpack- Tutorial zu diskutieren, mit dem ich irgendwann nicht einverstanden bin.

In diesem Beispiel die Aufnahmezeit tmax=N*T=0.75. Das Signal ist sin(50*2*pi*x)+0.5*sin(80*2*pi*x). Das Frequenzsignal sollte 2 Spitzen bei Frequenzen 50und 80mit Amplituden 1und enthalten 0.5. Wenn das analysierte Signal jedoch keine ganzzahlige Anzahl von Perioden aufweist, kann aufgrund der Verkürzung des Signals eine Diffusion auftreten:

• Pike 1: 50*tmax=37.5=> Frequenz 50ist kein Vielfaches von 1/tmax=> Vorhandensein von Diffusion aufgrund von Signalabschneidung bei dieser Frequenz.
• Pike 2: 80*tmax=60=> Frequenz 80ist ein Vielfaches von 1/tmax=> Keine Diffusion aufgrund von Signalabschneidung bei dieser Frequenz.

Hier ist ein Code, der das gleiche Signal wie im Tutorial ( sin(50*2*pi*x)+0.5*sin(80*2*pi*x)) analysiert, jedoch mit geringfügigen Unterschieden:

1. Das ursprüngliche Beispiel für scipy.fftpack.
2. Das ursprüngliche Beispiel für scipy.fftpack mit einer ganzzahligen Anzahl von Signalperioden ( tmax=1.0anstatt eine Kürzungsdiffusion 0.75zu vermeiden).
3. Das ursprüngliche Beispiel für scipy.fftpack mit einer ganzzahligen Anzahl von Signalperioden und wobei die Daten und Frequenzen der FFT-Theorie entnommen sind.

Der Code:

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

# 1. Linspace
N = 600
# sample spacing
tmax = 3/4
T = tmax / N # =1.0 / 800.0
x1 = np.linspace(0.0, N*T, N)
y1 = np.sin(50.0 * 2.0*np.pi*x1) + 0.5*np.sin(80.0 * 2.0*np.pi*x1)
yf1 = scipy.fftpack.fft(y1)
xf1 = np.linspace(0.0, 1.0/(2.0*T), N//2)

# 2. Integer number of periods
tmax = 1
T = tmax / N # sample spacing
x2 = np.linspace(0.0, N*T, N)
y2 = np.sin(50.0 * 2.0*np.pi*x2) + 0.5*np.sin(80.0 * 2.0*np.pi*x2)
yf2 = scipy.fftpack.fft(y2)
xf2 = np.linspace(0.0, 1.0/(2.0*T), N//2)

# 3. Correct positionning of dates relatively to FFT theory (arange instead of linspace)
tmax = 1
T = tmax / N # sample spacing
x3 = T * np.arange(N)
y3 = np.sin(50.0 * 2.0*np.pi*x3) + 0.5*np.sin(80.0 * 2.0*np.pi*x3)
yf3 = scipy.fftpack.fft(y3)
xf3 = 1/(N*T) * np.arange(N)[:N//2]

fig, ax = plt.subplots()
# Plotting only the left part of the spectrum to not show aliasing
ax.plot(xf1, 2.0/N * np.abs(yf1[:N//2]), label='fftpack tutorial')
ax.plot(xf2, 2.0/N * np.abs(yf2[:N//2]), label='Integer number of periods')
ax.plot(xf3, 2.0/N * np.abs(yf3[:N//2]), label='Correct positionning of dates')
plt.legend()
plt.grid()
plt.show()

Ausgabe:

Wie es hier sein kann, bleibt auch bei Verwendung einer ganzzahligen Anzahl von Perioden eine gewisse Diffusion bestehen. Dieses Verhalten ist auf eine schlechte Positionierung von Daten und Häufigkeiten im Tutorial scipy.fftpack zurückzuführen. Daher in der Theorie der diskreten Fourier-Transformationen:

• Das Signal sollte an Daten ausgewertet werden, an t=0,T,...,(N-1)*Tdenen T die Abtastperiode und die Gesamtdauer des Signals ist tmax=N*T. Beachten Sie, dass wir bei anhalten tmax-T.
• Die zugehörigen Frequenzen sind dort, f=0,df,...,(N-1)*dfwo df=1/tmax=1/(N*T)die Abtastfrequenz ist. Alle Harmonischen des Signals sollten ein Vielfaches der Abtastfrequenz sein, um eine Diffusion zu vermeiden.

Im obigen Beispiel können Sie sehen, dass die Verwendung von arangeanstelle von linspaceermöglicht, eine zusätzliche Diffusion im Frequenzspektrum zu vermeiden. Darüber hinaus führt die Verwendung der linspaceVersion auch zu einem Versatz der Spitzen, die sich bei etwas höheren Frequenzen befinden als sie sein sollten, wie auf dem ersten Bild zu sehen ist, wo sich die Spitzen etwas rechts von den Frequenzen 50und befinden 80.

Ich komme nur zu dem Schluss, dass das Verwendungsbeispiel durch den folgenden Code ersetzt werden sollte (der meiner Meinung nach weniger irreführend ist):

import numpy as np
from scipy.fftpack import fft
# Number of sample points
N = 600
T = 1.0 / 800.0
x = T*np.arange(N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
yf = fft(y)
xf = 1/(N*T)*np.arange(N//2)
import matplotlib.pyplot as plt
plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
plt.grid()
plt.show()

Ausgabe (die zweite Spitze ist nicht mehr diffus):

Ich denke, diese Antwort bringt noch einige zusätzliche Erklärungen, wie man eine korrekt diskrete Fourier-Transformation anwendet. Offensichtlich ist meine Antwort zu lang und es gibt immer zusätzliche Dinge zu sagen (@ewerlopes hat zum Beispiel kurz über Aliasing gesprochen und es kann viel über Fensterung gesagt werden ), also werde ich aufhören. Ich denke, dass es sehr wichtig ist, die Prinzipien der diskreten Fourier-Transformation bei der Anwendung genau zu verstehen, da wir alle so viele Leute kennen, die hier und da Faktoren hinzufügen, wenn sie angewendet werden, um das zu erhalten, was sie wollen.