고스 함수는 정적 분포의 문법으로 평균 $x0달러 및 분산 $\sigma 달러 사용

f(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp \Bigl({-\frac{(x-x_0)^2}{2 \sigma^2}} \Bigr)

이번에는 고스 함수의 부립엽 변환에 대해 이론 공식의 계산 결과를 FFT(빠른 부립엽 변환)의 결과와 비교한다.

카탈로그

・이론 공식에 기초한 부립엽 변환
・FFT 기반 부립엽 변환
• 양자의 비교
• FFT 실행 코드

이론 공식에 기초한 부립엽 변환

계산 과정이 생략되었고 고스 함수의 부립엽 변환 $F가 고스 함수로 바뀌었습니다.

\begin{eqnarray}
F(k)
&=& \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(x)e^{-ikx}dx \\
&=& \sigma e^{-\frac{\sigma^2k^2}{2}}e^{-ikx_0} \\ \\
\end{eqnarray}

진폭 스펙트럼 $|F|=\sigmae^{-\rac{\sigma^2k^2} 및,
$k\geq0달러로 위상(편각의 주값) ${rm Arg} (F) = {rmArg} (e^{ikx 0}) $를 그립니다.

위쪽은 진폭 스펙트럼이고, 아래쪽은 위상이다.0달러라고 쓰여 있지만 이번에는 $[-\pi\pi]의 주치를 채택하여 가장자리를 넘어 뛰어내렸다. 또한 진폭 스펙트럼이 떨어진 후의 위상 정보는 의미가 없기 때문에 무시했다.
(단, 이 도표의 가로축은 파수로 직접 각주파수로 바꿀 수 있으므로 진동주파수로 읽으려면 $2\pi달러를 나누어야 한다.)

FFT 기반 부립엽 변환

FFT의 집행 방법에 대한 설명은 별도로 설명할 것이다. 여기서 이론식과 FFT의 결과가 일치하는지 검증할 것이다. 나는 파이톤의 NumPy 모듈로 집행했다(마지막으로 참고 코드를 기재했다).
방금 FFT 결과에서 얻은 진폭 스펙트럼과 위상도를 표시합니다.

스펙트럼의 최대치는 아까와 달리 대체적으로 모양이 비슷하다. 참고로 데이터를 많이 찾지 않으면 좋은 결과가 나오지 않는다.

양자의 비교

두 결과가 일치하는지 검증하다.
진폭 스펙트럼 추출비, 위상 추출차 비교.
이론값과 FFT에서 얻은 값은 각각 아래 표 1, 2로 표시한다.

$k\leq4에서 비례가 고정되어 있고 차이는 0입니다. FFT에서 얻은 진폭 스펙트럼의 최대치는 다르지만 전체적으로 몇 배만 정해져 있기 때문에 정성적으로는 문제가 없다고 할 수 있습니다.
이에 따라 고스 함수 부립엽 변환의 이론식과 FFT를 비교한 결과 일치가 확인됐다.

FFT 실행을 위한 코드(Python)

참고로 이번에 FFT를 실행할 때 사용한 코드를 게재했습니다.

import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl

sigma = 1.0
dx = 0.01
x0 = 5
x = np.arange(0, 1000, dx)
N = len(x)

def y(x):
    return 1 / (np.sqrt(2 * np.pi) * sigma) \
           * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))

fourier = np.fft.fft(y(x))
fourier_abs = np.abs(fourier)
fourier_phase = np.arctan2(np.real(fourier), np.imag(fourier))
omega = np.linspace(0, 1/dx, N) * 2 * np.pi

fig = plt.figure()

ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)

ax1.plot(k, np.abs(fourier), linewidth=1, color='black')
ax1.set_ylabel('Amplitude')

ax2.plot(k, np.arctan2(np.real(fourier), np.imag(fourier)) / np.pi, linewidth=1, color='black')
ax2.set_xlabel('k')
ax2.set_ylabel('Phase[rad]')

ax1.set_xlim([0.0, 6.0])
ax2.set_xlim([0.0, 6.0])

mpl.rcParams['axes.xmargin'] = 0
mpl.rcParams['axes.ymargin'] = 0

plt.show()