GSL Fast-Fourier Transform - Imaginaire non nul pour Gaussian transformé?

Question

Comme une extension à this question que j'ai demandé. La transformée de Fourier d'un vrai Gaussien est un vrai Gaussien. Maintenant, bien sûr, une DFT d'un ensemble de points qui ne ressemblent qu'à un Gaussien ne sera pas toujours un Gaussien parfait, mais il devrait certainement être proche. Dans le code ci-dessous, je prends cette transformation de Fourier [discrète] en utilisant GSL. Mis à part le problème des composants réels retournés/transformés (décrits dans la question liée), j'obtiens un résultat bizarre pour le composant imaginaire (qui devrait être identique à zéro). Certes, c'est de très petite taille, mais c'est quand même bizarre. Quelle est la cause de cette sortie asymétrique &?

#include <gsl/gsl_fft_complex.h> 
#include <gsl/gsl_errno.h> 
#include <fstream> 
#include <iostream> 
#include <iomanip> 

#define REAL(z,i) ((z)[2*(i)]) //complex arrays stored as [Re(z0),Im(z0),Re(z1),Im(z1),...] 
#define IMAG(z,i) ((z)[2*(i)+1]) 
#define MODU(z,i) ((z)[2*(i)])*((z)[2*(i)])+((z)[2*(i)+1])*((z)[2*(i)+1]) 
#define PI 3.14159265359 

using namespace std; 

int main(){ 

    int n = pow(2,9); 
    double data[2*n]; 
    double N = (double) n; 

    ofstream file_out("out.txt"); 

    double xmin=-10.; 
    double xmax=10.; 
    double dx=(xmax-xmin)/N; 
    double x=xmin; 

    for (int i=0; i<n; ++i){ 
     REAL(data,i)=exp(-100.*x*x); 
     IMAG(data,i)=0.; 
     x+=dx; 
    } 

    gsl_fft_complex_radix2_forward(data, 1, n); 

    for (int i=0; i<n; ++i){ 
     file_out<<(i-n/2)<<" "<<IMAG(data,((i+n/2)%n))<<'\n'; 
    } 

    file_out.close(); 
} 

Answer 1

Votre résultat pour la partie imaginaire est correcte et attendue. La différence à zéro (10^-15) est inférieure à la précision que vous donnez à pi (12 chiffres, pi est utilisé dans la FFT, mais je ne peux pas savoir si vous remplacez le pi dans le routine).

La FFT d'une fonction réelle n'est généralement pas une fonction réelle. Lorsque vous faites le calcul analytique vous intégrez sur l'expression suivante:

f(t) e^{i w t} = f(t) cos wt + i f(t) sin wt, 

donc que si la fonction f (t) est real and even sera la partie imaginaire (qui est par ailleurs étrange) disparaissent lors de l'intégration. Cela a peu de sens cependant, puisque la partie réelle et la partie imaginaire n'ont de sens physique que dans des cas particuliers.

La signification physique directe est dans la valeur abs (magnitude spectrum), l'abs. valeur au carré (intensity spectrum) et la phase ou l'angle().

Un décalage plus important de zéro dans la partie imaginaire se produirait s'il n'était pas centré au centre de votre vecteur de temps. Essayez de déplacer le vecteur x par une fraction de dx. Voir ci-dessous comment le décalage de l'entrée par dx/2 (colonne de droite) affecte la partie imaginaire, mais pas la grandeur (exemple écrit en Python, Numpy).

from __future__ import division 
import numpy as np 
import matplotlib.pyplot as p 
%matplotlib inline 

n=512 # number of samples 2**9 
x0,x1=-10,10 
dx=(x1-x0)/n 

x= np.arange(-10,10,dx) # even number, asymmetric range [-10, 10-dx] 

#make signal 
s1= np.exp(-100*x**2) 
s2= np.exp(-100*(x+dx/2)**2) 

#make ffts 
f1=np.fft.fftshift(np.fft.fft(s1)) 
f2=np.fft.fftshift(np.fft.fft(s2)) 

#plots 
p.figure(figsize=(16,12)) 
p.subplot(421) 
p.title('gaussian (just ctr shown)') 
p.plot(s1[250:262]) 
p.subplot(422) 
p.title('same, shifted by dx/2') 
p.plot(s2[250:262]) 

p.subplot(423) 
p.plot(np.imag(f1)) 
p.title('imaginary part of FFT') 
p.subplot(424) 
p.plot(np.imag(f2)) 

p.subplot(425) 
p.plot(np.real(f1)) 
p.title('real part of FFT') 
p.subplot(426) 
p.plot(np.real(f2)) 

p.subplot(427) 
p.plot(np.abs(f1)) 
p.title('abs. value of FFT') 
p.subplot(428) 
p.plot(np.abs(f2))