[python] Python에서 고속 푸리에 변환 플로팅

나는 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)) 와 동일한 신호를 분석하는 코드 이지만 약간의 차이가 있습니다.

1. 원본 scipy.fftpack 예제.
2. 정수 신호주기가있는 원본 scipy.fftpack 예제 ( 잘림 확산 방지 tmax=1.0대신 0.75).
3. 정수 수의 신호주기와 날짜 및 빈도가 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 간단히 이야기 에 대해 앨리어싱 예를 들어와 많은에 대해 말할 수있다 윈도 내가 중단됩니다, 그래서는).

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