Merencanakan transformasi Fourier cepat dengan Python

Saya memiliki akses ke NumPy dan SciPy dan ingin membuat FFT sederhana dari kumpulan data. Saya punya dua daftar, satu adalah ynilai dan yang lainnya adalah cap waktu untuk ynilai - nilai itu.

Apa cara termudah untuk memasukkan daftar ini ke dalam metode SciPy atau NumPy dan memplot FFT yang dihasilkan?

Saya telah mencari contoh, tetapi semuanya bergantung pada pembuatan satu set data palsu dengan sejumlah titik data, frekuensi, dll. Dan tidak benar-benar menunjukkan cara melakukannya hanya dengan satu set data dan stempel waktu yang sesuai .

Saya telah mencoba contoh berikut:

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

Tetapi ketika saya mengubah argumen fftke kumpulan data saya dan memplotnya, saya mendapatkan hasil yang sangat aneh, dan tampaknya penskalaan untuk frekuensi mungkin tidak aktif. Saya tidak yakin

Ini adalah pastebin dari data yang saya coba FFT

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

Ketika saya menggunakan fft()semuanya, itu hanya memiliki lonjakan besar di nol dan tidak ada yang lain.

Ini kode saya:

## 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]')

Spasi sama dengan xInterp[1]-xInterp[0].

Jadi saya menjalankan bentuk kode Anda yang secara fungsional setara di notebook IPython:

%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()

Saya mendapatkan apa yang saya yakini sebagai hasil yang sangat wajar.

Sudah lebih lama dari yang saya akui sejak saya masih di sekolah teknik memikirkan tentang pemrosesan sinyal, tetapi lonjakan pada 50 dan 80 persis seperti yang saya harapkan. Jadi apa masalahnya?

Menanggapi data mentah dan komentar yang diposting

Masalahnya di sini adalah Anda tidak memiliki data berkala. Anda harus selalu memeriksa data yang Anda masukkan ke dalam algoritme apa pun untuk memastikannya sesuai.

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)

Hal penting tentang fft adalah bahwa ini hanya dapat diterapkan pada data yang stempel waktunya seragam ( yaitu pengambilan sampel seragam dalam waktu, seperti yang Anda tunjukkan di atas).

Dalam kasus pengambilan sampel yang tidak seragam, gunakan fungsi untuk menyesuaikan data. Ada beberapa tutorial dan fungsi yang bisa dipilih:

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

Jika pemasangan bukan suatu pilihan, Anda dapat langsung menggunakan beberapa bentuk interpolasi untuk menginterpolasi data ke pengambilan sampel yang seragam:

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

Jika Anda memiliki sampel yang seragam, Anda hanya perlu mengkhawatirkan delta waktu (t[1] - t[0] ) sampel Anda. Dalam hal ini, Anda dapat langsung menggunakan fungsi fft

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

Ini seharusnya menyelesaikan masalah Anda.

Lonjakan tinggi yang Anda miliki adalah karena bagian DC (tidak bervariasi, yaitu freq = 0) dari sinyal Anda. Ini masalah skala. Jika Anda ingin melihat konten frekuensi non-DC, untuk visualisasi, Anda mungkin perlu membuat plot dari offset 1 bukan dari offset 0 FFT sinyal.

Mengubah contoh yang diberikan di atas oleh @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:])

Plot keluaran:

Cara lain, adalah dengan memvisualisasikan data dalam skala log:

Menggunakan:

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

Akan menunjukkan:

Hanya sebagai pelengkap dari jawaban yang telah diberikan, saya ingin menunjukkan bahwa sering kali penting untuk bermain dengan ukuran nampan untuk FFT. Masuk akal untuk menguji sekumpulan nilai dan memilih yang lebih masuk akal untuk aplikasi Anda. Seringkali, besarnya sama dengan jumlah sampel. Ini seperti yang diasumsikan oleh sebagian besar jawaban yang diberikan, dan menghasilkan hasil yang bagus dan masuk akal. Jika seseorang ingin menjelajahinya, berikut adalah versi kode saya:

%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')

plot keluaran:

Saya telah membangun sebuah fungsi yang berhubungan dengan memplot FFT dari sinyal nyata. Bonus tambahan dalam fungsi saya relatif terhadap pesan di atas adalah Anda mendapatkan amplitudo sinyal yang SEBENARNYA. Selain itu, karena asumsi sinyal nyata, FFT simetris sehingga kita hanya dapat memplot sisi positif sumbu x:

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)

Sudah ada solusi hebat di halaman ini, tetapi semua berasumsi bahwa kumpulan data diambil sampelnya secara seragam / merata / didistribusikan. Saya akan mencoba memberikan contoh yang lebih umum dari data sampel secara acak. Saya juga akan menggunakan tutorial MATLAB ini sebagai contoh:

Menambahkan modul yang diperlukan:

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

Menghasilkan data sampel:

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

Menyortir kumpulan data:

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

Pengambilan sampel ulang:

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)

merencanakan data dan data sampel ulang:

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

sekarang menghitung 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 Saya akhirnya punya waktu untuk menerapkan algoritma yang lebih kanonik untuk mendapatkan transformasi Fourier dari data yang tidak terdistribusi secara merata. Anda mungkin melihat kode, deskripsi, dan contoh notebook Jupyter di sini .

Saya menulis jawaban tambahan ini untuk menjelaskan asal mula difusi paku saat menggunakan fft dan terutama membahas scipy.fftpack tutorial yang saya tidak setuju di beberapa titik.

Dalam contoh ini, waktu perekaman tmax=N*T=0.75. Sinyalnya adalah sin(50*2*pi*x)+0.5*sin(80*2*pi*x). Sinyal frekuensi harus berisi 2 lonjakan pada frekuensi 50dan 80dengan amplitudo 1dan 0.5. Namun, jika sinyal yang dianalisis tidak memiliki bilangan bulat periode, difusi dapat muncul karena pemotongan sinyal:

• Pike 1: 50*tmax=37.5=> frekuensi 50bukan kelipatan 1/tmax=> Adanya difusi akibat pemotongan sinyal pada frekuensi ini.
• Pike 2: 80*tmax=60=> frekuensi 80adalah kelipatan 1/tmax=> Tidak ada difusi akibat pemotongan sinyal pada frekuensi ini.

Berikut adalah kode yang menganalisis sinyal yang sama seperti di tutorial ( sin(50*2*pi*x)+0.5*sin(80*2*pi*x)) tetapi dengan sedikit perbedaan:

1. Contoh asli scipy.fftpack.
2. Contoh scipy.fftpack asli dengan nomor integer periode sinyal ( tmax=1.0bukan0.75 untuk menghindari difusi pemotongan).
3. Contoh asli scipy.fftpack dengan bilangan bulat periode sinyal dan di mana tanggal dan frekuensi diambil dari teori FFT.

Kode:

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

Keluaran:

Seperti yang bisa terjadi di sini, bahkan dengan menggunakan periode bilangan bulat beberapa difusi masih tetap ada. Perilaku ini disebabkan oleh posisi tanggal dan frekuensi yang buruk dalam tutorial scipy.fftpack. Karenanya, dalam teori transformasi Fourier diskrit:

• sinyal harus dievaluasi pada tanggal di t=0,T,...,(N-1)*Tmana T adalah periode pengambilan sampel dan total durasi sinyal tmax=N*T. Perhatikan bahwa kami berhenti di tmax-T.
• frekuensi terkait adalah di f=0,df,...,(N-1)*dfmana df=1/tmax=1/(N*T)frekuensi sampling. Semua harmonisa sinyal harus multipel dari frekuensi sampling untuk menghindari difusi.

Dalam contoh di atas, Anda dapat melihat bahwa penggunaan arangealih-alih linspacememungkinkan untuk menghindari difusi tambahan dalam spektrum frekuensi. Selain itu, menggunakan linspaceversi juga mengarah pada offset paku yang terletak pada frekuensi yang sedikit lebih tinggi dari yang seharusnya seperti yang dapat dilihat pada gambar pertama di mana paku sedikit di sebelah kanan frekuensi 50dan 80.

Saya hanya akan menyimpulkan bahwa contoh penggunaan harus diganti dengan kode berikut (yang menurut saya kurang menyesatkan):

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

Output (lonjakan kedua tidak tersebar lagi):

Saya rasa jawaban ini masih memberikan beberapa penjelasan tambahan tentang bagaimana menerapkan transformasi Fourier diskrit dengan benar. Jelas sekali, jawaban saya terlalu panjang dan selalu ada hal tambahan untuk dikatakan (@ewerlopes berbicara singkat tentang aliasing misalnya dan banyak yang bisa dikatakan tentang windowing ) jadi saya akan berhenti. Menurut saya, sangat penting untuk memahami secara mendalam prinsip transformasi Fourier diskrit saat menerapkannya karena kita semua tahu banyak orang menambahkan faktor di sana-sini saat menerapkannya untuk mendapatkan apa yang mereka inginkan.