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를 시도하는 데이터의 pastebin입니다.
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()
이것은 당신의 문제를 해결할 것입니다.
답변
높은 스파이크는 신호의 DC (비 변함, 즉 주파수 = 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 플로팅을 처리하는 함수를 만들었습니다. 이전 답변과 관련된 내 기능의 추가 보너스 는 신호 의 실제 진폭을 얻는 것 입니다.
또한 실제 신호를 가정하기 때문에 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) +\
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)
답변
이 페이지에는 이미 훌륭한 솔루션이 있지만 모두 데이터 세트가 균일하게 / 균등하게 샘플링 / 분산된다고 가정했습니다. 무작위로 샘플링 된 데이터의보다 일반적인 예를 제공하려고합니다. 또한 이 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)
추신 : 드디어 고르지 않게 분산 된 데이터의 푸리에 변환을 얻기 위해보다 표준적인 알고리즘을 구현할 시간을 얻었습니다. 코드, 설명 및 예제 Jupyter 노트북은 여기에서 볼 수 있습니다 .
답변
나는 FFT를 사용할 때 스파이크 확산의 기원을 설명하고 특히 어느 시점에서 동의하지 않는 scipy.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)) 와 동일한 신호를 분석하는 코드 이지만 약간의 차이가 있습니다.
- 원본 scipy.fftpack 예제.
- 정수 신호주기가있는 원본 scipy.fftpack 예제 ( 잘림 확산 방지
tmax=1.0대신0.75). - 정수 수의 신호주기와 날짜 및 빈도가 FFT 이론에서 가져온 원래 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 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()
산출:
여기에서 볼 수 있듯이 정수 기간을 사용하더라도 일부 확산은 여전히 남아 있습니다. 이 동작은 scipy.fftpack 자습서에서 날짜 및 빈도의 잘못된 위치 지정 때문입니다. 따라서 이산 푸리에 변환 이론에서 :
- 신호는
t=0,T,...,(N-1)*TT가 샘플링 기간이고 신호의 총 지속 시간이 인 날짜에 평가되어야합니다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 간단히 이야기 에 대해 앨리어싱 예를 들어와 많은에 대해 말할 수있다 윈도 내가 중단됩니다, 그래서는).
이산 푸리에 변환을 적용 할 때 이산 푸리에 변환의 원리를 깊이 이해하는 것이 매우 중요하다고 생각합니다. 우리 모두가 원하는 것을 얻기 위해 적용 할 때 여기 저기 요인을 추가하는 많은 사람들을 알고 있기 때문입니다.
