wake-up-neo.com

Eine schnelle Fourier-Transformation in Python zeichnen

Ich habe Zugriff auf Numpy und Scipy und möchte eine einfache FFT eines Datensatzes erstellen. Ich habe zwei Listen, eine ist Y-Werte und die andere Zeitstempel für diese Y-Werte. 

Was ist der einfachste Weg, um diese Listen in eine scipy- oder numpy-Methode einzugeben und die resultierende FFT zu zeichnen?

Ich habe Beispiele nachgeschlagen, aber alle sind darauf angewiesen, einen Satz gefälschter Daten mit einer bestimmten Anzahl von Datenpunkten, Häufigkeit, usw. zu erstellen, und zeigt nicht wirklich, wie dies mit einem Datensatz und den entsprechenden Zeitstempeln getan wird .

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 fft in meinen Datensatz umwandle und es plottiere, bekomme ich äußerst ungerade Ergebnisse. Es scheint, dass die Skalierung für die Frequenz deaktiviert ist. ich bin mir nicht sicher.

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

http://Pastebin.com/0WhjjMkbhttp://Pastebin.com/ksM4FvZS

Wenn ich die ganze Sache mit einem Fft mache, hat es bei Null einfach einen riesigen Spike 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 nur gleich xInterp[1]-xInterp[0] 

49
user3123955

Ich führe also eine funktional äquivalente Form Ihres Codes in einem IPython-Notebook 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, was ich für sehr vernünftige Ergebnisse halte.

enter image description here

Es ist länger, als ich zugeben möchte, seit ich an 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 Reaktion auf die Veröffentlichung der Rohdaten und Kommentare

Das Problem hier ist, dass Sie keine periodischen Daten haben. Sie sollten immer die Daten überprüfen, die Sie in den any -Algorithmus eingeben, um sicherzustellen, dass sie richtig 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)

enter image description here

62
Paul H

Das Wichtigste an fft ist, dass es nur auf Daten angewendet werden kann, bei denen der Zeitstempel einheitlich ist (d. H. Gleichmäßige zeitliche Abtastung, wie Sie es oben gezeigt haben). 

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

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

Wenn die Anpassung nicht möglich ist, können Sie eine Interpolationsform direkt verwenden, um Daten zu einer einheitlichen Abtastung zu interpolieren: 

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

Wenn Sie einheitliche Proben haben, müssen Sie sich nur noch um das Zeit-Delta (t[1] - t[0]) Ihrer Proben 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. 

16
ssm

Die hohe Spitze, die Sie haben, ist auf den Anteil DC (nicht variierend, d. H. Frequenz = 0) Ihres Signals zurückzuführen. Es ist eine Frage der Größenordnung. Wenn Sie für die Visualisierung einen Nicht-DC-Frequenzinhalt anzeigen möchten, müssen Sie möglicherweise den Offset 1 und nicht den Offset 0 der FFT des Signals darstellen.

Ändern des oben angegebenen Beispiels durch @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 Ausgabeplots: Ploting FFT signal with DC and then when removing it (skipping freq = 0)

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

Mit:

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

Wird zeigen:enter image description here

7
hesham_EE

Ergänzend 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 die Werte auszuwählen, die für Ihre Anwendung sinnvoller sind. Oft liegt es in derselben Größenordnung wie die Anzahl der Abtastungen. Dies wurde von den meisten Antworten angenommen und führt zu guten und vernünftigen Ergebnissen. Für den Fall, dass man das erforschen möchte, hier 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 Ausgabeplots:  enter image description here

4
ewerlopes

Es gibt bereits großartige Lösungen auf dieser Seite, aber alle haben davon ausgegangen, dass der Datensatz einheitlich/gleichmäßig erfasst/verteilt ist. Ich werde versuchen, ein allgemeineres Beispiel für zufällig ausgewählte Daten bereitzustellen. Ich werde auch dieses MATLAB-Tutorial 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

Datensatz sortieren:

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

 enter image description here 

jetzt berechne ich 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)

 enter image description here 

1
Foad

Ich muss eine Funktion erstellen, die sich mit dem Plotten von FFT von echten Signalen befasst. In meiner Funktion gibt es die IST-Amplitude des Signals (wiederum wegen der Annahme eines echten Signals, was Symmetrie bedeutet ...):

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

def fftPlot(sig, dt=None, block=False, 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(block=block)

    return sigFFTPos, freqAxisPos


if __== "__main__":
    dt = 1 / 1000
    f0 = 1 / dt / 4

    t = np.arange(0, 1 + dt, dt)
    sig = np.sin(2 * np.pi * f0 * t)
    fftPlot(sig, dt=dt)
    fftPlot(sig)

    t = np.arange(0, 1 + dt, dt)
    sig = np.sin(2 * np.pi * f0 * t) + 10 * np.sin(2 * np.pi * f0 / 2 * t)
    fftPlot(sig, dt=dt, block=True)
0
asaflotz