Traçando uma transformação rápida de Fourier em Python

Tenho acesso a NumPy e SciPy e desejo criar um FFT simples de um conjunto de dados. Tenho duas listas, uma de yvalores e a outra de timestamps para esses yvalores.

Qual é a maneira mais simples de alimentar essas listas em um método SciPy ou NumPy e plotar o FFT resultante?

Eu pesquisei exemplos, mas todos eles dependem da criação de um conjunto de dados falsos com um certo número de pontos de dados e frequência, etc. e não mostram realmente como fazer isso com apenas um conjunto de dados e os carimbos de data / hora correspondentes .

Tentei o seguinte exemplo:

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

Mas quando eu mudo o argumento de fftpara meu conjunto de dados e o ploto, obtenho resultados extremamente estranhos e parece que a escala da frequência pode estar errada. Não tenho certeza.

Aqui está um pastebin dos dados que estou tentando FFT

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

Quando eu uso fft()tudo isso, ele tem um grande pico de zero e nada mais.

Aqui está o meu código:

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

O espaçamento é igual a xInterp[1]-xInterp[0].

Então, executo uma forma funcionalmente equivalente do seu código em um 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()

Recebo o que acredito ser uma saída bastante razoável.

Faz mais tempo do que gostaria de admitir desde que estava na faculdade de engenharia pensando em processamento de sinal, mas picos de 50 e 80 são exatamente o que eu esperava. Qual é o problema?

Em resposta aos dados brutos e comentários postados

O problema aqui é que você não tem dados periódicos. Você sempre deve inspecionar os dados que alimenta em qualquer algoritmo para certificar-se de que sejam apropriados.

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)

O importante sobre fft é que ele só pode ser aplicado a dados nos quais o carimbo de data / hora é uniforme ( isto é, amostragem uniforme no tempo, como o que você mostrou acima).

No caso de amostragem não uniforme, use uma função para ajustar os dados. Existem vários tutoriais e funções para escolher:

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

Se o ajuste não for uma opção, você pode usar diretamente alguma forma de interpolação para interpolar dados para uma amostragem uniforme:

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

Quando você tem amostras uniformes, só precisa se preocupar com o delta de tempo ( t[1] - t[0]) de suas amostras. Neste caso, você pode usar diretamente as funções 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()

Isso deve resolver seu problema.

O pico alto que você tem é devido à porção DC (não variável, isto é, freq = 0) do seu sinal. É uma questão de escala. Se você quiser ver o conteúdo de frequência não DC, para visualização, pode ser necessário traçar a partir do deslocamento 1 e não do deslocamento 0 do FFT do sinal.

Modificando o exemplo dado acima por @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:])

A saída plota:

Outra forma é visualizar os dados em escala logarítmica:

Usando:

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

Mostrará:

Apenas como um complemento às respostas já dadas, gostaria de salientar que muitas vezes é importante brincar com o tamanho das caixas para o FFT. Faria sentido testar vários valores e escolher aquele que faz mais sentido para sua aplicação. Freqüentemente, é da mesma magnitude do número de amostras. Isso foi presumido pela maioria das respostas dadas e produz resultados excelentes e razoáveis. Caso alguém queira explorar isso, aqui está minha versão de código:

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

os gráficos de saída:

Eu construí uma função que trata da plotagem de FFT de sinais reais. O bônus extra em minha função em relação às respostas anteriores é que você obtém a amplitude real do sinal.

Além disso, devido ao pressuposto de um sinal real, o FFT é simétrico, portanto, podemos representar graficamente apenas o lado positivo do eixo 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)

Já existem ótimas soluções nesta página, mas todos assumiram que o conjunto de dados é uniformemente / uniformemente amostrado / distribuído. Tentarei fornecer um exemplo mais geral de dados de amostra aleatória. Também usarei este tutorial do MATLAB como exemplo:

Adicionando os módulos necessários:

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

Gerando dados de amostra:

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

Classificando o conjunto de dados:

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

Reamostragem:

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)

Traçando os dados e dados reamostrados:

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

Agora calculando o 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 Finalmente consegui tempo para implementar um algoritmo mais canônico para obter uma transformação de Fourier de dados distribuídos de maneira desigual. Você pode ver o código, a descrição e o bloco de notas Jupyter de exemplo aqui .

Escrevo esta resposta adicional para explicar as origens da difusão dos picos ao usar FFT e, especialmente, discuto o tutorial scipy.fftpack do qual discordo em algum ponto.

Neste exemplo, o tempo de gravação tmax=N*T=0.75. O sinal é sin(50*2*pi*x) + 0.5*sin(80*2*pi*x). O sinal de frequência deve conter dois picos em frequências 50e 80com amplitudes 1e 0.5. No entanto, se o sinal analisado não tiver um número inteiro de períodos, a difusão pode aparecer devido ao truncamento do sinal:

• Pike 1: 50*tmax=37.5=> frequência 50não é um múltiplo de 1/tmax=> Presença de difusão devido ao truncamento do sinal nesta frequência.
• Pike 2: 80*tmax=60=> frequência 80é um múltiplo de 1/tmax=> Sem difusão devido ao truncamento do sinal nesta frequência.

Aqui está um código que analisa o mesmo sinal do tutorial ( sin(50*2*pi*x) + 0.5*sin(80*2*pi*x)), mas com pequenas diferenças:

1. O exemplo original do scipy.fftpack.
2. O exemplo scipy.fftpack original com um número inteiro de períodos de sinal (em tmax=1.0vez de 0.75para evitar a difusão do truncamento).
3. O exemplo scipy.fftpack original com um número inteiro de períodos de sinal e onde as datas e frequências são tiradas da teoria FFT.

O código:

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

Resultado:

Como pode ser aqui, mesmo usando um número inteiro de períodos, alguma difusão ainda permanece. Esse comportamento é devido a um posicionamento incorreto de datas e frequências no tutorial scipy.fftpack. Portanto, na teoria das transformadas discretas de Fourier:

• o sinal deve ser avaliado em datas em t=0,T,...,(N-1)*Tque T é o período de amostragem e a duração total do sinal é tmax=N*T. Observe que paramos em tmax-T.
• as frequências associadas são f=0,df,...,(N-1)*dfonde df=1/tmax=1/(N*T)está a frequência de amostragem. Todos os harmônicos do sinal devem ser múltiplos da frequência de amostragem para evitar a difusão.

No exemplo acima, você pode ver que o uso de em arangevez de linspacepermite evitar difusão adicional no espectro de frequência. Além disso, o uso da linspaceversão também leva a um deslocamento dos picos que estão localizados em frequências ligeiramente mais altas do que deveriam ser, como pode ser visto na primeira foto, onde os picos estão um pouco à direita das frequências 50e 80.

Vou apenas concluir que o exemplo de uso deve ser substituído pelo seguinte código (que é menos enganoso na minha opinião):

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

Resultado (o segundo pico não é mais difundido):

Acho que essa resposta ainda traz algumas explicações adicionais sobre como aplicar corretamente a transformada discreta de Fourier. Obviamente, minha resposta é muito longa e sempre há coisas adicionais a dizer ( ewerlopes falou brevemente sobre aliasing, por exemplo, e muito pode ser dito sobre janelamento ), então vou parar.

Eu acho que é muito importante entender profundamente os princípios da transformada discreta de Fourier ao aplicá-la, porque todos nós sabemos que muitas pessoas adicionam fatores aqui e ali ao aplicá-la para obter o que desejam.