как следует пересчитывать точки fft, чтобы получить те же результаты, что и аналитическое решение?

Question

Я хочу работать с быстрым преобразованием Фурье, используя пакет nfty fft, и затем я пытаюсь сравнить результаты между аналитическим решением и быстрым преобразованием Фурье, и хотя я могу видеть на графиках, которые я сделал, кривыепохожи, довольно очевидно, что масштабы разные.

Я испробовал несколько разных версий частоты (угловая частота, частота и волновое число), но все мои попытки не сработали, и в пустой документации неясно, насколько точно быстрое преобразование Фурьеопределены.Например, я хочу работать с преобразованием Фурье экспоненциальной функции во времени в угловую частотную область, f (t) = Exp (-a | t |), F (w) = a / pi * (a² + w²) (существует несколько версий этого аналитического решения в зависимости от того, какое частотное пространство мы рассматриваем)

def e(t):
    return np.exp(-0.5*abs(t))
def F(w):
    return 0.5/(np.pi)*(1/(((0.5)**2)+((w)**2)))

t=np.linspace(0,100,1000)

w=np.fft.fftfreq(len(t))
plt.plot(w,F(w),'o',label='F(w)')
plt.legend()
plt.show()

fourier=np.fft.fft(e(t))
plt.plot(w,fourier,'o')
plt.show()

Я пробовал несколько различных вариантов вышеуказанного кода специально для частоты, но я все еще не дошел доточка, где БПФ и аналитическое решение похожи.Может ли кто-нибудь помочь мне?

Answer 1

Преобразование Фурье может применяться к интегрируемым функциям, таким как np.exp(-0.5*abs(t)). Но Дискретное преобразование Фурье вычисляет преобразование Фурье периодических сигналов. См. https://dsp.stackexchange.com/questions/26884/about-fourier-transform-of-periodic-signal и Что на самом деле вычисляет FFTW .

Следовательно, ДПФ кадра длины T соответствует преобразованию Фурье периодизированного кадра. Поскольку кадр начинается с 0, вычисляется преобразование Фурье периодического экспоненциального затухания справа: Как видите, половина функции np.exp(-0.5*abs(t)) не отображается. Давайте исправим ее и добавим периодическую возрастающую часть двустороннего экспоненциального затухания. Я использовал частоту в качестве параметра:

import matplotlib.pyplot as plt
import numpy as np

def e(t):
    return np.exp(-0.5*abs(t))
def F(w):
    return 0.5/(np.pi)*(1/(((0.5)**2)+((w)**2)))

def Fc(xi):
    #ok , that's sourced from https://en.wikipedia.org/wiki/Fourier_transform ... Square-integrable functions, one-dimensional, line 207
    return 2*0.5/(((0.5)**2)+(4*(np.pi*xi)**2))


framelength=100.
nbsample=1000
def ep(t):
    #the periodized negative part is added at the end of the frame.
    return np.maximum(np.exp(-0.5*abs(t)),np.exp(-0.5*abs(t-framelength)))

t=np.linspace(0,framelength,nbsample, endpoint=False)

#plotting the periodized signal, to see what's happening
ein=ep(t)
tp=np.linspace(0,10*framelength,10*nbsample, endpoint=False)
periodized=np.zeros(10*nbsample)
for i in range(10):
    for j in range(nbsample):
       periodized[i*nbsample+j]=ein[j]

plt.plot(tp,periodized,'k-',label='periodized frame')
plt.legend()
plt.show()

fourier=np.fft.fft(ep(t))/np.size(ep(t))*framelength

#comparing the mean is useful to check the overall scaling
print np.mean(ep(t))*framelength
print fourier[0]
print Fc(0)

#the frenquencies of the DFT of a frame of length T are 1/T, 2/T ... and negative for the second part of the array.
xi=np.fft.fftfreq(len(t), framelength/len(t))

# comparison between analytical Fourier transform and dft.
plt.plot(xi,Fc(xi),'o',label='F(xi)')
plt.plot(xi,np.real(fourier),'k-', lw=3, color='red', label='DTF')
plt.legend()
plt.show()

А вот и результат:

Для экспериментального непериодического сигнала появляется искусственный разрыв при периодизации кадра. Он вызывает спектральную утечку и окна применяются для ослабления разрыва и его эффектов. Одно из потенциальных окон, называемое окном Пуассона, является двусторонним экспоненциальным затуханием!