fft - 我可以使用 FFTW 沿 3D 数组的第一个和最后一个轴进行变换吗？

标签 fft fftw

使用 fftw_plan_many_dft 我可以沿 x,y 和 y,z 轴进行变换:

vector<complex<float_type>> yz_fft(vector<complex<float_type>> input, int N_X, int N_Y, int N_Z){
    vector<complex<float_type>> result(input.size());
    int rank = 2;
    int n[] = {N_Y,N_Z};
    int *inembed = n;
    int *onembed = n;
    int istride = 1;
    int ostride = 1;
    int idist = N_Y*N_Z;
    int odist = N_Y*N_Z;
    int howmany = N_X;
    fftw_plan plan = fftw_plan_many_dft(
            rank,
            n,
            howmany,
            reinterpret_cast<fftw_complex *>(input.data()),
            inembed,
            istride,
            idist,
            reinterpret_cast<fftw_complex *>(result.data()),
            onembed,
            ostride,
            odist,
            FFTW_FORWARD,
            FFTW_ESTIMATE);
    fftw_execute(plan);
    return result;
}

vector<complex<float_type>> xy_fft(vector<complex<float_type>> input, int N_X, int N_Y, int N_Z){
    vector<complex<float_type>> result(input.size());
    int rank = 2;
    int n[] = {N_X,N_Y};
    int *inembed = n;
    int *onembed = n;
    int istride = N_Z;
    int ostride = N_Z;
    int idist = 1;
    int odist = 1;
    int howmany = N_Z;
    fftw_plan plan = fftw_plan_many_dft(
            rank,
            n,
            howmany,
            reinterpret_cast<fftw_complex *>(input.data()),
            inembed,
            istride,
            idist,
            reinterpret_cast<fftw_complex *>(result.data()),
            onembed,
            ostride,
            odist,
            FFTW_FORWARD,
            FFTW_ESTIMATE);
    fftw_execute(plan);
    return result;
}

但我不知道如何进行 x,z 变换。我该怎么做？

最佳答案

所以有一种方法可以使用fftw_plan_many_dft来进行xz变换。反对票可能表明人们对此不感兴趣，但我还是决定分享它。如需解决方案，请检查下面的struct xz_fft_many。

#include <iostream>
#include <numeric>
#include <complex>
#include <fftw3.h>

#include <benchmark/benchmark.h>


using namespace std;

using float_type = double;
using index_type = unsigned long;

vector<complex<float_type>> get_data(index_type N){

    std::vector<complex<float_type>> data(N);
    iota(data.begin(), data.end(),0);

    return data;
}

void print(vector<complex<float_type>> data,index_type N_X,index_type N_Y,index_type N_Z){
    for(int i=0; i!=N_X; ++i){
        for(int j=0; j!=N_Y; ++j){
            for(int k=0; k!=N_Z; ++k){
                index_type idx = i*(N_Y*N_Z)+j*N_Z+k;
                cout<<"[ "<<i<<", "<<j<<", "<<k<<" ] = "<<data.data()[idx]<<endl;
            }
        }
    }
}

struct x_fft {
    vector<complex<float_type>>& data;
    vector<complex<float_type>> result;
    fftw_plan fft_plan;
    index_type N_X;
    index_type N_Y;
    index_type N_Z;


    x_fft(vector<complex<float_type>>& data,index_type N_X,index_type N_Y,index_type N_Z)
            : data(data), N_X(N_X), N_Y(N_Y), N_Z(N_Z)
    {
        result = vector<complex<float_type>>(data.size());
        int rank = 1;
        int n[] = {(int)N_X};
        int *inembed = n;
        int *onembed = n;
        int istride = N_Y*N_Z;
        int ostride = istride;
        int idist = 1;
        int odist = idist;
        int howmany = N_Y*N_Z;
        fft_plan = fftw_plan_many_dft(
                rank,
                n,
                howmany,
                reinterpret_cast<fftw_complex *>(data.data()),
                inembed,
                istride,
                idist,
                reinterpret_cast<fftw_complex *>(result.data()),
                onembed,
                ostride,
                odist,
                FFTW_FORWARD,
                FFTW_MEASURE);
    }

    const vector<complex<float_type>> &getResult() const {
        return result;
    }

    vector<complex<float_type>>& run(){
        fftw_execute(fft_plan);
        return result;
    }

};

struct z_fft {
    vector<complex<float_type>>& data;
    vector<complex<float_type>> result;
    fftw_plan fft_plan;
    index_type N_X;
    index_type N_Y;
    index_type N_Z;


    z_fft(vector<complex<float_type>>& data,index_type N_X,index_type N_Y,index_type N_Z)
            : data(data), N_X(N_X), N_Y(N_Y), N_Z(N_Z)
    {
        result = vector<complex<float_type>>(data.size());
        int rank = 1;
        int n[] = {(int)N_Z};
        int *inembed = n;
        int *onembed = n;
        int istride = 1;
        int ostride = istride;
        int idist = N_Z;
        int odist = idist;
        int howmany = N_X*N_Y;
        fft_plan = fftw_plan_many_dft(
                rank,
                n,
                howmany,
                reinterpret_cast<fftw_complex *>(data.data()),
                inembed,
                istride,
                idist,
                reinterpret_cast<fftw_complex *>(result.data()),
                onembed,
                ostride,
                odist,
                FFTW_FORWARD,
                FFTW_MEASURE);
    }

    vector<complex<float_type>>& run(){
        fftw_execute(fft_plan);
        return result;
    }

};


struct xz_fft_many {
    vector<complex<float_type>>& data;
    vector<complex<float_type>> result;
    fftw_plan fft_plan;
    index_type N_X;
    index_type N_Y;
    index_type N_Z;


    xz_fft_many(vector<complex<float_type>>& data,index_type N_X,index_type N_Y,index_type N_Z)
            : data(data), N_X(N_X), N_Y(N_Y), N_Z(N_Z)
    {
        result = vector<complex<float_type>>(data.size());
        int rank = 2;
        int n[] = {(int) N_X, (int) N_Z};
        int inembed[] = {(int) N_X, (int) (N_Z*N_Y)};
        int *onembed = inembed;
        int istride = 1;
        int ostride = 1;
        int idist = N_Z;
        int odist = N_Z;
        int howmany = N_Y;
        fft_plan = fftw_plan_many_dft(
                rank,
                n,
                howmany,
                reinterpret_cast<fftw_complex *>(data.data()),
                inembed,
                istride,
                idist,
                reinterpret_cast<fftw_complex *>(result.data()),
                onembed,
                ostride,
                odist,
                FFTW_FORWARD,FFTW_MEASURE);
    }

    vector<complex<float_type>>& run(){
        fftw_execute(fft_plan);
        return result;
    }

};

struct xz_fft_composition {
    vector<complex<float_type>>& data;
    index_type N_X;
    index_type N_Y;
    index_type N_Z;
    x_fft* xFft;
    z_fft* zFft;


    xz_fft_composition(vector<complex<float_type>>& data,index_type N_X,index_type N_Y,index_type N_Z)
            : data(data), N_X(N_X), N_Y(N_Y), N_Z(N_Z)
    {
        xFft = new x_fft(data,N_X,N_Y,N_Z);
        zFft = new z_fft(xFft->result,N_X,N_Y,N_Z);
    }

    vector<complex<float_type>>& run(){
        xFft->run();
        return zFft->run();
    }

};

struct TestData{
    index_type N_X = 512;
    index_type N_Y = 16;
    index_type N_Z = 16;

    index_type ARRAY_SIZE = N_X * N_Y * N_Z;
    std::vector<complex<float_type>> data = get_data(ARRAY_SIZE);

    TestData() {
//        print(data,N_X,N_Y,N_Z);
    }
};

TestData testData;

struct SanityTest{
    SanityTest() {
        xz_fft_many fft_many(testData.data, testData.N_X, testData.N_Y, testData.N_Z);
        xz_fft_composition fft_composition(testData.data, testData.N_X, testData.N_Y, testData.N_Z);
        std::vector<complex<float_type>> fft_many_result =  fft_many.run();
        std::vector<complex<float_type>> fft_composition_result =  fft_composition.run();

        bool equal = std::equal(fft_composition_result.begin(), fft_composition_result.end(), fft_many_result.begin());
        assert(equal);
        if(equal){
            cout << "ok" << endl;
        }
    }
};

SanityTest sanityTest;

static void XZ_test_many(benchmark::State& state) {
    xz_fft_many fft(testData.data, testData.N_X, testData.N_Y, testData.N_Z);
    for (auto _ : state) {
        auto result = fft.run();
    }
}

static void XZ_test_composition(benchmark::State& state) {
    xz_fft_composition fft(testData.data, testData.N_X, testData.N_Y, testData.N_Z);
    for (auto _ : state) {
        auto result = fft.run();
    }
}

BENCHMARK(XZ_test_many)->Iterations(1000);
BENCHMARK(XZ_test_composition)->Iterations(1000);

BENCHMARK_MAIN();

如果我正确地进行了基准测试，对于不同的 N_X、N_Y、N_Z 组合，fftw_plan_many_dft 和组合方法之间存在一些显着差异。例如使用

    index_type N_X = 512;
    index_type N_Y = 16;
    index_type N_Z = 16;

我有近两倍的差异，有利于 fftw_plan_many_dft，但对于其他输入参数集，我经常发现组合方法更快，但没有那么多.

------------------------------------------------------------------------------
Benchmark                                    Time             CPU   Iterations
------------------------------------------------------------------------------
XZ_test_many/iterations:1000           1412647 ns      1364813 ns         1000
XZ_test_composition/iterations:1000    2619807 ns      2542472 ns         1000

关于fft - 我可以使用 FFTW 沿 3D 数组的第一个和最后一个轴进行变换吗？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/73145984/

上一篇：react-native - Tailwind rn 已安装，但未设置样式 View

下一篇：javascript - 如何在 JavaScript 中 chop UTF8 字符串而不破坏多字节字符？

相关文章：

python - IFFT 在 Python 中返回什么？

linux - 为什么我不能用 LD_LIBRARY_PATH 覆盖动态库的搜索路径？

multithreading - 共享内存计算机上的多线程 FFTW 3.1.2

fftw - fftw 输出取决于输入的大小吗？

c++ - 基于 OpenCV DFT 的卷积被移位

iphone - 使用 Apple FFT 和 Accelerate 框架

iphone - 用于乐器调音器的 AurioTouch 和 FFT

C# 库做 fft 和 ifft？

julia - 如何在 Julia 中可视化信号的 FFT？

build - 如何在 Windows 上将 FFTW 与 cmake 一起使用？

©2024 IT工具网  联系我们