GSL高速フーリエ変換 - 変換されたガウスの非ゼロ虚数？

Question

私が尋ねたthis questionの拡張です。実際のガウス分布のフーリエ変換は実際のガウス分布である。もちろん、Gaussianに似ている点集合のDFTは必ずしも完全なGaussianではありませんが、確かに近いはずです。下のコードでは、GSLを使ってこの[離散]フーリエ変換を行っています。返された/変換された実際のコンポーネント（リンクされた質問で概説されている）の問題とは別に、私は虚数コンポーネント（これは全く同じでなければならない）に対して奇妙な結果を得ています。確かに、それは規模は非常に小さいですが、それはまだ変です。 この非対称の原因は何ですか？&ファンキーな出力？

#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

虚部のためのあなたの結果が正しいと期待されています。

ゼロ（10^-15）との差は、pi（12桁、piはFFTで使用されますが、私はあなたの中のpiをオーバーライドしているかどうかを知ることができません。ルーチン）。

実関数のFFTは一般的に実関数ではありません。あなたが解析的に計算を行うときは、次の式を積分：

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

ので、関数f（t）がreal and evenである場合にのみ（それ以外の奇数である）、虚部は​​、統合時に消えます。実数部と虚数部は特別な場合にのみ物理的な意味を持つため、意味がありません。

直接的な物理的意味は、abs値（magnitude spectrum）、abs。値の二乗（intensity spectrum）と位相または角度（phase spectrum）です。

虚数部のゼロからのより重要なオフセットは、時間ベクトルの中央に中央揃えされていないと発生します。 xベクターをdxの何分の1かだけシフトしてみてください。

dx/2（右の列）による入力のシフトが虚数部にどのように影響するかを示します（例：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))