GSL-Fast-Fourier-Transformation - Nicht-Null-Imaginär für transformierte Gauß-Transformation?

Question

Als eine Erweiterung zu this question, die ich fragte. Die Fourier-Transformation eines realen Gaußschen ist eine echte Gaußsche. Natürlich wird eine DFT einer Menge von Punkten, die nur einer Gaußschen ähneln, nicht immer eine perfekte Gauß'sche sein, aber sie sollte sicher nahe sein. Im folgenden Code nehme ich diese [diskrete] Fourier-Transformation unter Verwendung von GSL. Abgesehen von der Frage der zurückgegebenen/transformierten realen Komponenten (die in der verknüpften Frage umrissen sind), bekomme ich ein merkwürdiges Ergebnis für die imaginäre Komponente (die identisch Null sein sollte). Zugegeben, es ist sehr klein in der Größenordnung, aber es ist immer noch seltsam. Was ist die Ursache für diesen asymmetrischen & funky Ausgang?

#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

Ihr Ergebnis für den imaginären Teil richtig ist und zu erwarten.

Der Unterschied zu Null (10^-15) ist weniger als Genauigkeit, die Sie Pi geben (12 Ziffern, Pi wird in der FFT verwendet, aber ich kann nicht wissen, ob Sie die Pi innerhalb der Routine).

Die FFT einer realen Funktion ist im Allgemeinen keine reelle Funktion. Wenn Sie die Mathematik analytisch integrieren Sie über den folgenden Ausdruck:

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

so nur, wenn die Funktion f (t) ist real and even wird der imaginäre Teil (die sonst ungerade ist) bei der Integration verschwinden. Dies hat jedoch wenig Bedeutung, da der Realteil und der Imaginärteil nur in speziellen Fällen physikalische Bedeutung haben.

Direkte physikalische Bedeutung ist in der Abs-Wert (magnitude spectrum), die abs. Wert im Quadrat (intensity spectrum) und die Phase oder Winkel (phase spectrum).

Ein signifikanterer Versatz von Null im Imaginärteil würde auftreten, wenn er nicht in der Mitte des Zeitvektors zentriert wäre. Versuchen Sie, den Vektor x um einen Bruchteil von dx zu verschieben.

Siehe unten, wie sich die Verschiebung der Eingabe um dx/2 (rechte Spalte) auf den Imaginärteil auswirkt, aber nicht auf die Größe (Beispiel in Python, Numpy geschrieben).

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))