是否可以在 C++ 中使用 FFTW 对我的数据集执行 FFT？

问：

我正在编写一个程序来计算一组复杂矩阵的傅里叶变换，这些复矩阵位于 3D k 空间中沿一条线的 N 个离散点处，从 k=（0， 0， 1/2） 到 k=（1/2， 1/3， 1/2）。

这是我试图计算的公式。

这是我目前的工作实现，使用 Eigen C++ 包进行矩阵运算： （k_coordinates 是 3xN 向量，R_coordinates 是 3x1155 向量）

// Perform Fourier transform on H(R) to get H(k)
    Matrix18 H_k_data[k_coordinates.size()]; // 2*N array of 18x18 matrices
    for (int k_index=0; k_index<k_coordinates.size(); k_index++) { // Loop over k points
        for (int R=0; R<1155; R++) { // Loop over R points
            // Calculate dot product k.R
            double k_dot_R = 0;
            for (int i=0; i<3; i++) {
                k_dot_R += (k_coordinates[k_index][i] * R_coordinates_data[R][i] * 2 * M_PI);
            }
            std::complex<double> phase_factor = std::polar(1.0, -k_dot_R); // Compute phase factor as e^(-ik.r)

            // Add final value to H(k) at each index
            H_k_data[k_index] += (H_R_data[R] * phase_factor);
        }
    }

这有效但非常慢，所以我正在尝试使用 FFTW C++ 库实现 FFT 算法，但是我无法弄清楚如何正确设置我的数据以使用它。我相信我需要使用 fftw_plan_dft_3d（） 函数（3D 离散傅里叶变换），但看不到我如何实际将数据输入到该函数中。

从本质上讲，我需要的是k_coordinates中每个点 k=（kx， ky， kz），以计算 （H（R） * e^（-2pi * （kxRx + kyRy + kzRz）） 的所有 R 坐标的总和。

下面是上下文的完整代码：

#include<iostream>
#include<iomanip>
#include<fstream>
#include<vector>
#include<string>
#include<complex>
#include<Eigen/Dense>
#include<math.h>
#include<chrono>
#include<fftw3.h>

using namespace Eigen;
typedef Matrix<std::complex<double>, 18, 18> Matrix18;

template<typename T>
std::vector<double> linspace(T start_in, T end_in, int num_in)
{
    std::vector<double> linspaced;

    double start = static_cast<double>(start_in);
    double end = static_cast<double>(end_in);
    double num = static_cast<double>(num_in);

    if (num == 0) {
        return linspaced;
    }
    if (num == 1) {
        linspaced.push_back(start);
        return linspaced;
    }

    double delta = (end - start) / (num - 1);

    for(int i=0; i < num-1; ++i) {
        linspaced.push_back(start + delta * i);
    }
    linspaced.push_back(end);
    return linspaced;
}

std::vector<Matrix18> convert_to_eigen(std::vector<std::vector<std::vector<std::complex<double>>>> data)
{
    std::vector<Matrix18> data_eigen;
    for (int k=0; k<data.size(); k++){
        Matrix18 H;
        for (int i=0; i<18; i++) {
            for (int j=0; j<18; j++) {
                H(i,j) = data[k][i][j];
            }
        }
        data_eigen.push_back(H);
    }
    return data_eigen;
}

void print_message(std::string message, std::chrono::steady_clock::time_point starting_time)
{
    std::chrono::steady_clock::time_point end_convert = std::chrono::steady_clock::now(); 
    std::cout << message << "(" << std::chrono::duration_cast<std::chrono::seconds>(end_convert - starting_time).count() << "s)" << std::endl;
}

std::vector<std::vector<double>> get_k_coordinates()
{
    // Create vector of k coordinates 
    int num_slices = 1000;
    std::vector<std::vector<double>> k_coordinates; // Coordinates of all points along line
    std::vector<double> kx = linspace(0.0, 0.5, num_slices);
    std::vector<double> ky = linspace(0.0, 1.0/3.0, num_slices);
    std::vector<double> kz = linspace(0.5, 0.5, num_slices);

    for (int i=0; i<num_slices; i++) { // Points along L-A
        std::vector<double> slice_coords; // 3x1 vector containing the coordinates of one point along L-A
        slice_coords.push_back(kx[i]);
        slice_coords.push_back(ky[i]);
        slice_coords.push_back(kz[i]);
        k_coordinates.push_back(slice_coords);
    }
    return k_coordinates;
}


int main()
{
    std::chrono::steady_clock::time_point starting_time = std::chrono::steady_clock::now();

    print_message("Reading in data...", starting_time);

    std::string filename = "BiTeI_hr_topological.dat";
    std::ifstream data_stream;
    data_stream.open(filename);
    std::string line;
    
    // Read in the R weights data
    std::vector<int> R_weights_data; // 1D Vector to store weights of lattice vectors (1155)
    int lines_to_skip = 3;
    int count = 0;
    while ((count < 77) && std::getline(data_stream, line)) {
        if (lines_to_skip > 0) { // Skip first 3 lines
            lines_to_skip--;
        } else {
            std::istringstream iss(line);
            int num;
            while (iss >> num) {
                R_weights_data.push_back(num);
            }
            count++;
        }
    }

    // Read in hamiltonians and corresponding R coordinate data
    std::vector<std::vector<int>> R_coordinates_data; // 2D Vector to store R coefficients l,m,n (3x1155)
    std::vector<std::vector<std::vector<std::complex<double>>>> H_R_data; // 3D vector to store hamiltonians H(R)
    std::vector<std::vector<std::complex<double>>> hamiltonian; // 2D vector to store complex hopping parameters
    std::vector<std::complex<double>> hamiltonian_line; // 1D vector to store single line in hamiltonian
    int line_count = 0;
    while (std::getline(data_stream, line)) {
        std::istringstream iss(line);
        int num;
        
        std::vector<int> R_coordinates_line;
        if (line_count % 324 == 0) {
            
            for (int i=0; i<3; i++) {
                iss >> num;
                R_coordinates_line.push_back(num);
            }
            R_coordinates_data.push_back(R_coordinates_line);
        } else {
            for (int i=0; i<3; i++) { // Still skip coefficient characters even when not needed
                iss >> num;
            }
        }
        for (int i=0; i<2; i++){ // Skip index characters
            iss >> num;
        }

        // Read in hamiltonian matrix values
        double real, imag;
        iss >> real;
        iss >> imag;
        std::complex<double> matrix_value(real, imag);
        if (hamiltonian_line.size() == 18) {
            hamiltonian.push_back(hamiltonian_line);
            hamiltonian_line.clear();
        }
        hamiltonian_line.push_back(matrix_value);
        if (hamiltonian.size() == 18) {
            H_R_data.push_back(hamiltonian);
            hamiltonian.clear();
        }

        line_count++;
    }
    hamiltonian.push_back(hamiltonian_line);
    hamiltonian_data.push_back(hamiltonian);

    std::vector<std::vector<double>> k_coordinates = get_k_coordinates(); // Create set of k coordinate vectors for which to perform fourier transform

    // Convert H_R data into Eigen matrices
    std::vector<Matrix18> H_R_data_eigen; // 1155 length vector of 18x18 matrices
    for (int k=0; k<H_R_data.size(); k++){
        Matrix18 H;
        for (int i=0; i<18; i++) {
            for (int j=0; j<18; j++) {
                H(i,j) = H_R_data[k][i][j];
            }
        }
        H_R_data_eigen.push_back(H);
    }

    print_message("Data processed. Performing Fourier transform...", starting_time);

    // Perform Fourier transform on H(R) to get H(k)
    Matrix18 H_k_data[k_coordinates.size()]; // 2*N array of 18x18 matrices
    for (int k_index=0; k_index<k_coordinates.size(); k_index++) { // Loop over k points
        for (int R=0; R<1155; R++) { // Loop over R points
            // Calculate dot product k.R
            double k_dot_R = 0;
            for (int i=0; i<3; i++) {
                k_dot_R += (k_coordinates[k_index][i] * R_coordinates_data[R][i] * 2 * M_PI);
            }
            std::complex<double> phase_factor = std::polar(1.0, -k_dot_R); // Compute phase factor as e^(-ik.r)

            // Add final value to H(k) at each index
            H_k_data[k_index] += (H_R_data[R] * phase_factor);
        }
    }
    
    print_message("Done.", starting_time);
    
    return 0;
}

我尝试在我的数据集上实现 FFT 算法，但我无法弄清楚如何操作我的数据以在 FFTW 函数中使用。