FourierTransform.h

Go to the documentation of this file.

/* Copyright 2016 Kristofer Björnson
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef COM_DAFER45_TBTK_FOURIER_TRANSFORM
#define COM_DAFER45_TBTK_FOURIER_TRANSFORM

#include "TBTK/CArray.h"
#include "TBTK/Index.h"

#include <fftw3.h>

#include <complex>
#include <vector>

namespace TBTK{

class FourierTransform{
public:
    template<typename DataType>
    class Plan{
    public:
        Plan(
            const CArray<DataType> &in,
            CArray<DataType> &out,
            const std::vector<unsigned int> &ranges,
            int sign
        );

        Plan(const Plan &plan) = delete;

        Plan(Plan &&plan);

        ~Plan();

        Plan& operator=(const Plan &plan) = delete;

        Plan& operator=(Plan &&plan);

        void setNormalizationFactor(double normalizationFactor);

        double getNormalizationFactor() const;
    private:
        fftw_plan *plan;

        double normalizationFactor;

        unsigned int size;

        const CArray<DataType> &input;

        CArray<DataType> &output;

        fftw_plan& getFFTWPlan();

        unsigned int getSize() const;

        CArray<DataType>& getInput();

        CArray<DataType>& getOutput();

        friend class FourierTransform;
    };

    template<typename DataType>
    class ForwardPlan : public Plan<DataType>{
    public:
        ForwardPlan(
            const CArray<DataType> &in,
            CArray<DataType> &out,
            const std::vector<unsigned int> &ranges
        ) : Plan<DataType>(in, out, ranges, -1){}
    };

    template<typename DataType>
    class InversePlan : public Plan<DataType>{
    public:
        InversePlan(
            const CArray<DataType> &in,
            CArray<DataType> &out,
            const std::vector<unsigned int> &ranges
        ) : Plan<DataType>(in, out, ranges, 1
        ){}
    };

    static void transform(
        const CArray<std::complex<double>> &in,
        CArray<std::complex<double>> &out,
        const std::vector<unsigned int> &ranges,
        int sign
    );

    template<typename DataType>
    static void transform(Plan<DataType> &plan);

    static void forward(
        const CArray<std::complex<double>> &in,
        CArray<std::complex<double>> &out,
        const std::vector<unsigned int> &ranges
    );

    static void inverse(
        const CArray<std::complex<double>> &in,
        CArray<std::complex<double>> &out,
        const std::vector<unsigned int> &ranges
    );

private:
};

template<typename DataType>
inline void FourierTransform::transform(Plan<DataType> &plan){
    fftw_execute(plan.getFFTWPlan());

    double normalizationFactor = plan.getNormalizationFactor();
    if(normalizationFactor != 1.){
        CArray<DataType> &output = plan.getOutput();
        for(unsigned int n = 0; n < plan.getSize(); n++)
            output[n] /= normalizationFactor;
    }
}

inline void FourierTransform::forward(
    const CArray<std::complex<double>> &in,
    CArray<std::complex<double>> &out,
    const std::vector<unsigned int> &ranges
){
    transform(in, out, ranges, -1);
}

inline void FourierTransform::inverse(
    const CArray<std::complex<double>> &in,
    CArray<std::complex<double>> &out,
    const std::vector<unsigned int> &ranges
){
    transform(in, out, ranges, 1);
}

template<typename DataType>
inline FourierTransform::Plan<DataType>::Plan(Plan &&plan){
    this->plan = plan.plan;
    plan.plan = nullptr;

    normalizationFactor = plan.normalizationFactor;
    size = plan.size;
    input = plan.input;
    output = plan.output;
}

template<typename DataType>
inline FourierTransform::Plan<DataType>::~Plan(){
    if(plan != nullptr){
        #pragma omp critical (TBTK_FOURIER_TRANSFORM)
        fftw_destroy_plan(*plan);

        delete plan;
    }

}

template<typename DataType>
inline FourierTransform::Plan<DataType>& FourierTransform::Plan<
    DataType
>::operator=(Plan &&rhs){
    if(this != &rhs){
        if(this->plan != nullptr){
            #pragma omp critical (TBTK_FOURIER_TRANSFORM)
            fftw_destroy_plan(*this->plan);

            delete this->plan;

            this->plan = rhs.plan;

            normalizationFactor = rhs.normalizationFactor;
            size = rhs.size;
            input = rhs.input;
            output = rhs.output;
        }
    }

    return *this;
}

template<typename DataType>
inline void FourierTransform::Plan<DataType>::setNormalizationFactor(
    double normalizationFactor
){
    this->normalizationFactor = normalizationFactor;
}

template<typename DataType>
inline double FourierTransform::Plan<DataType>::getNormalizationFactor() const{
    return normalizationFactor;
}

template<typename DataType>
inline fftw_plan& FourierTransform::Plan<DataType>::getFFTWPlan(){
    return *plan;
}

template<typename DataType>
inline unsigned int FourierTransform::Plan<DataType>::getSize() const{
    return size;
}

template<typename DataType>
inline CArray<DataType>& FourierTransform::Plan<DataType>::getInput(){
    return input;
}

template<typename DataType>
inline CArray<DataType>& FourierTransform::Plan<DataType>::getOutput(){
    return output;
}

};  //End of namespace TBTK

#endif