作为 this question 的扩展我问的。实高斯的傅立叶变换是实高斯。当然,仅类似于高斯分布的一组点的 DFT 并不总是完美的高斯分布,但它肯定应该接近。在下面的代码中,我使用 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();
}
最佳答案
您的虚部结果是正确且符合预期的。
与零 (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) 为实数且偶数时,虚部(否则为奇数)才会在积分过程中消失。但这没有什么意义,因为实部和虚部仅在特殊情况下才具有物理意义。
直接的物理意义在于abs值(幅度谱
),即abs。值的平方(强度谱
)和相位或角度(相位谱
)。
如果虚数部分不以时间向量的中心为中心,则会发生更明显的零偏移。尝试将 x
向量移动 dx
的一部分。
请参阅下文,输入的移位 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))
关于fft - GSL 快速傅立叶变换 - 变换高斯的非零虚数?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/37015321/