Построение быстрого преобразования Фурье в Python

У меня есть доступ к NumPy и SciPy, и я хочу создать простой БПФ набора данных. У меня есть два списка, один - это yзначения, а другой - временные метки для этих yзначений.

Каков самый простой способ передать эти списки в метод SciPy или NumPy и построить результат БПФ?

Я просмотрел примеры, но все они полагаются на создание набора поддельных данных с некоторым определенным количеством точек данных, частотой и т. Д. И на самом деле не показывают, как это сделать с помощью только набора данных и соответствующих временных меток. .

Я пробовал следующий пример:

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

Но когда я меняю аргумент fftсвоего набора данных и строю его график, я получаю чрезвычайно странные результаты, и, похоже, масштабирование частоты может быть отключено. Я не уверен.

Вот вставка данных, которые я пытаюсь выполнить БПФ

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

Когда я использую fft()все это, у него просто огромный всплеск на нуле и ничего больше.

Вот мой код:

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

Интервал просто равен xInterp[1]-xInterp[0].

Поэтому я запускаю функционально эквивалентную форму вашего кода в записной книжке 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()

Я получаю то, что считаю очень разумным.

С тех пор, как я учился в инженерной школе и думал об обработке сигналов, прошло больше времени, чем я хотел бы признать, но всплески на уровне 50 и 80 - это именно то, что я ожидал. Так в чем проблема?

В ответ на публикацию необработанных данных и комментариев

Проблема здесь в том, что у вас нет периодических данных. Вы всегда должны проверять данные, которые вы вводите в любой алгоритм, чтобы убедиться, что они подходят.

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)

Важным моментом в fft является то, что он может применяться только к данным, в которых временная метка единообразна ( т.е. единообразная выборка по времени, как вы показали выше).

В случае неоднородной выборки используйте функцию подбора данных. Есть несколько обучающих программ и функций на выбор:

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

Если подгонка не подходит, вы можете напрямую использовать некоторую форму интерполяции для интерполяции данных до однородной выборки:

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

Когда у вас есть однородные образцы, вам нужно будет беспокоиться только о временной дельте ( t[1] - t[0]) ваших образцов. В этом случае вы можете напрямую использовать функции 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()

Это должно решить вашу проблему.

Высокий пик, который у вас есть, связан с постоянным током (неизменным, т.е. freq = 0) вашего сигнала. Это вопрос масштаба. Если вы хотите видеть частотное содержимое, отличное от постоянного тока, для визуализации вам может потребоваться построить график от смещения 1, а не от смещения 0 БПФ сигнала.

Изменение приведенного выше примера с помощью @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:])

Выходные графики:

Другой способ - визуализировать данные в масштабе журнала:

С помощью:

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

Покажет:

В качестве дополнения к уже приведенным ответам я хотел бы отметить, что часто важно поиграть с размером ячеек для БПФ. Было бы разумно протестировать набор значений и выбрать то, которое больше подходит для вашего приложения. Часто оно равно количеству отсчетов. Это предполагалось в большинстве представленных ответов и дает отличные и разумные результаты. Если кто-то хочет это изучить, вот моя версия кода:

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

графики вывода:

Я построил функцию, которая занимается построением БПФ реальных сигналов. Дополнительным преимуществом моей функции по сравнению с сообщениями выше является то, что вы получаете ФАКТИЧЕСКУЮ амплитуду сигнала. Кроме того, из-за предположения о реальном сигнале БПФ является симметричным, поэтому мы можем построить только положительную сторону оси 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)

На этой странице уже есть отличные решения, но все предполагают, что набор данных равномерно / равномерно отобран / распределен. Я постараюсь привести более общий пример данных, выбранных случайным образом. Я также буду использовать это руководство по MATLAB в качестве примера:

Добавление необходимых модулей:

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

Создание выборки данных:

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

Сортировка набора данных:

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

Передискретизация:

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)

построение графика данных и данных передискретизации:

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

теперь вычисляем 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 Наконец-то у меня появилось время реализовать более канонический алгоритм для получения преобразования Фурье неравномерно распределенных данных. Вы можете увидеть код, описание и пример Jupyter ноутбук здесь .

Я пишу этот дополнительный ответ, чтобы объяснить происхождение распространения шипов при использовании fft, и особенно обсуждаю руководство scipy.fftpack, с которым я в какой-то момент не согласен.

В этом примере время записи tmax=N*T=0.75. Сигнал есть sin(50*2*pi*x)+0.5*sin(80*2*pi*x). Частотный сигнал должен содержать 2 пика на частотах 50и 80с амплитудами 1и 0.5. Однако, если анализируемый сигнал не имеет целого числа периодов, может появиться диффузия из-за усечения сигнала:

• Пайк 1: 50*tmax=37.5=> частота 50не кратна 1/tmax=> Наличие диффузии из-за усечения сигнала на этой частоте.
• Пайк 2: 80*tmax=60=> частота 80кратна 1/tmax=> Нет диффузии из-за усечения сигнала на этой частоте.

Вот код, который анализирует тот же сигнал, что и в учебнике ( sin(50*2*pi*x)+0.5*sin(80*2*pi*x)), но с небольшими отличиями:

1. Исходный пример scipy.fftpack.
2. Исходный пример scipy.fftpack с целым числом периодов сигнала ( tmax=1.0вместо того, 0.75чтобы избежать распространения усечения).
3. Исходный пример scipy.fftpack с целым числом периодов сигнала, где даты и частоты взяты из теории БПФ.

Код:

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

Выход:

Как это может быть здесь, даже при использовании целого числа периодов некоторая диффузия все же остается. Такое поведение связано с неправильным позиционированием дат и частот в учебнике scipy.fftpack. Следовательно, в теории дискретных преобразований Фурье:

• сигнал должен оцениваться в те даты, t=0,T,...,(N-1)*Tгде T - период выборки, а общая длительность сигнала равна tmax=N*T. Обратите внимание, что мы останавливаемся на tmax-T.
• соответствующие частоты , f=0,df,...,(N-1)*dfгде df=1/tmax=1/(N*T)частота дискретизации. Все гармоники сигнала должны быть кратны частоте дискретизации, чтобы избежать диффузии.

В приведенном выше примере вы можете видеть, что использование arangeвместо linspaceпозволяет избежать дополнительной диффузии в частотном спектре. Более того, использование linspaceверсии также приводит к смещению выбросов, которые расположены на несколько более высоких частотах, чем они должны быть, как это видно на первом рисунке, где выбросы находятся немного справа от частот 50и 80.

Я просто заключу, что пример использования следует заменить следующим кодом (который, на мой взгляд, менее вводит в заблуждение):

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

Вывод (второй шип уже не разлетается):

Я думаю, что этот ответ по-прежнему дает некоторые дополнительные объяснения о том, как правильно применять дискретное преобразование Фурье. Очевидно, мой ответ слишком длинный, и всегда есть что сказать (@ewerlopes, например, кратко рассказал о псевдониме, и многое можно сказать об оконном режиме ), поэтому я остановлюсь. Я думаю, что очень важно глубоко понимать принципы дискретного преобразования Фурье при его применении, потому что мы все знаем, как много людей добавляют факторы здесь и там, применяя его, чтобы получить то, что они хотят.