Pregunta La transformada de Fourier de un gaussiano no es gaussiana, ¡pero eso está mal! - Python


Estoy tratando de utilizar la función fft de Numpy, sin embargo, cuando le doy a la función una función gausiana simple, la función gausiana no es gausiana, está cerca pero está dividida a la mitad, de modo que cada mitad está en cada extremo del eje x.

La función gaussiana que estoy calculando es y = exp (-x ^ 2)

Aquí está mi código:

from cmath import *
from numpy import multiply
from numpy.fft import fft
from pylab import plot, show

""" Basically the standard range() function but with float support """
def frange (min_value, max_value, step):
    value = float(min_value)
    array = []
    while value < float(max_value):
        array.append(value)
        value += float(step)
    return array


N = 256.0 # number of steps
y = []
x = frange(-5, 5, 10/N)

# fill array y with values of the Gaussian function   
cache = -multiply(x, x)
for i in cache: y.append(exp(i))

Y = fft(y)

# plot the fft of the gausian function
plot(x, abs(Y))
show()

El resultado no es del todo correcto, porque la FFT de una función gaussiana debería ser una función gaussiana en sí misma ...


8
2018-03-22 21:52


origen


Respuestas:


np.fft.fft devuelve un resultado en el llamado "orden estándar": (de los documentos)

Si A = fft(a, n), entonces A[0]   contiene el término de frecuencia cero (el   media de la señal), que siempre es   puramente real para entradas reales. Entonces    A[1:n/2] contiene el   términos de frecuencia positiva, y    A[n/2+1:] contiene el   términos de frecuencia negativa, en orden de   frecuencia decreciente negativa.

La función np.fft.fftshift reorganiza el resultado en el orden que la mayoría de los humanos espera (y que es bueno para el trazado):

La rutina np.fft.fftshift(A)   desplaza transforma y su   frecuencias para poner la frecuencia cero   componentes en el medio ...

Entonces usando np.fft.fftshift:

import matplotlib.pyplot as plt
import numpy as np

N = 128
x = np.arange(-5,5,10./(2*N))
y = np.exp(-x*x)
y_fft = np.fft.fftshift(np.abs(np.fft.fft(y)))/ np.sqrt(2 * N)
plt.plot(x,y)
plt.plot(x,y_fft)
plt.show()

enter image description here


13
2018-03-22 22:59



Se muestra con el centro (es decir, la media) en el índice de coeficientes cero. Es por eso que parece que la mitad derecha está a la izquierda y viceversa.

EDITAR: explore el siguiente código:

import scipy
import scipy.signal as sig
import pylab
x = sig.gaussian(2048, 10)
X = scipy.absolute(scipy.fft(x))
pylab.plot(x)
pylab.plot(X)
pylab.plot(X[range(1024, 2048)+range(0, 1024)])

La última línea trazará X comenzando desde el centro del vector, luego envuelve al principio.


3
2018-03-22 21:54



Su resultado ni siquiera está cerca de un gaussiano, ni siquiera uno dividido en dos mitades.

Para obtener el resultado que espera, tendrá que posicionar su propio Gaussian con el centro en el índice 0, y el resultado también se posicionará de esa manera. Pruebe el siguiente código:

from pylab import *
N = 128
x = r_[arange(0, 5, 5./N), arange(-5, 0, 5./N)]
y = exp(-x*x)
y_fft = fft(y) / sqrt(2 * N)
plot(r_[y[N:], y[:N]])
plot(r_[y_fft[N:], y_fft[:N]])
show()

Los comandos de trazado dividen las matrices en dos partes y las intercambian para obtener una mejor imagen.

Plot


3
2018-03-22 22:05



Una transformada de Fourier se repite implícitamente indefinidamente, ya que es una transformación de una señal que se repite implícitamente indefinidamente. Tenga en cuenta que cuando pasa y ser transformado, el x los valores no se suministran, por lo tanto, de hecho, el gaussiano que se transforma se centra en el valor mediano entre 0 y 256, entonces 128.

Recuerde también que la traducción de f (x) es un cambio de fase de F (x).


2
2018-03-22 22:08



Siguiendo con la respuesta de Sven Marnach, una versión más simple sería esta:

from pylab import *
N = 128

x = ifftshift(arange(-5,5,5./N))

y = exp(-x*x)
y_fft = fft(y) / sqrt(2 * N)

plot(fftshift(y))
plot(fftshift(y_fft))

show()

Esto produce una trama idéntica a la anterior.

La clave (y esto me parece extraño) es que el supuesto de datos asumido por NumPy --- en ambos frecuencia y los dominios de tiempo --- es tener primero el valor "cero". Esto no es lo que esperaría de otras implementaciones de FFT, como las bibliotecas FFTW3 en C.

Esto fue un poco fudge en las respuestas de unutbu y Steve Tjoa arriba, porque están tomando el valor absoluto de la FFT antes de trazarlo, borrando así los problemas de fase resultantes de no usar el "orden estándar" en el tiempo.


2
2018-05-17 18:27