Python에서 빠른 푸리에 변환 플롯팅

Python에서 빠른 푸리에 변환 플롯팅

NumPy 및 SciPy에 액세스할 수 있으며 데이터 세트의 간단한 FFT를 만들고 싶습니다.y입니다에 입니다.y가치.

이러한 목록을 SciPy 또는 NumPy 방법으로 공급하고 결과 FFT를 표시하는 가장 간단한 방법은 무엇입니까?

예를 찾아봤지만, 이들은 모두 일정한 수의 데이터 포인트와 빈도 등으로 가짜 데이터 집합을 만드는 데 의존하고 있으며, 데이터 집합과 해당 타임스탬프만으로 수행하는 방법을 보여주지는 않습니다.

저는 다음 예시를 시도해 보았습니다.

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제 데이터 세트와 플롯에 따르면, 저는 매우 이상한 결과를 얻었으며, 주파수에 대한 스케일링이 꺼져 있을 수 있습니다.확신이 서지 않습니다.

FFT를 시도하는 데이터의 페이스트빈입니다.

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

가 사용할 때.fft()전체적으로 볼 때 0에 큰 스파이크가 있을 뿐 다른 것은 없습니다.

여기 내 코드가 있습니다.

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

그러면 문제가 해결될 겁니다.

저는 실제 신호의 FFT 플롯을 다루는 함수를 만들었습니다.앞의 답변에 비해 제 기능에서 추가적인 이점은 신호의 실제 진폭을 얻을 수 있다는 것입니다.

또한 실제 신호를 가정하기 때문에 FFT는 대칭이므로 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)
        + 10 * np.sin(2 * np.pi * (f0 * 2 + 0.5) * t)  # <--- not sampled on grid so the peak will not be actual height
    )
    # Result in frequencies
    fftPlot(sig, dt=dt)
    # Result in samples (if the frequencies axis is unknown)
    fftPlot(sig)

높은 스파이크는 신호의 DC(비 varying, 즉 주파수 = 0) 부분 때문입니다.규모의 문제입니다.비DC 주파수 컨텐츠를 보려면 시각화를 위해 신호의 FFT의 오프셋 0이 아닌 오프셋 1부터 플롯해야 할 수도 있습니다.

@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]))

:

이미 제시된 답변을 보완하는 것처럼, FFT를 위해 빈의 크기를 가지고 플레이하는 것이 중요하다는 점을 지적하고자 합니다.여러 값을 테스트하고 응용 프로그램에 더 적합한 값을 선택하는 것이 타당할 것입니다.표본 수와 동일한 크기인 경우가 많습니다.이는 제시된 대부분의 답변에서 가정한 것으로, 훌륭하고 합리적인 결과를 도출합니다.만약 누군가 그것을 탐색하고 싶다면, 여기 제 코드 버전이 있습니다.

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

출력 그림:

FFT를 사용할 때 스파이크가 확산되는 원인을 설명하고 특히 스키피에 대해 논의하기 위해 이 추가 답변을 씁니다.어느 시점에서 내가 동의하지 않는 fftpack 튜토리얼.

tmax=N*T=0.75는. 는.sin(50*2*pi*x) + 0.5*sin(80*2*pi*x)에서 두 . .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. 원래 스키피.fft pack 예제.
2. 원래 스키피.신호 주기가 정수인 fft pack 예제(tmax=1.0대신에0.75절단 확산을 방지합니다.
3. 원래 스키피.신호 주기의 정수를 포함하고 FFT 이론에서 날짜와 주파수를 가져온 FFT 팩 예제.

코드:

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

출력:

여기에 있을 수 있듯이 정수 개의 주기를 사용하더라도 일부 확산은 여전히 남아 있습니다.이 동작은 스키피에서 날짜와 빈도의 잘못된 위치 때문입니다.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()

출력(두 번째 스파이크는 더 이상 확산되지 않음):

이 답은 여전히 이산 푸리에 변환을 정확하게 적용하는 방법에 대한 추가적인 설명을 가져온다고 생각합니다.분명히 제 대답은 너무 길고 항상 추가적인 할 말이 있습니다(예를 들어 에일리어싱에 대해 짧은 이야기를 했고 윈도잉에 대해 많은 이야기를 할 수 있음). 그래서 그만 하겠습니다.

이산 푸리에 변환을 적용할 때 이산 푸리에 변환의 원리를 깊이 이해하는 것은 매우 중요하다고 생각합니다. 왜냐하면 우리는 사람들이 원하는 것을 얻기 위해 여기저기 인자를 추가하는 것을 매우 많이 알고 있기 때문입니다.

이 페이지에는 이미 훌륭한 솔루션이 있지만 모두 데이터셋이 균일하게/균일하게 샘플링/분산되어 있다고 가정했습니다.무작위로 추출한 데이터의 보다 일반적인 예를 들어보겠습니다.또한 이 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)

추신. 저는 마침내 불균등하게 분포된 데이터의 푸리에 변환을 얻기 위해 보다 표준적인 알고리즘을 구현할 시간을 가졌습니다.당신은 여기서 코드와 설명 그리고 쥬피터 노트의 예시를 볼 수 있습니다.

언급URL : https://stackoverflow.com/questions/25735153/plotting-a-fast-fourier-transform-in-python