Search code examples
c++fftquaternionsfftw

Quaternion FFT using 1D FFT in C++


I am currently implemented research paper that use a Quaternion FFT in c++. However, I could not find any C++ library that support Quaternion FFT. After some surveys, I have found someone on the Internet said that it is possible to convert Quaternion FFT process to several 1D complex-to-complex FFT. Does anyone know how to do that?

I try to use FFTW++ library which supports some basic FFT methods to implement it. I would be really appreciate if anyone could help.


Solution

  • Thanks for Severin's help, I finally follow the instruction mentioned in this paper to seperate the quaternion FFT into two Complex-to-Complex 2D FFT and has successfully reproduce the result that shown in the paper.

    Something like this: ( Please tell me if I am wrong :) )

    #include <Array.h>
    #include <fftw++.h>
    
    using namespace std;
    using namespace utils;
    using namespace Array;
    using namespace fftwpp;
    
    void SaliencyMapHandler::quaternionFourierTransform(int dim1, int dim2, double* d1, double* d2, double* d3, double* d4) {
        // dim1 is the 1-st dimension of data, dim2 is the 2-nd dimension of data
        fftw::maxthreads = get_max_threads();
        size_t align = sizeof(Complex);
    
        array2<Complex> f1(dim1, dim2, align);
        array2<Complex> f2(dim1, dim2, align);
    
        fft2d forward_1(-1, f1);
        fft2d backward_1(1, f1);
        fft2d forward_2(-1, f2);
        fft2d backward_2(1, f2);
    
        for(int j=0; j<dim1; j++) {
            for(int i=0; i<dim2; i++) {
                f1(i,j) = Complex(d1[j*dim2 + i], d2[j*dim2 + i]);
                f2(i,j) = Complex(d3[j*dim2 + i], d4[j*dim2 + i]);
            }
        }
    
        forward_1.fft(f1);
        forward_2.fft(f2);
    
        // Do something on frequency domain
    
        backward_1.fftNormalized(f1);
        backward_2.fftNormalized(f2);
    
        for(int j=0; j<dim1; j++) {
            for(int i=0; i<dim2; i++) {
                double p1 = real(f1(i,j));
                double p2 = imag(f1(i,j));
                double p3 = real(f2(i,j));
                double p4 = imag(f2(i,j));
    
                // Do something after inverse transform
            }
        }
    }