Wykreślanie szybkiej transformaty Fouriera w Pythonie

Mam dostęp do NumPy i SciPy i chcę stworzyć prosty FFT zestawu danych. Mam dwie listy, jedną z ywartościami, a drugą z sygnaturami czasowymi tych ywartości.

Jaki jest najprostszy sposób wprowadzenia tych list do metody SciPy lub NumPy i wykreślenia wynikowego FFT?

Sprawdziłem przykłady, ale wszystkie polegają na utworzeniu zestawu fałszywych danych z pewną liczbą punktów danych, częstotliwością itp. I tak naprawdę nie pokazują, jak to zrobić, mając tylko zestaw danych i odpowiadające im znaczniki czasu .

Wypróbowałem następujący przykład:

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

Ale kiedy zmieniam argument z fftna mój zestaw danych i wykreślam go, otrzymuję bardzo dziwne wyniki i wydaje się, że skalowanie częstotliwości może być wyłączone. Nie jestem pewien.

Oto zestawienie danych, które próbuję wykonać w FFT

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

Kiedy używam fft()całości, ma po prostu ogromny skok na zero i nic więcej.

Oto mój kod:

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

Rozstaw jest po prostu równy xInterp[1]-xInterp[0].

Więc uruchamiam funkcjonalnie równoważną formę twojego kodu w notatniku 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()

Mam to, co uważam za bardzo rozsądną produkcję.

Minęło więcej czasu, niż chciałbym się przyznać, odkąd studiowałem inżynierię, myśląc o przetwarzaniu sygnału, ale skoki na poziomie 50 i 80 są dokładnie tym, czego bym się spodziewał. Więc o co chodzi?

W odpowiedzi na zamieszczone surowe dane i komentarze

Problem polega na tym, że nie masz okresowych danych. Zawsze należy sprawdzać dane, które wprowadzasz do dowolnego algorytmu, aby upewnić się, że są one odpowiednie.

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)

Ważną rzeczą w fft jest to, że można go zastosować tylko do danych, w których znacznik czasu jest jednolity ( tj. Jednolite próbkowanie w czasie, tak jak pokazano powyżej).

W przypadku nierównomiernego pobierania próbek należy użyć funkcji dopasowania danych. Do wyboru jest kilka samouczków i funkcji:

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

Jeśli dopasowanie nie jest opcją, możesz bezpośrednio użyć jakiejś formy interpolacji do interpolacji danych do jednolitego próbkowania:

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

Gdy masz jednolite próbki, będziesz musiał się tylko martwić o deltę czasu ( t[1] - t[0]) swoich próbek. W takim przypadku możesz bezpośrednio skorzystać z funkcji 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()

To powinno rozwiązać twój problem.

Wysoki skok, który masz, wynika z części DC (niezmiennej, tj. Freq = 0) sygnału. To kwestia skali. Jeśli chcesz zobaczyć zawartość częstotliwości innej niż DC, w celu wizualizacji może być konieczne wykreślenie od przesunięcia 1, a nie od przesunięcia 0 FFT sygnału.

Modyfikacja podanego powyżej przykładu przez @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:])

Wykresy wyjściowe:

Innym sposobem jest wizualizacja danych w skali log:

Za pomocą:

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

Pokaże:

W ramach uzupełnienia już udzielonych odpowiedzi, chciałbym zwrócić uwagę, że często ważne jest, aby bawić się rozmiarem pojemników FFT. Warto przetestować kilka wartości i wybrać tę, która ma większy sens w przypadku Twojej aplikacji. Często jest to ta sama wielkość, co liczba próbek. Tak było w przypadku większości udzielonych odpowiedzi i daje świetne i rozsądne wyniki. Na wypadek, gdyby ktoś chciał to zbadać, oto moja wersja kodu:

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

wykresy wyjściowe:

Zbudowałem funkcję, która zajmuje się wykreślaniem FFT rzeczywistych sygnałów. Dodatkową zaletą mojej funkcji w stosunku do poprzednich odpowiedzi jest to, że otrzymujesz rzeczywistą amplitudę sygnału.

Ponadto, ze względu na założenie rzeczywistego sygnału, FFT jest symetryczny, więc możemy wykreślić tylko dodatnią stronę osi 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 preferred 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 frequence 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)

    # Result in frequencies
    fftPlot(sig, dt=dt)
    # Result in samples (if the frequencies axis is unknown)
    fftPlot(sig)

Na tej stronie są już świetne rozwiązania, ale wszyscy założyli, że zbiór danych jest równomiernie / równomiernie próbkowany / dystrybuowany. Postaram się podać bardziej ogólny przykład losowo próbkowanych danych. Jako przykładu użyję również tego samouczka MATLAB :

Dodawanie wymaganych modułów:

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

Generowanie przykładowych danych:

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

Sortowanie zbioru danych:

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

Ponowne próbkowanie:

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)

Wykreślanie danych i ponowne próbkowanie danych:

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

Teraz obliczając 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 W końcu miałem czas na zaimplementowanie bardziej kanonicznego algorytmu, aby uzyskać transformację Fouriera nierównomiernie rozłożonych danych. Możesz zobaczyć kod, opis i przykładowy notatnik Jupyter tutaj .

Piszę tę dodatkową odpowiedź, aby wyjaśnić pochodzenie dyfuzji skoków podczas korzystania z FFT, a zwłaszcza omówić samouczek scipy.fftpack, z którym w pewnym momencie się nie zgadzam.

W tym przykładzie czas nagrywania tmax=N*T=0.75. Sygnał jest sin(50*2*pi*x) + 0.5*sin(80*2*pi*x). Sygnał częstotliwości powinien zawierać dwa skoki o częstotliwościach 50i 80amplitudach 1i 0.5. Jeśli jednak analizowany sygnał nie ma całkowitej liczby okresów, może wystąpić dyfuzja z powodu obcięcia sygnału:

• Pik 1: 50*tmax=37.5=> częstotliwość 50nie jest wielokrotnością 1/tmax=> Obecność dyfuzji z powodu obcięcia sygnału na tej częstotliwości.
• Pik 2: 80*tmax=60=> częstotliwość 80jest wielokrotnością 1/tmax=> Brak dyfuzji z powodu obcięcia sygnału na tej częstotliwości.

Oto kod analizujący ten sam sygnał co w tutorialu ( sin(50*2*pi*x) + 0.5*sin(80*2*pi*x)), ale z niewielkimi różnicami:

1. Oryginalny przykład scipy.fftpack.
2. Oryginalny przykład scipy.fftpack z całkowitą liczbą okresów sygnału ( tmax=1.0zamiast w 0.75celu uniknięcia dyfuzji obcięcia).
3. Oryginalny przykład scipy.fftpack z całkowitą liczbą okresów sygnału oraz datami i częstotliwościami zaczerpniętymi z teorii FFT.

Kod:

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 positioning 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 positioning of dates')
plt.legend()
plt.grid()
plt.show()

Wynik:

Jak to może być tutaj, nawet przy użyciu całkowitej liczby okresów nadal istnieje pewne rozproszenie. To zachowanie jest spowodowane złym ustawieniem dat i częstotliwości w samouczku scipy.fftpack. Stąd w teorii dyskretnych przekształceń Fouriera:

• sygnał należy oceniać w datach, w t=0,T,...,(N-1)*Tktórych T oznacza okres próbkowania, a całkowity czas trwania sygnału wynosi tmax=N*T. Zauważ, że zatrzymujemy się na tmax-T.
• powiązane częstotliwości to f=0,df,...,(N-1)*dfgdzie df=1/tmax=1/(N*T)jest częstotliwość próbkowania. Wszystkie harmoniczne sygnału powinny być wielokrotnością częstotliwości próbkowania, aby uniknąć dyfuzji.

W powyższym przykładzie widać, że użycie arangezamiast linspacepozwala uniknąć dodatkowej dyfuzji w widmie częstotliwości. Co więcej, używanie tej linspacewersji prowadzi również do przesunięcia kolców, które znajdują się na nieco wyższych częstotliwościach niż powinny, jak widać na pierwszym zdjęciu, gdzie kolce znajdują się trochę po prawej stronie częstotliwości 50i 80.

Wnioskuję tylko, że przykład użycia należy zastąpić następującym kodem (który moim zdaniem jest mniej mylący):

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

Wyjście (drugi skok nie jest już rozproszony):

Myślę, że ta odpowiedź nadal przynosi dodatkowe wyjaśnienia dotyczące prawidłowego stosowania dyskretnej transformaty Fouriera. Oczywiście, moja odpowiedź jest zbyt długa i nie zawsze jest więcej rzeczy do powiedzenia ( ewerlopes rozmawiałem krótko o aliasing na przykład i wiele można powiedzieć o okienkowy ), więc będę się zatrzymać.

Myślę, że bardzo ważne jest, aby dogłębnie zrozumieć zasady dyskretnej transformaty Fouriera podczas jej stosowania, ponieważ wszyscy wiemy, że tak wiele osób dodaje czynniki tu i tam, stosując ją, aby uzyskać to, czego chcą.