Построение преобразования Фурье для функции Гаусса с python, но результат был неверным - PullRequest
1 голос
/ 29 января 2020

Я долго об этом думал, но не понимаю, в чем проблема. Надеюсь, вы можете помочь мне, спасибо.

F(s) Gaussian function
F(s)=1/(√2π s) e^(-(w-μ)^2/(2s^2 ))

Код:

import numpy as np
from matplotlib import pyplot as plt
from math import pi
from scipy.fft import fft

def F_S(w, mu, sig):
    return (np.exp(-np.power(w-mu, 2)/(2 * np.power(sig, 2))))/(np.power(2*pi, 0.5)*sig)

w=np.linspace(-5,5,100)
plt.plot(w, np.real(np.fft.fft(F_S(w, 0, 1))))
plt.show()

Результат:

enter image description here

Ответы [ 4 ]

0 голосов
/ 30 января 2020

Как упоминалось ранее, вы хотите абсолютное значение, а не действительную часть. Минимальный пример, показывающий re / im, abs / phase спектры.

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

n=1001                       # add 1 to keep the interval a round number when using linspace
t = np.linspace(-5, 5, n )   # presumed to be time
dt=t[1]-t[0]                 # time resolution
print(f'sampling every  {dt:.3f} sec , so at {1/dt:.1f} Sa/sec, max. freq will be {1/2/dt:.1f} Hz')


y = np.exp(-(t**2)/0.01)      # signal in time

fr= np.fft.fftshift(np.fft.fftfreq(n, dt))  # shift helps with sorting the frequencies for better plotting
ft=np.fft.fftshift(np.fft.fft(y))           # fftshift only necessary for plotting in sequence

p.figure(figsize=(20,12))
p.subplot(231)
p.plot(t,y,'.-')
p.xlabel('time (secs)')
p.title('signal in time')

p.subplot(232)
p.plot(fr,np.abs(ft), '.-',lw=0.3)      
p.xlabel('freq (Hz)')
p.title('spectrum, abs');

p.subplot(233)
p.plot(fr,np.real(ft), '.-',lw=0.3)     
p.xlabel('freq (Hz)')
p.title('spectrum, real');

p.subplot(235)
p.plot(fr,np.angle(ft), '.-', lw=0.3)    
p.xlabel('freq (Hz)')
p.title('spectrum, phase');

p.subplot(236)
p.plot(fr,np.imag(ft), '.-',lw=0.3)      
p.xlabel('freq (Hz)')
p.title('spectrum, imag');

enter image description here

0 голосов
/ 29 января 2020

Когда вы делаете БПФ, вы получите преобразование simetri c, то есть зеркало кривой от положительного к отрицательному. Как правило, вы только посмотрите на положительную сторону. Кроме того, вам следует позаботиться о частоте дискретизации, так как FFT предназначен для преобразования временной области входных данных в частотную область , время или частота выборки входной информации имеет значение. Так что добавьте временной шаг в np.fft.fftfreq(n, d=timestep) для вашей частоты дискретизации.

Если вы просто хотите сделать частичку нормального сигнала dist, вот еще один вопрос с ним и несколько хороших объяснений того, почему вы получаете такое поведение:

Фурье-преобразование гауссиана не является гауссовским, но это неправильно! - Python

0 голосов
/ 29 января 2020

В вашем коде есть две ошибки:

  • Не принимайте действительного, возьмите абсолютное значение при построении.
  • Из документов:

Если A = fft(a, n), то A[0] содержит член нулевой частоты (среднее значение сигнала), который всегда является чисто действительным для реальных входов. Тогда A[1:n/2] содержит члены с положительной частотой, а A[n/2+1:] содержит члены с отрицательной частотой в порядке убывания отрицательной частоты.

Элементы можно переставить с помощью np.fft.fftshift .

Рабочий код:

import numpy as np
from matplotlib import pyplot as plt
from math import pi
from scipy.fftpack import fft, fftshift

def F_S(w, mu, sig):
    return (np.exp(-np.power(w-mu, 2)/(2 * np.power(sig, 2))))/(np.power(2*pi, 0.5)*sig)

w=np.linspace(-5,5,100)
plt.plot(w, fftshift(np.abs(np.fft.fft(F_S(w, 0, 1)))))
plt.show()

Кроме того, вы также можете рассмотреть возможность масштабирования оси x.

0 голосов
/ 29 января 2020

Вы должны перейти от шкалы времени к шкале частоты

...