Aproximación numérica de la transformada continua de Fourier

Question

Dada una función $F(k)$ en el espacio de frecuencias (suficientemente agradable, por ejemplo, una gaussiana), me gustaría calcular su inversa de Fourier \begin {Edición}f(x) = \frac {1}{2 \pi } \int_ {- \infty }^{ \infty }e^{ikx}F(k)dk \end {ecuación} numéricamente (no hay una fórmula analítica explícita), en algunos puntos específicos distribuidos uniformemente $x_n$ .

Supongamos que $F(k)$ es simétrica respecto a $k=0$ y "esencialmente cero" fuera del intervalo $-a/2\leq k \leq a/2$ .

Ingenuamente, una forma de proceder sería dividir el intervalo $(-a/2,a/2)$ en $N$ subintervalos (elija $N$ para que sea una potencia de dos) y luego aproximar a trozos la integral: $$ f(x_n) \approx \frac{1}{2\pi}\frac{a}{N}\sum_{m=0}^{N-1} e^{ik_mx_n}F(k_m)$$ donde, por ejemplo, podríamos tomar $k_m = (m-N/2)\frac{a}{N}$ para $0\leq m\leq N-1$ .

Por eficiencia, me gustaría utilizar la transformada rápida de Fourier (FFT) para aproximar $f(x_n)$ . Así que tengo que poner $f(x)$ en el formato correcto para la Transformada Discreta de Fourier (DFT) $$L_k = \frac{1}{N}\sum_{m=0}^{N-1}G_me^{(2\pi i) km/N},\,\,\,\,k = 0,\ldots,N-1$$ donde $G_m$ es el $m$ de una función determinada $G$ .

Para proceder, dividimos el intervalo $(-a/2,a/2)$ en $N$ piezas ( $N$ una potencia de dos) y que $k_m = (m-N/2)\frac{a}{N}$ para $0\leq m\leq N-1$ como antes. Entonces $$f(x) \approx \frac{1}{2\pi}\frac{a}{N}\sum_{m=0}^{N-1} e^{ik_mx}F(k_m)$$ Sin embargo, para que esta fórmula tenga el formato correcto para utilizar la DFT, no podemos evaluar f en cualquier $x$ : Creo que deben ser de la forma $x_k = \frac{2\pi}{a}(k-N/2), 0\leq k \leq N-1$ .

Así que para utilizar la FFT, los puntos en los que puedo "evaluar" la función $f$ está determinada por la forma en que muestreo mi función $F(k)$ y el número de puntos de muestreo.

Pero en el método ingenuo, puedo, en principio, encontrar $f(x)$ para cualquier x dada, independientemente del número de puntos de muestra que tome de la función $F(k)$ .

Como puedes ver, mis conocimientos numéricos al respecto son, en el mejor de los casos, escasos.

Mis preguntas son:

1. ¿Es posible utilizar la FFT, si necesito calcular f en puntos específicos $x_n$ ?

2. Si no es así, ¿hay otros métodos numéricos que sean más adecuados?

Answer 1

Los siguientes dos artículos pueden ser de utilidad. No permiten hacer la transformación en cualquier valor de x. Pero sí permiten hacer la transformación desde m puntos igualmente espaciados en el dominio w a m puntos igualmente espaciados en el dominio x. El artículo de Bailey y Swarztrauber también utiliza una función gaussiana como la tuya como ejemplo. Siento no poder ser más específico pero estoy aprendiendo todo esto por mi cuenta.

Bailey, D., y P. Swarztrauber. 1994. A Fast Method for the Numerical Evaluation of Continuous Fourier and Laplace Transforms'. SIAM Journal on Scientific Computing 15 (5): 1105-10. doi:10.1137/0915067.

Inverarity, G. 2002. Fast Computation of Multidimensional Fourier Integrals'. SIAM Journal on Scientific Computing 24 (2): 645-51. doi:10.1137/S106482750138647X.

Answer 2

Una forma que se me ocurre sería calcular ambos aproximaciones para $f(x_n)$ y para $f'(x_n)$ un poco de rejilla $x_n$ que permite utilizar la FFT. Para calcular las aproximaciones de $f'(x_n)$ se puede utilizar la identidad $\mathcal{F}(f')(k) = ik\mathcal{F}(f)(k)$ (es decir, se obtiene la transformada de Fourier de la derivada multiplicando por $ik$ ).

Para evaluar $f(x)$ en $y$ , entonces encontrarías el más cercano $x_n,x_{n+1}$ con $x_n \leq y < x_{n+1}$ y utilizar la interpolación polinómica (con polinomios de grado 3) para encontrar $f(y)$ de $f(x_n)$ , $f(x_{n+1})$ , $f'(x_n)$ , $f'(x_{n+1})$ .