вычисление dFT на частотах БПФ

Я вычисляю dFT функции f(x), выбранной при x_i, i=0,1,...,N (с известным dx) на частотах u_j, j=0,1,...,N, где u_j — частоты, которые генерирует np.fft.fftfreq(N, dx), и сравнивают их с результатом np.fft.fft(f(x)). Я вижу, что эти двое не согласны...

Я что-то упускаю? Разве они по определению не должны быть одинаковыми? (Разница еще хуже, когда я смотрю на части изображения dFT / FFT).

Я прилагаю сценарий, который я использовал, который генерирует этот график, который сравнивает реальную и воображаемую части dFT и FFT.

import numpy as np
import matplotlib.pyplot as plt
from astropy import units

def func_1D(x, sigma_x):
    return np.exp(-(x**2.0 / (2.0 * sigma_x**2)))


n_pixels = int(2**5.0)
pixel_scale = 0.05 # units of arcsec

x_rad = np.linspace(
    -n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) + pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    +n_pixels * pixel_scale / 2.0 * (units.arcsec).to(units.rad) - pixel_scale / 2.0 * (units.arcsec).to(units.rad),
    n_pixels)


sigma_x = 0.5 # in units of arcsec
image = func_1D(
    x=x_rad,
    sigma_x=sigma_x * units.arcsec.to(units.rad),
)
image_FFT = np.fft.fftshift(np.fft.fft(np.fft.fftshift(image)))
u_grid = np.fft.fftshift(np.fft.fftfreq(n_pixels, d=pixel_scale * units.arcsec.to(units.rad)))

image_dFT = np.zeros(shape=n_pixels, dtype="complex")
for i in range(u_grid.shape[0]):
    for j in range(n_pixels):
        image_dFT[i] += image[j] * np.exp(
            -2.0
            * np.pi
            * 1j
            * (u_grid[i] * x_rad[j])
        )

value = 0.23

figure, axes = plt.subplots(nrows=1,ncols=3,figsize=(14,6))
axes[0].plot(x_rad * 10**6.0, image, marker="o")
for x_i in x_rad:
    axes[0].axvline(x_i * 10**6.0, linestyle="--", color="black")
axes[0].set_xlabel(r"x ($\times 10^6$; rad)")
axes[0].set_title("x-plane")

for u_grid_i in u_grid:
    axes[1].axvline(u_grid_i / 10**6.0, linestyle="--", color="black")
axes[1].plot(u_grid / 10**6.0, image_FFT.real, color="b")
axes[1].plot(u_grid / 10**6.0, image_dFT.real, color="r", linestyle="None", marker="o")
axes[1].set_title("u-plane (real)")
axes[1].set_xlabel(r"u ($\times 10^{-6}$; rad$^{-1}$)")
axes[1].plot(u_grid / 10**6.0, image_FFT.real - image_dFT.real, color="black", label="difference")

axes[2].plot(u_grid / 10**6.0, image_FFT.imag, color="b")
axes[2].plot(u_grid / 10**6.0, image_dFT.imag, color="r", linestyle="None", marker="o")
axes[2].set_title("u-plane (imag)")
axes[2].set_xlabel(r"u ($\times 10^{-6}$; rad$^{-1}$)")
#axes[2].plot(u_grid / 10**6.0, image_FFT.imag - image_dFT.imag, color="black", label="difference")
axes[1].legend()
plt.show()

Sketos 06.02.2020 источник

comment

Если ваш ввод реален (я бы позволил гауссиане перейти к низким значениям), вам следует искать абсолютное значение в спектре. Пожалуйста, прочитайте минимальный пример. Я бы пропустил единицы, некоторые из этих манипуляций не так просто выполнить. - roadrunner66 08.02.2020

Ответы (2)

arrow_upward
1
arrow_downward

Я сделал минимальный пример (надеюсь). Я получаю практически те же числа для БПФ и наивного интеграла Фурье (вычисленного для тех же значений частоты).

import numpy as np
import matplotlib.pyplot as p 
%matplotlib inline

def signal(x, sigma_x):
    return np.exp(-(x**2.0 / (2.0 * sigma_x**2)))

t=np.linspace(-10,10,1000)
sigma=.3
sig=np.exp(-(t**2.0 / (2.0 * sigma **2)))

p.subplot(311)
p.plot(t,sig);

ft=np.fft.fftshift(np.fft.fft(sig))
freq=np.fft.fftshift(np.fft.fftfreq(1000,0.02))
p.subplot(312)
p.plot(freq,np.abs(ft))
print(np.abs(ft)[500:505])
# naive fourier integral 
fi=[]
for f in freq: 
  i=np.sum( sig* np.exp(- 1j* 2 *np.pi*f*t ))
  fi.append(np.abs(i))

p.subplot(313)
p.plot(freq,fi)

print(np.abs(fi)[500:505])

roadrunner66 08.02.2020

comment

Вам нужно построить как действительные, так и мнимые компоненты. Получение величины скрывает проблемы, вызванные смещением сигнала (например, если использование fftshift неверно). - Cris Luengo; 08.02.2020

comment

@roadrunner66 спасибо за ваш пример, это было очень полезно. Я обновил его ниже, чтобы вместо этого показать реальную и мнимую части БПФ. - Sketos; 09.02.2020

arrow_upward
0
arrow_downward

Я обновил пример @roadrunner66, чтобы вместо этого показать реальную и мнимую части FT, а не величину, поскольку приложение, для которого я хочу его использовать, включает в себя работу с реальными и мнимыми частями FT (обычно называемые видимостью в интерферометрия).

Ниже приведен немного обновленный пример.

import numpy as np
import matplotlib.pyplot as plt

t=np.linspace(-10,10,1000)
sigma=.3
sig=np.exp(-(t**2.0 / (2.0 * sigma **2)))

ft=np.fft.fftshift(np.fft.fft(sig))
freq=np.fft.fftshift(np.fft.fftfreq(len(t),abs(t[0] - t[1])))

# naive fourier integral
fi_real=[]
fi_imag=[]
for f in freq:
  i=np.sum( sig* np.exp(- 1j* 2 *np.pi*f*t ))
  fi_real.append(i.real)
  fi_imag.append(i.imag)

figure, axes = plt.subplots(nrows=1,ncols=2)
axes[0].plot(freq,ft.real, color="b", label="np.fft.fft")
axes[0].plot(freq,fi_real, color="r", label="exact")
axes[0].set_xlim(-5.0, 5.0)
axes[0].set_title("real")
axes[0].legend()
axes[1].plot(freq,ft.imag, color="b", label="np.fft.fft")
axes[1].plot(freq,fi_imag, color="r", label="exact")
axes[1].set_xlim(-5.0, 5.0)
axes[1].set_title("imag")
axes[1].legend()
plt.show()

Глядя на выходной рисунок, я думаю, становится ясно, что np.fft не подходит, когда вы хотите работать с реальной и мнимой частями БПФ.

Sketos 09.02.2020

comment

Это интересно. Не могли бы вы дать ссылку? Есть, конечно, случаи, когда интересными физическими величинами являются только действительная или мнимая часть (зависит только от определения), но я не знал об этом для реальной или мнимой части электрического поля, где приложения Im осведомлены представляют собой смесь действительных и мнимых частей после любого БПФ (поле, мощность, FT (линейная автокорреляция во времени) = спектр интенсивности, нелинейная автокорреляция во времени, видимость пространственной полосы ‹-› длина когерентности и многие другие). - roadrunner66; 09.02.2020

comment

Единственное, что вам нужно изменить в моем примере, это построить real и imag из abs. Нет необходимости генерировать чисто реальные или чисто мнимые векторы в numpy, пока вам это не понадобится. Сама идея сложного FT состоит в том, чтобы не писать все дважды. - roadrunner66; 09.02.2020

вычисление dFT на частотах БПФ

Ответы (2)

Похожие вопросы