FFTW 库:验证傅立叶变换导数 属性
FFTW library: verifying fourier transform derivative property
我正在尝试验证与 fftw 库的这种关系:
因此,我选择f为高斯分布,计算其导数的傅里叶变换,并将其与高斯乘以ik的傅里叶变换进行比较。这是我得到的:
这很奇怪,尤其是因为高斯导数的傅里叶变换图(即红色的)在原点处不是 0,而它应该是(我检查了解析图)。
代码对我来说似乎没问题,无论如何它在这里(我正在使用 C):
int main() {
int i, N = 100;
double v[N], x[N], k[N/2+1], vd[N];
double dx = 2*pi/N, dk=2*pi/(N*dx), tmp;
fftw_complex *out;
fftw_plan forward,inverse;
out = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex )*( N/2 + 1 ));
forward = fftw_plan_dft_r2c_1d(N, v, out, FFTW_ESTIMATE);
inverse = fftw_plan_dft_c2r_1d(N, out, vd, FFTW_ESTIMATE);
//Initialise arrays
for( i = 0; i < N; i++ ) {
x[i] = dx*i;
v[i] = -2*x[i]*exp( -pow( x[i], 2) );
printf( " %le %le \n ", x[i], v[i] );
}
for( i = 0; i < N; i++ ) {
k[i]=i*dk;
}
k[N/2]=0.;
//Compute fft
fftw_execute( forward );
//Print the results
for( i = 0; i < N/2 + 1 ; i++ ) {
printf( "%le %le %le \n", i*dk, out[i][0], out[i][1] );
}
//Multiply by ik
for( i = 0; i < ( N/2 + 1 ); i++ ) {
tmp=out[i][0];
out[i][0]=-k[i]*out[i][1];
out[i][1]=k[i]*tmp;
printf( "%le %le %le \n", i*dk, out[i][0], out[i][1] );
}
fftw_destroy_plan(forward);
fftw_destroy_plan(inverse);
fftw_free(out);
return 0;
}
谁能告诉我我做错了什么?
导出的高斯离散化信号应该是:
x[i] = dx*i;
v[i] = -2*x[i]*exp( -pow( x[i], 2) );
然而,离散傅里叶变换对应于周期性信号的傅里叶变换。实际上,下划线离散化函数被写为正弦波的无限加权和。
因此,当应用 DFT 时,上面的离散化信号对应于高斯导数的周期化一半。实际上,它的平均值(零频率)不为零,因为所有值都是负数。
要模拟高斯导数(或 "finite" 范围内的任何其他信号),必须覆盖信号的整个范围。因此,必须选择 dx
使得 dx*N>>sigma
,其中 sigma
是高斯的标准偏差。并且必须覆盖函数的所有支持,包括导数的积极方面。
你能试试这样的东西吗:
double sigma=dx*N*0.1;
x[i] = dx*i-dx*(N/2);
v[i] = -2*x[i]*exp( -pow( x[i]/sigma, 2) );
由于标准偏差的值,必须留有缩放比例。
DFT 对 non-periodic 函数仍然有用,但是 non-periodic 函数要通过使用 window 映射到周期函数。这里观察到的是应用矩形 window、周期化然后推导与推导、应用矩形 window 最后周期化不同。虽然都是线性的,但这些运算符不会通勤!
我正在尝试验证与 fftw 库的这种关系:
因此,我选择f为高斯分布,计算其导数的傅里叶变换,并将其与高斯乘以ik的傅里叶变换进行比较。这是我得到的:
这很奇怪,尤其是因为高斯导数的傅里叶变换图(即红色的)在原点处不是 0,而它应该是(我检查了解析图)。
代码对我来说似乎没问题,无论如何它在这里(我正在使用 C):
int main() {
int i, N = 100;
double v[N], x[N], k[N/2+1], vd[N];
double dx = 2*pi/N, dk=2*pi/(N*dx), tmp;
fftw_complex *out;
fftw_plan forward,inverse;
out = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex )*( N/2 + 1 ));
forward = fftw_plan_dft_r2c_1d(N, v, out, FFTW_ESTIMATE);
inverse = fftw_plan_dft_c2r_1d(N, out, vd, FFTW_ESTIMATE);
//Initialise arrays
for( i = 0; i < N; i++ ) {
x[i] = dx*i;
v[i] = -2*x[i]*exp( -pow( x[i], 2) );
printf( " %le %le \n ", x[i], v[i] );
}
for( i = 0; i < N; i++ ) {
k[i]=i*dk;
}
k[N/2]=0.;
//Compute fft
fftw_execute( forward );
//Print the results
for( i = 0; i < N/2 + 1 ; i++ ) {
printf( "%le %le %le \n", i*dk, out[i][0], out[i][1] );
}
//Multiply by ik
for( i = 0; i < ( N/2 + 1 ); i++ ) {
tmp=out[i][0];
out[i][0]=-k[i]*out[i][1];
out[i][1]=k[i]*tmp;
printf( "%le %le %le \n", i*dk, out[i][0], out[i][1] );
}
fftw_destroy_plan(forward);
fftw_destroy_plan(inverse);
fftw_free(out);
return 0;
}
谁能告诉我我做错了什么?
导出的高斯离散化信号应该是:
x[i] = dx*i;
v[i] = -2*x[i]*exp( -pow( x[i], 2) );
然而,离散傅里叶变换对应于周期性信号的傅里叶变换。实际上,下划线离散化函数被写为正弦波的无限加权和。
因此,当应用 DFT 时,上面的离散化信号对应于高斯导数的周期化一半。实际上,它的平均值(零频率)不为零,因为所有值都是负数。
要模拟高斯导数(或 "finite" 范围内的任何其他信号),必须覆盖信号的整个范围。因此,必须选择 dx
使得 dx*N>>sigma
,其中 sigma
是高斯的标准偏差。并且必须覆盖函数的所有支持,包括导数的积极方面。
你能试试这样的东西吗:
double sigma=dx*N*0.1;
x[i] = dx*i-dx*(N/2);
v[i] = -2*x[i]*exp( -pow( x[i]/sigma, 2) );
由于标准偏差的值,必须留有缩放比例。
DFT 对 non-periodic 函数仍然有用,但是 non-periodic 函数要通过使用 window 映射到周期函数。这里观察到的是应用矩形 window、周期化然后推导与推导、应用矩形 window 最后周期化不同。虽然都是线性的,但这些运算符不会通勤!