问题描述
目标
我有表示 实 标量场 f(x,y) 的二维数据，我想计算 fourier (q) space 中的一维频谱。
我的项目主要是用 C++ 编写的，但有一些代码块是用 C 语言编写的。我也在使用 FFTW3 库并执行就地转换。
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <vector>
double PI = 4*atan(1.0);
int WrapAround(int k,int N)
{
if (k <= N/2)
return k;
else
return k-N;
}
int main(int argc,char* argv[])
{
if (argc <= 1 or argc >= 3)
{
std::cout << "One argument is needed: degrees of freedom per dimension. Preferably power of 2. Exiting..." << std::endl;
exit(-1);
}
/* Generate x and y axis */
/* ----------------------------------------------------------- */
int N = strtol(argv[1],nullptr,0);
double Lx = 2*PI; //length of Box
double Ly = Lx;
double dx = Lx/N;
double dy = Ly/N;
std::vector<double> x,y;
for (int i=0; i<N; i++)
x.push_back(i*dx);
y = x;
/* ----------------------------------------------------------- */
/* Define data for which the power spectrum will be calculated */
/* ----------------------------------------------------------- */
double* field = new double[N*N]{};
std::cout << "Field is: " << std::endl;
for (int j=0; j<N; j++)
{
for (int i=0; i<N; i++)
{
/* test different cases here */
//field[ i + N*j ] = 1.;
field[ i + N*j ] = 4.0 * sin(3.0*x[i]);
//field[ i + N*j ] = 4.0 * cos(1.0*y[j]);
//field[ i + N*j ] = 5.0 * sin(-x[i] + 3.0*y[j]) + 3.0 * sin(x[i] + 3.0*y[j]);
//field[ i + N*j ] = 2.0 * sin( 3.0*(2*PI*x[i]/Lx) + 3.0*(2*PI*y[j]/Ly) );
std::cout << field[ i + N*j ] << "\t";
}
std::cout << std::endl;
}std::cout << std::endl;
/* ----------------------------------------------------------- */
/* Prepare in-place fft */
/* ----------------------------------------------------------- */
int blockSize = N*(N+2);
fftw_plan x2k,k2x;
double* IN = (double*) fftw_malloc( sizeof(double)*blockSize );
fftw_complex* OUT = (fftw_complex*) IN;
x2k = fftw_plan_dft_r2c_2d(N,N,IN,OUT,FFTW_MEASURE);
k2x = fftw_plan_dft_c2r_2d(N,FFTW_MEASURE);
/* ----------------------------------------------------------- */
/* Fill the allocated block pointed by the pointer "in" with the field data */
/* ----------------------------------------------------------- */
for (int j=0; j<N; j++){
for(int i=0; i<N; i++)
IN[ i + j*(N+2) ] = field[ i + N*j ]; //note the (N+2) term
}
/* ----------------------------------------------------------- */
/* Execute FFT */
fftw_execute(x2k);
/* Fill a buffer with the power spectrum values */
/* ----------------------------------------------------------- */
double* powSpec = new double[N*(N/2+1)]{};
double Re,Im,val,q1,q2;
std::cout << "Power spectrum values: " << std::endl;
for (int k2=0; k2<N; k2++){
q2 = WrapAround(k2,N);
for (int k1=0; k1<N/2+1; k1++){
q1 = k1;
Re = IN[ (N+2)*k2 + 2*k1 ];
Im = IN[ (N+2)*k2 + 2*k1+1 ];
val = (Re*Re + Im*Im) / pow(N,4); // normalize
powSpec[ (N/2+1)*k2 + k1 ] = val;
// Now print the buffer
if (val < 1e-16)
val = 0.0;
std::cout << "(" << q1 << "," << q2 << ")\t";
std::cout << "(" << k1 << "," << k2 << ")\t" << val << "\t";
}
std::cout << std::endl;
}std::cout << std::endl;
/* ----------------------------------------------------------- */
/* Execute inverse FFT */
fftw_execute(k2x);
/* Print result of inverse FFT */
/* ----------------------------------------------------------- */
std::cout << "After applying inverse FFT: " << std::endl;
for (int k2=0; k2<N; k2++){
for (int k1=0; k1<N; k1++){
std::cout << IN[ (N+2)*k2 + k1 ] << "\t";
}
std::cout << std::endl;
}std::cout << std::endl;
/* ----------------------------------------------------------- */
/* Radial averaging */
/* ----------------------------------------------------------- */
double dk = 2*PI/Lx;
double factor = pow(dk,2);
int qdiagMax = floor( sqrt( pow(dk*N/2,2) + pow(dk*N/2,2) ) ) + 1;
int qbin;
double qSquared;
double* S1D = new double[qdiagMax]{}; // buffer to contain 1D spectrum
double* DoS = new double[qdiagMax]{}; // buffer to contain the density of states
for (int k2=0; k2<N; k2++){
q2 = WrapAround(k2,N);
for (int k1=0; k1<N/2; k1++){
q1 = WrapAround(k1,N);
qSquared = factor * (q1*q1 + q2*q2);
qbin = floor(sqrt(qSquared));
S1D[qbin] += powSpec[ N/2*k2 + k1 ];
DoS[qbin] ++;
}
}
DoS[0] = 1; //the first grid point on its own bin
DoS[qdiagMax-1] = 1; // the last grid point on its own bin
//the two lines above are wrong and inner loop must go up to N/2+1. (EDIT)
for (int i=0; i<qdiagMax; i++)
S1D[i] /= DoS[i];
/* ----------------------------------------------------------- */
/* Print results of radial averaging */
/* ----------------------------------------------------------- */
double printOut;
std::cout << "Radially averaged power spectrum: " << std::endl;
for (int i=0; i<qdiagMax; i++)
{
printOut = S1D[i];
if ( S1D[i] < 1e-16 ) printOut = 0;
std::cout << printOut << "\t";
}std::cout << "\n" <<std::endl;
std::cout << "Density of states: " << std::endl;
for (int i=0; i<qdiagMax; i++)
std::cout << DoS[i] << "\t";
std::cout << std::endl;
/* ----------------------------------------------------------- */
/* Destroy plan & release resources */
fftw_destroy_plan(x2k);
fftw_free(IN);
OUT = nullptr;
IN = nullptr;
return 0;
}
让我们看看 field[ i + N*j ] = 4.0 * sin(3.0*x[i]); 的情况。
我们期望位置 (3,0) 中有一个非零值。注意第一个括号
引用傅立叶空间中的 q 值，而 k 只是索引。 y
方向被缠绕。让我们执行 ./a.out 8
Field is:
0 2.82843 -4 2.82843 1.46958e-15 -2.82843 4 -2.82843
0 2.82843 -4 2.82843 1.46958e-15 -2.82843 4 -2.82843
0 2.82843 -4 2.82843 1.46958e-15 -2.82843 4 -2.82843
0 2.82843 -4 2.82843 1.46958e-15 -2.82843 4 -2.82843
0 2.82843 -4 2.82843 1.46958e-15 -2.82843 4 -2.82843
0 2.82843 -4 2.82843 1.46958e-15 -2.82843 4 -2.82843
0 2.82843 -4 2.82843 1.46958e-15 -2.82843 4 -2.82843
0 2.82843 -4 2.82843 1.46958e-15 -2.82843 4 -2.82843
Power spectrum values:
(0,0) (0,0) 0 (1,0) (1,0) 0 (2,0) (2,0) 0 (3,0) (3,0) 4 (4,0) (4,0) 0
(0,1) (0,1) 0 (1,1) (1,1) 0 (2,1) (2,1) 0 (3,1) (3,1) 0 (4,1) (4,1) 0
(0,2) (0,2) 0 (1,2) (1,2) 0 (2,2) (2,2) 0 (3,2) (3,2) 0 (4,2) (4,2) 0
(0,3) (0,3) 0 (1,3) (1,3) 0 (2,3) (2,3) 0 (3,3) (3,3) 0 (4,3) (4,3) 0
(0,4) (0,4) 0 (1,4) (1,4) 0 (2,4) (2,4) 0 (3,4) (3,4) 0 (4,4) (4,4) 0
(0,-3) (0,5) 0 (1,-3) (1,5) 0 (2,-3) (2,5) 0 (3,-3) (3,5) 0 (4,-3) (4,5) 0
(0,-2) (0,6) 0 (1,-2) (1,6) 0 (2,-2) (2,6) 0 (3,-2) (3,6) 0 (4,-2) (4,6) 0
(0,-1) (0,7) 0 (1,-1) (1,7) 0 (2,-1) (2,7) 0 (3,-1) (3,7) 0 (4,-1) (4,7) 0
After applying inverse FFT:
6.31089e-30 181.019 -256 181.019 9.40529e-14 -181.019 256 -181.019
6.31089e-30 181.019 -256 181.019 9.40529e-14 -181.019 256 -181.019
6.31089e-30 181.019 -256 181.019 9.40529e-14 -181.019 256 -181.019
6.31089e-30 181.019 -256 181.019 9.40529e-14 -181.019 256 -181.019
6.31089e-30 181.019 -256 181.019 9.40529e-14 -181.019 256 -181.019
6.31089e-30 181.019 -256 181.019 9.40529e-14 -181.019 256 -181.019
6.31089e-30 181.019 -256 181.019 9.40529e-14 -181.019 256 -181.019
6.31089e-30 181.019 -256 181.019 9.40529e-14 -181.019 256 -181.019
Radially averaged power spectrum:
0 0 0 0.363636 0 0
Density of states:
1 5 9 11 5 1
我得到的值是 4，而我希望有 16，因为功率是幅度的平方。我必须乘以 4 倍吗？它是否与 1D 情况类似，例如在时域中没有负频率并且由于 exp(-it) 项而我们乘以 2？
径向平均一维光谱
我的问题是它是否正确实施，因为我担心我遗漏了一个因素。例如，由于变换是就地的，我只有第一和第四象限 傅立叶空间。最后我必须乘以二吗？频率似乎出现在正确的 bin 中，但值是否正确？
代码很长，但我认为足够干净，所以如果有人在信号分析方面有经验可以花一些时间来帮助正确实现这个一维频谱，我将不胜感激。
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）