文档章节

离散傅里叶变换

Pulsar-V
 Pulsar-V
发布于 2017/07/24 16:43
字数 6710
阅读 56
收藏 2

离散傅里叶变换

    #include<iostream>  
    #include<math.h>  
    using namespace std;  
    #define PI 3.14159265354  
    struct complex{  
        double r,i;  
    };  
    complex multi(complex a,complex b){  
        complex tmp;  
        tmp.r=a.r*b.r-a.i*b.i;  
        tmp.i=a.r*b.i+a.i*b.r;  
        return tmp;  
    }  
    int fi(double in){  
        if((in-(int)in)>0.5) return (int)in+1;  
        else return (int)in;  
    }  
    /* 离散傅立叶正变换,输出[][2]数组实部在前,采样容量n可以任意正整数 */  
    void DFT(int *in,double **out,const int &n)  
    {  
        int i,j;  
        complex **W=new complex*[n];  
        for(i=0;i<n;i++){  
            W[i]=new complex[n];  
        }  
        complex *lis=new complex[(n-1)*(n-1)+1];  
        lis[0].r=1;lis[0].i=0;  
        lis[1].r=cos(2.0*PI/n);  
        lis[1].i=-1.0*sin(2.0*PI/n);  
        for(i=2;i<=(n-1)*(n-1);i++){  
            lis[i]=multi(lis[1],lis[i-1]);  
        }  
        for(i=0;i<n;i++){  
            for(j=0;j<n;j++){  
                W[i][j]=lis[i*j];  
            }  
        }  
        complex sum;  
        for(i=0;i<n;i++){  
            sum.r=0;sum.i=0;  
            for(j=0;j<n;j++){  
                sum.r+=in[j]*W[i][j].r;  
                sum.i+=in[j]*W[i][j].i;  
            }  
            out[i][0]=sum.r;  
            out[i][1]=sum.i;  
        }  
        for(i=0;i<n;i++) delete []W[i];  
        delete []W;  
        delete []lis;  
    }  
    /* 离散傅立叶逆变换 */  
    void IDFT(double **in,int *out,const int &n)  
    {  
        int i,j;  
        complex **W=new complex*[n];  
        for(i=0;i<n;i++){  
            W[i]=new complex[n];  
        }  
        complex *lis=new complex[(n-1)*(n-1)+1];  
        lis[0].r=1;lis[0].i=0;  
        lis[1].r=cos(2.0*PI/n);  
        lis[1].i=sin(2.0*PI/n);  
        for(i=2;i<=(n-1)*(n-1);i++){  
            lis[i]=multi(lis[1],lis[i-1]);  
        }  
        for(i=0;i<n;i++){  
            for(j=0;j<n;j++){  
                W[i][j]=lis[i*j];  
            }  
        }  
        complex sum;  
        for(i=0;i<n;i++){  
            sum.r=0;sum.i=0;  
            for(j=0;j<n;j++){  
                sum.r+=W[i][j].r*in[j][0]-W[i][j].i*in[j][1];  
                sum.i+=W[i][j].i*in[j][0]+W[i][j].r*in[j][1];  
            }  
            out[i]=fi(sum.r/n);  
        }  
        for(i=0;i<n;i++) delete []W[i];  
        delete []W;  
        delete []lis;  
    }  

OpenCV内部DFT变换实现(dft.hpp文件)

// fbc_cv is free software and uses the same licence as OpenCV  
// Email: fengbingchun@163.com  
#ifndef FBC_CV_DFT_HPP_  
#define FBC_CV_DFT_HPP_  
  
/* reference: include/opencv2/core.hpp 
              modules/core/src/dxt.cpp 
*/  
  
#include "core/mat.hpp"  
#include "core/core.hpp"  
  
namespace fbc {  
  
template<typename dump> static int DFTFactorize(int n, int* factors);  
template<typename dump> static void DFTInit(int n0, int nf, int* factors, int* itab, int elem_size, void* _wave, int inv_itab);  
template<typename T> static void DFT_32f(const Complex<T>* src, Complex<T>* dst, int n, int nf, const int* factors, const int* itab, const Complex<T>* wave, int tab_size, const void* spec, Complex<T>* buf, int flags, double _scale);  
template<typename T> static void RealDFT_32f(const T* src, T* dst, int n, int nf, int* factors, const int* itab, const Complex<T>* wave, int tab_size, const void* spec, Complex<T>* buf, int flags, double _scale);  
template<typename T> static void CCSIDFT_32f(const T* src, T* dst, int n, int nf, int* factors, const int* itab, const Complex<T>* wave, int tab_size, const void* spec, Complex<T>* buf, int flags, double _scale);  
template<typename _Tp, int chs> static void complementComplexOutput(Mat_<_Tp, chs>& dst, int len, int dft_dims);  
template<typename dump> static void CopyColumn(const uchar* _src, size_t src_step, uchar* _dst, size_t dst_step, int len, size_t elem_size);  
template<typename dump> static void CopyFrom2Columns(const uchar* _src, size_t src_step, uchar* _dst0, uchar* _dst1, int len, size_t elem_size);  
template<typename dump> static void CopyTo2Columns(const uchar* _src0, const uchar* _src1, uchar* _dst, size_t dst_step, int len, size_t elem_size);  
template<typename dump> static void ExpandCCS(uchar* _ptr, int n, int elem_size);  
  
static unsigned char bitrevTab[] =  
{  
    0x00, 0x80, 0x40, 0xc0, 0x20, 0xa0, 0x60, 0xe0, 0x10, 0x90, 0x50, 0xd0, 0x30, 0xb0, 0x70, 0xf0,  
    0x08, 0x88, 0x48, 0xc8, 0x28, 0xa8, 0x68, 0xe8, 0x18, 0x98, 0x58, 0xd8, 0x38, 0xb8, 0x78, 0xf8,  
    0x04, 0x84, 0x44, 0xc4, 0x24, 0xa4, 0x64, 0xe4, 0x14, 0x94, 0x54, 0xd4, 0x34, 0xb4, 0x74, 0xf4,  
    0x0c, 0x8c, 0x4c, 0xcc, 0x2c, 0xac, 0x6c, 0xec, 0x1c, 0x9c, 0x5c, 0xdc, 0x3c, 0xbc, 0x7c, 0xfc,  
    0x02, 0x82, 0x42, 0xc2, 0x22, 0xa2, 0x62, 0xe2, 0x12, 0x92, 0x52, 0xd2, 0x32, 0xb2, 0x72, 0xf2,  
    0x0a, 0x8a, 0x4a, 0xca, 0x2a, 0xaa, 0x6a, 0xea, 0x1a, 0x9a, 0x5a, 0xda, 0x3a, 0xba, 0x7a, 0xfa,  
    0x06, 0x86, 0x46, 0xc6, 0x26, 0xa6, 0x66, 0xe6, 0x16, 0x96, 0x56, 0xd6, 0x36, 0xb6, 0x76, 0xf6,  
    0x0e, 0x8e, 0x4e, 0xce, 0x2e, 0xae, 0x6e, 0xee, 0x1e, 0x9e, 0x5e, 0xde, 0x3e, 0xbe, 0x7e, 0xfe,  
    0x01, 0x81, 0x41, 0xc1, 0x21, 0xa1, 0x61, 0xe1, 0x11, 0x91, 0x51, 0xd1, 0x31, 0xb1, 0x71, 0xf1,  
    0x09, 0x89, 0x49, 0xc9, 0x29, 0xa9, 0x69, 0xe9, 0x19, 0x99, 0x59, 0xd9, 0x39, 0xb9, 0x79, 0xf9,  
    0x05, 0x85, 0x45, 0xc5, 0x25, 0xa5, 0x65, 0xe5, 0x15, 0x95, 0x55, 0xd5, 0x35, 0xb5, 0x75, 0xf5,  
    0x0d, 0x8d, 0x4d, 0xcd, 0x2d, 0xad, 0x6d, 0xed, 0x1d, 0x9d, 0x5d, 0xdd, 0x3d, 0xbd, 0x7d, 0xfd,  
    0x03, 0x83, 0x43, 0xc3, 0x23, 0xa3, 0x63, 0xe3, 0x13, 0x93, 0x53, 0xd3, 0x33, 0xb3, 0x73, 0xf3,  
    0x0b, 0x8b, 0x4b, 0xcb, 0x2b, 0xab, 0x6b, 0xeb, 0x1b, 0x9b, 0x5b, 0xdb, 0x3b, 0xbb, 0x7b, 0xfb,  
    0x07, 0x87, 0x47, 0xc7, 0x27, 0xa7, 0x67, 0xe7, 0x17, 0x97, 0x57, 0xd7, 0x37, 0xb7, 0x77, 0xf7,  
    0x0f, 0x8f, 0x4f, 0xcf, 0x2f, 0xaf, 0x6f, 0xef, 0x1f, 0x9f, 0x5f, 0xdf, 0x3f, 0xbf, 0x7f, 0xff  
};  
  
static const double DFTTab[][2] =  
{  
    { 1.00000000000000000, 0.00000000000000000 },  
    { -1.00000000000000000, 0.00000000000000000 },  
    { 0.00000000000000000, 1.00000000000000000 },  
    { 0.70710678118654757, 0.70710678118654746 },  
    { 0.92387953251128674, 0.38268343236508978 },  
    { 0.98078528040323043, 0.19509032201612825 },  
    { 0.99518472667219693, 0.09801714032956060 },  
    { 0.99879545620517241, 0.04906767432741802 },  
    { 0.99969881869620425, 0.02454122852291229 },  
    { 0.99992470183914450, 0.01227153828571993 },  
    { 0.99998117528260111, 0.00613588464915448 },  
    { 0.99999529380957619, 0.00306795676296598 },  
    { 0.99999882345170188, 0.00153398018628477 },  
    { 0.99999970586288223, 0.00076699031874270 },  
    { 0.99999992646571789, 0.00038349518757140 },  
    { 0.99999998161642933, 0.00019174759731070 },  
    { 0.99999999540410733, 0.00009587379909598 },  
    { 0.99999999885102686, 0.00004793689960307 },  
    { 0.99999999971275666, 0.00002396844980842 },  
    { 0.99999999992818922, 0.00001198422490507 },  
    { 0.99999999998204725, 0.00000599211245264 },  
    { 0.99999999999551181, 0.00000299605622633 },  
    { 0.99999999999887801, 0.00000149802811317 },  
    { 0.99999999999971945, 0.00000074901405658 },  
    { 0.99999999999992983, 0.00000037450702829 },  
    { 0.99999999999998246, 0.00000018725351415 },  
    { 0.99999999999999567, 0.00000009362675707 },  
    { 0.99999999999999889, 0.00000004681337854 },  
    { 0.99999999999999978, 0.00000002340668927 },  
    { 0.99999999999999989, 0.00000001170334463 },  
    { 1.00000000000000000, 0.00000000585167232 },  
    { 1.00000000000000000, 0.00000000292583616 }  
};  
  
#define BitRev(i,shift) \  
   ((int)((((unsigned)bitrevTab[(i)&255] << 24)+ \  
           ((unsigned)bitrevTab[((i)>> 8)&255] << 16)+ \  
           ((unsigned)bitrevTab[((i)>>16)&255] <<  8)+ \  
           ((unsigned)bitrevTab[((i)>>24)])) >> (shift)))  
  
enum { DFT_NO_PERMUTE = 256, DFT_COMPLEX_INPUT_OR_OUTPUT = 512 };  
  
// Performs a forward or inverse Discrete Fourier transform of a 1D or 2D floating-point array  
/* 
The function performs one of the following : 
-   Forward the Fourier transform of a 1D vector of N elements : 
    \f[Y = F^{ (N) }  \cdot X, \f] 
    where \f$F^{ (N) }_{ jk } = \exp(-2\pi i j k / N)\f$ and \f$i = \sqrt{ -1 }\f$ 
-   Inverse the Fourier transform of a 1D vector of N elements : 
    \f[\begin{ array }{l} X'=  \left (F^{(N)} \right )^{-1}  \cdot Y =  \left (F^{(N)} \right )^*  \cdot y  \\ X = (1/N)  \cdot X, \end{array}\f] 
    where \f$F^ *= \left(\textrm{ Re }(F^{ (N) }) - \textrm{ Im }(F^{ (N) })\right) ^ T\f$ 
-   Forward the 2D Fourier transform of a M x N matrix : 
    \f[Y = F^{ (M) }  \cdot X  \cdot F^{ (N) }\f] 
-   Inverse the 2D Fourier transform of a M x N matrix : 
    \f[
\begin{ array }{l} X'=  \left (F^{(M)} \right )^*  \cdot Y  \cdot \left (F^{(N)} \right )^* \\ X =  \frac{1}{M \cdot N} \cdot X' \end{ array }
\f] 
*/  
// support type: float, multi-channels  
template<typename _Tp, int chs1, int chs2>  
int dft(const Mat_<_Tp, chs1>& src0, Mat_<_Tp, chs2>& dst, int flags = 0, int nonzero_rows = 0)  
{  
    FBC_Assert(typeid(float).name() == typeid(_Tp).name());  
    FBC_Assert(chs1 == 1 || chs1 == 2 || chs2 == 1 || chs2 == 2);  
  
    AutoBuffer<uchar> buf;  
    int prev_len = 0, stage = 0;  
    bool inv = (flags & DFT_INVERSE) != 0;  
    int nf = 0, real_transform = src0.channels == 1 || (inv && (flags & DFT_REAL_OUTPUT) != 0);  
    int elem_size = (int)src0.elemSize1(), complex_elem_size = elem_size * 2;  
    int factors[34];  
    bool inplace_transform = false;  
  
    if (!inv && src0.channels == 1 && (flags & DFT_COMPLEX_OUTPUT)) {  
        FBC_Assert(chs2 == 2);  
  
        if (dst.empty()) {  
            dst = Mat_<_Tp, chs2>(src0.rows, src0.cols);  
        } else {  
            FBC_Assert(src0.rows == dst.rows && src0.cols == dst.cols);  
        }  
    } else if (inv && src0.channels == 2 && (flags & DFT_REAL_OUTPUT)) {  
        FBC_Assert(chs2 == 1);  
  
        if (dst.empty()) {  
            dst = Mat_<_Tp, chs2>(src0.rows, src0.cols);  
        } else {  
            FBC_Assert(src0.rows == dst.rows && src0.cols == dst.cols);  
        }  
    } else {  
        FBC_Assert(chs2 == chs1);  
  
        if (dst.empty()) {  
            dst = Mat_<_Tp, chs2>(src0.rows, src0.cols);  
        } else {  
            FBC_Assert(src0.rows == dst.rows && src0.cols == dst.cols);  
        }  
    }  
  
    if (!real_transform)  
        elem_size = complex_elem_size;  
  
    if (src0.cols == 1 && nonzero_rows > 0) {  
        FBC_Error("This mode (using nonzero_rows with a single-column matrix) breaks the function's logic, so it is prohibited.\n"  
            "For fast convolution/correlation use 2-column matrix or single-row matrix instead");  
    }  
  
    // determine, which transform to do first - row-wise  
    // (stage 0) or column-wise (stage 1) transform  
    if (!(flags & DFT_ROWS) && src0.rows > 1 &&  
        ((src0.cols == 1 && (!src0.isContinuous() || !dst.isContinuous())) ||  
        (src0.cols > 1 && inv && real_transform)))  
        stage = 1;  
  
    Mat_<_Tp, chs1> src = src0;  
    for (;;) {  
        double scale = 1;  
        uchar* wave = 0;  
        int* itab = 0;  
        uchar* ptr;  
        int i, len, count, sz = 0;  
        int use_buf = 0, odd_real = 0;  
  
        if (stage == 0) { // row-wise transform  
            len = !inv ? src.cols : dst.cols;  
            count = src.rows;  
            if (len == 1 && !(flags & DFT_ROWS)) {  
                len = !inv ? src.rows : dst.rows;  
                count = 1;  
            }  
            odd_real = real_transform && (len & 1);  
        } else {  
            len = dst.rows;  
            count = !inv ? src0.cols : dst.cols;  
            sz = 2 * len*complex_elem_size;  
        }  
  
        void *spec = 0;  
  
        {  
            if (len != prev_len)  
                nf = DFTFactorize<dump>(len, factors);  
  
            inplace_transform = factors[0] == factors[nf - 1];  
            sz += len*(complex_elem_size + sizeof(int));  
            i = nf > 1 && (factors[0] & 1) == 0;  
            if ((factors[i] & 1) != 0 && factors[i] > 5)  
                sz += (factors[i] + 1)*complex_elem_size;  
  
            if ((stage == 0 && ((src.data == dst.data && !inplace_transform) || odd_real)) ||  
                (stage == 1 && !inplace_transform)) {  
                use_buf = 1;  
                sz += len*complex_elem_size;  
            }  
        }  
  
        ptr = (uchar*)buf;  
        buf.allocate(sz + 32);  
        if (ptr != (uchar*)buf)  
            prev_len = 0; // because we release the buffer,  
        // force recalculation of twiddle factors and permutation table  
        ptr = (uchar*)buf;  
        if (!spec) {  
            wave = ptr;  
            ptr += len*complex_elem_size;  
            itab = (int*)ptr;  
            ptr = (uchar*)fbcAlignPtr<dump>(ptr + len*sizeof(int), 16);  
  
            if (len != prev_len || (!inplace_transform && inv && real_transform))  
                DFTInit<dump>(len, nf, factors, itab, complex_elem_size, wave, stage == 0 && inv && real_transform);  
            // otherwise reuse the tables calculated on the previous stage  
        }  
  
        if (stage == 0) {  
            uchar* tmp_buf = 0;  
            int dptr_offset = 0;  
            int dst_full_len = len*elem_size;  
            int _flags = (int)inv + (src.channels != dst.channels ? DFT_COMPLEX_INPUT_OR_OUTPUT : 0);  
  
            if (use_buf) {  
                tmp_buf = ptr;  
                ptr += len*complex_elem_size;  
                if (odd_real && !inv && len > 1 && !(_flags & DFT_COMPLEX_INPUT_OR_OUTPUT))  
                    dptr_offset = elem_size;  
            }  
  
            if (!inv && (_flags & DFT_COMPLEX_INPUT_OR_OUTPUT))  
                dst_full_len += (len & 1) ? elem_size : complex_elem_size;  
  
            if (count > 1 && !(flags & DFT_ROWS) && (!inv || !real_transform))  
                stage = 1;  
            else if (flags & FBC_DXT_SCALE)  
                scale = 1. / (len * (flags & DFT_ROWS ? 1 : count));  
  
            if (nonzero_rows <= 0 || nonzero_rows > count)  
                nonzero_rows = count;  
  
            for (i = 0; i < nonzero_rows; i++) {  
                const uchar* sptr = src.ptr(i);  
                uchar* dptr0 = dst.ptr(i);  
                uchar* dptr = dptr0;  
  
                if (tmp_buf)  
                    dptr = tmp_buf;  
  
                if (!real_transform) {  
                    DFT_32f<_Tp>(static_cast<const Complex<_Tp>*>((void*)sptr), static_cast<Complex<_Tp>*>((void*)dptr), len, nf, factors, itab, static_cast<const Complex<_Tp>*>((void*)wave), len, spec, static_cast<Complex<float>*>((void*)ptr), _flags, scale);  
                } else if (!inv) {  
                    RealDFT_32f<_Tp>(static_cast<const _Tp*>((void*)sptr), static_cast<_Tp*>((void*)dptr), len, nf, factors, itab, static_cast<const Complex<_Tp>*>((void*)wave), len, spec, static_cast<Complex<_Tp>*>((void*)ptr), _flags, scale);  
                } else if (inv) {  
                    CCSIDFT_32f<_Tp>(static_cast<const _Tp*>((void*)sptr), static_cast<_Tp*>((void*)dptr), len, nf, factors, itab, static_cast<const Complex<_Tp>*>((void*)wave), len, spec, static_cast<Complex<_Tp>*>((void*)ptr), _flags, scale);  
                } else {  
                    FBC_Error("no support type");  
                }  
  
                if (dptr != dptr0)  
                    memcpy(dptr0, dptr + dptr_offset, dst_full_len);  
            }  
  
            for (; i < count; i++) {  
                uchar* dptr0 = dst.ptr(i);  
                memset(dptr0, 0, dst_full_len);  
            }  
  
            if (stage != 1) {  
                if (!inv && real_transform && dst.channels == 2)  
                    complementComplexOutput(dst, nonzero_rows, 1);  
                break;  
            }  
  
            src = dst;  
        } else {  
            int a = 0, b = count;  
            uchar *buf0, *buf1, *dbuf0, *dbuf1;  
            const uchar* sptr0 = src.ptr();  
            uchar* dptr0 = dst.ptr();  
            buf0 = ptr;  
            ptr += len*complex_elem_size;  
            buf1 = ptr;  
            ptr += len*complex_elem_size;  
            dbuf0 = buf0, dbuf1 = buf1;  
  
            if (use_buf) {  
                dbuf1 = ptr;  
                dbuf0 = buf1;  
                ptr += len*complex_elem_size;  
            }  
  
            if (real_transform && inv && src.cols > 1)  
                stage = 0;  
            else if (flags & FBC_DXT_SCALE)  
                scale = 1. / (len * count);  
  
            if (real_transform) {  
                int even;  
                a = 1;  
                even = (count & 1) == 0;  
                b = (count + 1) / 2;  
                if (!inv) {  
                    memset(buf0, 0, len*complex_elem_size);  
                    CopyColumn<dump>(sptr0, src.step, buf0, complex_elem_size, len, elem_size);  
                    sptr0 += dst.channels*elem_size;  
                    if (even) {  
                        memset(buf1, 0, len*complex_elem_size);  
                        CopyColumn<dump>(sptr0 + (count - 2)*elem_size, src.step, buf1, complex_elem_size, len, elem_size);  
                    }  
                } else if (src.channels == 1) {  
                    CopyColumn<dump>(sptr0, src.step, buf0, elem_size, len, elem_size);  
                    ExpandCCS<dump>(buf0, len, elem_size);  
                    if (even) {  
                        CopyColumn<dump>(sptr0 + (count - 1)*elem_size, src.step, buf1, elem_size, len, elem_size);  
                        ExpandCCS<dump>(buf1, len, elem_size);  
                    }  
                    sptr0 += elem_size;  
                } else {  
                    CopyColumn<dump>(sptr0, src.step, buf0, complex_elem_size, len, complex_elem_size);  
                    if (even) {  
                        CopyColumn<dump>(sptr0 + b*complex_elem_size, src.step, buf1, complex_elem_size, len, complex_elem_size);  
                    }  
                    sptr0 += complex_elem_size;  
                }  
  
                if (even)  
                    DFT_32f<_Tp>(static_cast<const Complex<_Tp>*>((void*)buf1), static_cast<Complex<_Tp>*>((void*)dbuf1), len, nf, factors, itab, static_cast<const Complex<_Tp>*>((void*)wave), len, spec, static_cast<Complex<_Tp>*>((void*)ptr), inv, scale);  
                DFT_32f<_Tp>(static_cast<const Complex<_Tp>*>((void*)buf0), static_cast<Complex<_Tp>*>((void*)dbuf0), len, nf, factors, itab, static_cast<const Complex<_Tp>*>((void*)wave), len, spec, static_cast<Complex<_Tp>*>((void*)ptr), inv, scale);  
  
                if (dst.channels == 1) {  
                    if (!inv) {  
                        // copy the half of output vector to the first/last column.  
                        // before doing that, defgragment the vector  
                        memcpy(dbuf0 + elem_size, dbuf0, elem_size);  
                        CopyColumn<dump>(dbuf0 + elem_size, elem_size, dptr0, dst.step, len, elem_size);  
                        if (even) {  
                            memcpy(dbuf1 + elem_size, dbuf1, elem_size);  
                            CopyColumn<dump>(dbuf1 + elem_size, elem_size, dptr0 + (count - 1)*elem_size, dst.step, len, elem_size);  
                        }  
                        dptr0 += elem_size;  
                    } else {  
                        // copy the real part of the complex vector to the first/last column  
                        CopyColumn<dump>(dbuf0, complex_elem_size, dptr0, dst.step, len, elem_size);  
                        if (even)  
                            CopyColumn<dump>(dbuf1, complex_elem_size, dptr0 + (count - 1)*elem_size, dst.step, len, elem_size);  
                        dptr0 += elem_size;  
                    }  
                } else {  
                    assert(!inv);  
                    CopyColumn<dump>(dbuf0, complex_elem_size, dptr0, dst.step, len, complex_elem_size);  
                    if (even)  
                        CopyColumn<dump>(dbuf1, complex_elem_size, dptr0 + b*complex_elem_size, dst.step, len, complex_elem_size);  
                    dptr0 += complex_elem_size;  
                }  
            }  
  
            for (i = a; i < b; i += 2) {  
                if (i + 1 < b) {  
                    CopyFrom2Columns<dump>(sptr0, src.step, buf0, buf1, len, complex_elem_size);  
                    DFT_32f<_Tp>(static_cast<const Complex<_Tp>*>((void*)buf1), static_cast<Complex<_Tp>*>((void*)dbuf1), len, nf, factors, itab, static_cast<const Complex<_Tp>*>((void*)wave), len, spec, static_cast<Complex<_Tp>*>((void*)ptr), inv, scale);  
                } else  
                    CopyColumn<dump>(sptr0, src.step, buf0, complex_elem_size, len, complex_elem_size);  
  
                DFT_32f<_Tp>(static_cast<const Complex<_Tp>*>((void*)buf0), static_cast<Complex<_Tp>*>((void*)dbuf0), len, nf, factors, itab, static_cast<const Complex<_Tp>*>((void*)wave), len, spec, static_cast<Complex<_Tp>*>((void*)ptr), inv, scale);  
  
                if (i + 1 < b)  
                    CopyTo2Columns<dump>(dbuf0, dbuf1, dptr0, dst.step, len, complex_elem_size);  
                else  
                    CopyColumn<dump>(dbuf0, complex_elem_size, dptr0, dst.step, len, complex_elem_size);  
                sptr0 += 2 * complex_elem_size;  
                dptr0 += 2 * complex_elem_size;  
            }  
  
            if (stage != 0) {  
                if (!inv && real_transform && dst.channels == 2 && len > 1)  
                    complementComplexOutput(dst, len, 2);  
                break;  
            }  
  
            src = dst;  
        }  
    }  
  
    return 0;  
}  
  
template<typename T> struct DFT_VecR4  
{  
    int operator()(Complex<T>*, int, int, int&, const Complex<T>*) const { return 1; }  
};  
  
// mixed-radix complex discrete Fourier transform: double-precision version  
template<typename T>  
static void DFT_32f(const Complex<T>* src, Complex<T>* dst, int n, int nf, const int* factors, const int* itab, const Complex<T>* wave, int tab_size, const void* spec, Complex<T>* buf, int flags, double _scale)  
{  
    static const T sin_120 = (T)0.86602540378443864676372317075294;  
    static const T fft5_2 = (T)0.559016994374947424102293417182819;  
    static const T fft5_3 = (T)-0.951056516295153572116439333379382;  
    static const T fft5_4 = (T)-1.538841768587626701285145288018455;  
    static const T fft5_5 = (T)0.363271264002680442947733378740309;  
  
    int n0 = n, f_idx, nx;  
    int inv = flags & DFT_INVERSE;  
    int dw0 = tab_size, dw;  
    int i, j, k;  
    Complex<T> t;  
    T scale = (T)_scale;  
    int tab_step;  
  
    tab_step = tab_size == n ? 1 : tab_size == n * 2 ? 2 : tab_size / n;  
  
    // 0. shuffle data  
    if (dst != src) {  
        assert((flags & DFT_NO_PERMUTE) == 0);  
        if (!inv) {  
            for (i = 0; i <= n - 2; i += 2, itab += 2 * tab_step) {  
                int k0 = itab[0], k1 = itab[tab_step];  
                assert((unsigned)k0 < (unsigned)n && (unsigned)k1 < (unsigned)n);  
                dst[i] = src[k0]; dst[i + 1] = src[k1];  
            }  
  
            if (i < n)  
                dst[n - 1] = src[n - 1];  
        } else {  
            for (i = 0; i <= n - 2; i += 2, itab += 2 * tab_step) {  
                int k0 = itab[0], k1 = itab[tab_step];  
                assert((unsigned)k0 < (unsigned)n && (unsigned)k1 < (unsigned)n);  
                t.re = src[k0].re; t.im = -src[k0].im;  
                dst[i] = t;  
                t.re = src[k1].re; t.im = -src[k1].im;  
                dst[i + 1] = t;  
            }  
  
            if (i < n) {  
                t.re = src[n - 1].re; t.im = -src[n - 1].im;  
                dst[i] = t;  
            }  
        }  
    } else {  
        if ((flags & DFT_NO_PERMUTE) == 0) {  
            FBC_Assert(factors[0] == factors[nf - 1]);  
            if (nf == 1) {  
                if ((n & 3) == 0) {  
                    int n2 = n / 2;  
                    Complex<T>* dsth = dst + n2;  
  
                    for (i = 0; i < n2; i += 2, itab += tab_step * 2) {  
                        j = itab[0];  
                        assert((unsigned)j < (unsigned)n2);  
  
                        FBC_SWAP(dst[i + 1], dsth[j], t);  
                        if (j > i) {  
                            FBC_SWAP(dst[i], dst[j], t);  
                            FBC_SWAP(dsth[i + 1], dsth[j + 1], t);  
                        }  
                    }  
                }  
                // else do nothing  
            } else {  
                for (i = 0; i < n; i++, itab += tab_step) {  
                    j = itab[0];  
                    assert((unsigned)j < (unsigned)n);  
                    if (j > i)  
                        FBC_SWAP(dst[i], dst[j], t);  
                }  
            }  
        }  
  
        if (inv) {  
            for (i = 0; i <= n - 2; i += 2) {  
                T t0 = -dst[i].im;  
                T t1 = -dst[i + 1].im;  
                dst[i].im = t0; dst[i + 1].im = t1;  
            }  
  
            if (i < n)  
                dst[n - 1].im = -dst[n - 1].im;  
        }  
    }  
  
    n = 1;  
    // 1. power-2 transforms  
    if ((factors[0] & 1) == 0) {  
        if (factors[0] >= 4) {  
            DFT_VecR4<T> vr4;  
            n = vr4(dst, factors[0], n0, dw0, wave);  
        }  
  
        // radix-4 transform  
        for (; n * 4 <= factors[0];) {  
            nx = n;  
            n *= 4;  
            dw0 /= 4;  
  
            for (i = 0; i < n0; i += n) {  
                Complex<T> *v0, *v1;  
                T r0, i0, r1, i1, r2, i2, r3, i3, r4, i4;  
  
                v0 = dst + i;  
                v1 = v0 + nx * 2;  
  
                r0 = v1[0].re; i0 = v1[0].im;  
                r4 = v1[nx].re; i4 = v1[nx].im;  
  
                r1 = r0 + r4; i1 = i0 + i4;  
                r3 = i0 - i4; i3 = r4 - r0;  
  
                r2 = v0[0].re; i2 = v0[0].im;  
                r4 = v0[nx].re; i4 = v0[nx].im;  
  
                r0 = r2 + r4; i0 = i2 + i4;  
                r2 -= r4; i2 -= i4;  
  
                v0[0].re = r0 + r1; v0[0].im = i0 + i1;  
                v1[0].re = r0 - r1; v1[0].im = i0 - i1;  
                v0[nx].re = r2 + r3; v0[nx].im = i2 + i3;  
                v1[nx].re = r2 - r3; v1[nx].im = i2 - i3;  
  
                for (j = 1, dw = dw0; j < nx; j++, dw += dw0) {  
                    v0 = dst + i + j;  
                    v1 = v0 + nx * 2;  
  
                    r2 = v0[nx].re*wave[dw * 2].re - v0[nx].im*wave[dw * 2].im;  
                    i2 = v0[nx].re*wave[dw * 2].im + v0[nx].im*wave[dw * 2].re;  
                    r0 = v1[0].re*wave[dw].im + v1[0].im*wave[dw].re;  
                    i0 = v1[0].re*wave[dw].re - v1[0].im*wave[dw].im;  
                    r3 = v1[nx].re*wave[dw * 3].im + v1[nx].im*wave[dw * 3].re;  
                    i3 = v1[nx].re*wave[dw * 3].re - v1[nx].im*wave[dw * 3].im;  
  
                    r1 = i0 + i3; i1 = r0 + r3;  
                    r3 = r0 - r3; i3 = i3 - i0;  
                    r4 = v0[0].re; i4 = v0[0].im;  
  
                    r0 = r4 + r2; i0 = i4 + i2;  
                    r2 = r4 - r2; i2 = i4 - i2;  
  
                    v0[0].re = r0 + r1; v0[0].im = i0 + i1;  
                    v1[0].re = r0 - r1; v1[0].im = i0 - i1;  
                    v0[nx].re = r2 + r3; v0[nx].im = i2 + i3;  
                    v1[nx].re = r2 - r3; v1[nx].im = i2 - i3;  
                }  
            }  
        }  
  
        for (; n < factors[0];) {  
            // do the remaining radix-2 transform  
            nx = n;  
            n *= 2;  
            dw0 /= 2;  
  
            for (i = 0; i < n0; i += n) {  
                Complex<T>* v = dst + i;  
                T r0 = v[0].re + v[nx].re;  
                T i0 = v[0].im + v[nx].im;  
                T r1 = v[0].re - v[nx].re;  
                T i1 = v[0].im - v[nx].im;  
                v[0].re = r0; v[0].im = i0;  
                v[nx].re = r1; v[nx].im = i1;  
  
                for (j = 1, dw = dw0; j < nx; j++, dw += dw0) {  
                    v = dst + i + j;  
                    r1 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im;  
                    i1 = v[nx].im*wave[dw].re + v[nx].re*wave[dw].im;  
                    r0 = v[0].re; i0 = v[0].im;  
  
                    v[0].re = r0 + r1; v[0].im = i0 + i1;  
                    v[nx].re = r0 - r1; v[nx].im = i0 - i1;  
                }  
            }  
        }  
    }  
  
    // 2. all the other transforms  
    for (f_idx = (factors[0] & 1) ? 0 : 1; f_idx < nf; f_idx++) {  
        int factor = factors[f_idx];  
        nx = n;  
        n *= factor;  
        dw0 /= factor;  
  
        if (factor == 3) {  
            // radix-3  
            for (i = 0; i < n0; i += n) {  
                Complex<T>* v = dst + i;  
  
                T r1 = v[nx].re + v[nx * 2].re;  
                T i1 = v[nx].im + v[nx * 2].im;  
                T r0 = v[0].re;  
                T i0 = v[0].im;  
                T r2 = sin_120*(v[nx].im - v[nx * 2].im);  
                T i2 = sin_120*(v[nx * 2].re - v[nx].re);  
                v[0].re = r0 + r1; v[0].im = i0 + i1;  
                r0 -= (T)0.5*r1; i0 -= (T)0.5*i1;  
                v[nx].re = r0 + r2; v[nx].im = i0 + i2;  
                v[nx * 2].re = r0 - r2; v[nx * 2].im = i0 - i2;  
  
                for (j = 1, dw = dw0; j < nx; j++, dw += dw0) {  
                    v = dst + i + j;  
                    r0 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im;  
                    i0 = v[nx].re*wave[dw].im + v[nx].im*wave[dw].re;  
                    i2 = v[nx * 2].re*wave[dw * 2].re - v[nx * 2].im*wave[dw * 2].im;  
                    r2 = v[nx * 2].re*wave[dw * 2].im + v[nx * 2].im*wave[dw * 2].re;  
                    r1 = r0 + i2; i1 = i0 + r2;  
  
                    r2 = sin_120*(i0 - r2); i2 = sin_120*(i2 - r0);  
                    r0 = v[0].re; i0 = v[0].im;  
                    v[0].re = r0 + r1; v[0].im = i0 + i1;  
                    r0 -= (T)0.5*r1; i0 -= (T)0.5*i1;  
                    v[nx].re = r0 + r2; v[nx].im = i0 + i2;  
                    v[nx * 2].re = r0 - r2; v[nx * 2].im = i0 - i2;  
                }  
            }  
        } else if (factor == 5) {  
            // radix-5  
            for (i = 0; i < n0; i += n) {  
                for (j = 0, dw = 0; j < nx; j++, dw += dw0) {  
                    Complex<T>* v0 = dst + i + j;  
                    Complex<T>* v1 = v0 + nx * 2;  
                    Complex<T>* v2 = v1 + nx * 2;  
  
                    T r0, i0, r1, i1, r2, i2, r3, i3, r4, i4, r5, i5;  
  
                    r3 = v0[nx].re*wave[dw].re - v0[nx].im*wave[dw].im;  
                    i3 = v0[nx].re*wave[dw].im + v0[nx].im*wave[dw].re;  
                    r2 = v2[0].re*wave[dw * 4].re - v2[0].im*wave[dw * 4].im;  
                    i2 = v2[0].re*wave[dw * 4].im + v2[0].im*wave[dw * 4].re;  
  
                    r1 = r3 + r2; i1 = i3 + i2;  
                    r3 -= r2; i3 -= i2;  
  
                    r4 = v1[nx].re*wave[dw * 3].re - v1[nx].im*wave[dw * 3].im;  
                    i4 = v1[nx].re*wave[dw * 3].im + v1[nx].im*wave[dw * 3].re;  
                    r0 = v1[0].re*wave[dw * 2].re - v1[0].im*wave[dw * 2].im;  
                    i0 = v1[0].re*wave[dw * 2].im + v1[0].im*wave[dw * 2].re;  
  
                    r2 = r4 + r0; i2 = i4 + i0;  
                    r4 -= r0; i4 -= i0;  
  
                    r0 = v0[0].re; i0 = v0[0].im;  
                    r5 = r1 + r2; i5 = i1 + i2;  
  
                    v0[0].re = r0 + r5; v0[0].im = i0 + i5;  
  
                    r0 -= (T)0.25*r5; i0 -= (T)0.25*i5;  
                    r1 = fft5_2*(r1 - r2); i1 = fft5_2*(i1 - i2);  
                    r2 = -fft5_3*(i3 + i4); i2 = fft5_3*(r3 + r4);  
  
                    i3 *= -fft5_5; r3 *= fft5_5;  
                    i4 *= -fft5_4; r4 *= fft5_4;  
  
                    r5 = r2 + i3; i5 = i2 + r3;  
                    r2 -= i4; i2 -= r4;  
  
                    r3 = r0 + r1; i3 = i0 + i1;  
                    r0 -= r1; i0 -= i1;  
  
                    v0[nx].re = r3 + r2; v0[nx].im = i3 + i2;  
                    v2[0].re = r3 - r2; v2[0].im = i3 - i2;  
  
                    v1[0].re = r0 + r5; v1[0].im = i0 + i5;  
                    v1[nx].re = r0 - r5; v1[nx].im = i0 - i5;  
                }  
            }  
        } else {  
            // radix-"factor" - an odd number  
            int p, q, factor2 = (factor - 1) / 2;  
            int d, dd, dw_f = tab_size / factor;  
            Complex<T>* a = buf;  
            Complex<T>* b = buf + factor2;  
  
            for (i = 0; i < n0; i += n) {  
                for (j = 0, dw = 0; j < nx; j++, dw += dw0) {  
                    Complex<T>* v = dst + i + j;  
                    Complex<T> v_0 = v[0];  
                    Complex<T> vn_0 = v_0;  
  
                    if (j == 0) {  
                        for (p = 1, k = nx; p <= factor2; p++, k += nx) {  
                            T r0 = v[k].re + v[n - k].re;  
                            T i0 = v[k].im - v[n - k].im;  
                            T r1 = v[k].re - v[n - k].re;  
                            T i1 = v[k].im + v[n - k].im;  
  
                            vn_0.re += r0; vn_0.im += i1;  
                            a[p - 1].re = r0; a[p - 1].im = i0;  
                            b[p - 1].re = r1; b[p - 1].im = i1;  
                        }  
                    } else {  
                        const Complex<T>* wave_ = wave + dw*factor;  
                        d = dw;  
  
                        for (p = 1, k = nx; p <= factor2; p++, k += nx, d += dw) {  
                            T r2 = v[k].re*wave[d].re - v[k].im*wave[d].im;  
                            T i2 = v[k].re*wave[d].im + v[k].im*wave[d].re;  
  
                            T r1 = v[n - k].re*wave_[-d].re - v[n - k].im*wave_[-d].im;  
                            T i1 = v[n - k].re*wave_[-d].im + v[n - k].im*wave_[-d].re;  
  
                            T r0 = r2 + r1;  
                            T i0 = i2 - i1;  
                            r1 = r2 - r1;  
                            i1 = i2 + i1;  
  
                            vn_0.re += r0; vn_0.im += i1;  
                            a[p - 1].re = r0; a[p - 1].im = i0;  
                            b[p - 1].re = r1; b[p - 1].im = i1;  
                        }  
                    }  
  
                    v[0] = vn_0;  
  
                    for (p = 1, k = nx; p <= factor2; p++, k += nx) {  
                        Complex<T> s0 = v_0, s1 = v_0;  
                        d = dd = dw_f*p;  
  
                        for (q = 0; q < factor2; q++) {  
                            T r0 = wave[d].re * a[q].re;  
                            T i0 = wave[d].im * a[q].im;  
                            T r1 = wave[d].re * b[q].im;  
                            T i1 = wave[d].im * b[q].re;  
  
                            s1.re += r0 + i0; s0.re += r0 - i0;  
                            s1.im += r1 - i1; s0.im += r1 + i1;  
  
                            d += dd;  
                            d -= -(d >= tab_size) & tab_size;  
                        }  
  
                        v[k] = s0;  
                        v[n - k] = s1;  
                    }  
                }  
            }  
        }  
    }  
  
    if (scale != 1) {  
        T re_scale = scale, im_scale = scale;  
        if (inv)  
            im_scale = -im_scale;  
  
        for (i = 0; i < n0; i++) {  
            T t0 = dst[i].re*re_scale;  
            T t1 = dst[i].im*im_scale;  
            dst[i].re = t0;  
            dst[i].im = t1;  
        }  
    } else if (inv) {  
        for (i = 0; i <= n0 - 2; i += 2) {  
            T t0 = -dst[i].im;  
            T t1 = -dst[i + 1].im;  
            dst[i].im = t0;  
            dst[i + 1].im = t1;  
        }  
  
        if (i < n0)  
            dst[n0 - 1].im = -dst[n0 - 1].im;  
    }  
}  
  
/* FFT of real vector 
   output vector format: 
    re(0), re(1), im(1), ... , re(n/2-1), im((n+1)/2-1) [, re((n+1)/2)] OR 
    re(0), 0, re(1), im(1), ..., re(n/2-1), im((n+1)/2-1) [, re((n+1)/2), 0] */  
template<typename T>  
static void RealDFT_32f(const T* src, T* dst, int n, int nf, int* factors, const int* itab, const Complex<T>* wave, int tab_size, const void* spec, Complex<T>* buf, int flags, double _scale)  
{  
    int complex_output = (flags & DFT_COMPLEX_INPUT_OR_OUTPUT) != 0;  
    T scale = (T)_scale;  
    int j, n2 = n >> 1;  
    dst += complex_output;  
  
    assert(tab_size == n);  
  
    if (n == 1) {  
        dst[0] = src[0] * scale;  
    } else if (n == 2) {  
        T t = (src[0] + src[1])*scale;  
        dst[1] = (src[0] - src[1])*scale;  
        dst[0] = t;  
    } else if (n & 1) {  
        dst -= complex_output;  
        Complex<T>* _dst = (Complex<T>*)dst;  
        _dst[0].re = src[0] * scale;  
        _dst[0].im = 0;  
        for (j = 1; j < n; j += 2) {  
            T t0 = src[itab[j]] * scale;  
            T t1 = src[itab[j + 1]] * scale;  
            _dst[j].re = t0;  
            _dst[j].im = 0;  
            _dst[j + 1].re = t1;  
            _dst[j + 1].im = 0;  
        }  
        DFT_32f(_dst, _dst, n, nf, factors, itab, wave, tab_size, 0, buf, DFT_NO_PERMUTE, 1);  
        if (!complex_output)  
            dst[1] = dst[0];  
    } else {  
        T t0, t;  
        T h1_re, h1_im, h2_re, h2_im;  
        T scale2 = scale*(T)0.5;  
        factors[0] >>= 1;  
  
        DFT_32f((Complex<T>*)src, (Complex<T>*)dst, n2, nf - (factors[0] == 1), factors + (factors[0] == 1), itab, wave, tab_size, 0, buf, 0, 1);  
        factors[0] <<= 1;  
  
        t = dst[0] - dst[1];  
        dst[0] = (dst[0] + dst[1])*scale;  
        dst[1] = t*scale;  
  
        t0 = dst[n2];  
        t = dst[n - 1];  
        dst[n - 1] = dst[1];  
  
        for (j = 2, wave++; j < n2; j += 2, wave++) {  
            /* calc odd */  
            h2_re = scale2*(dst[j + 1] + t);  
            h2_im = scale2*(dst[n - j] - dst[j]);  
  
            /* calc even */  
            h1_re = scale2*(dst[j] + dst[n - j]);  
            h1_im = scale2*(dst[j + 1] - t);  
  
            /* rotate */  
            t = h2_re*wave->re - h2_im*wave->im;  
            h2_im = h2_re*wave->im + h2_im*wave->re;  
            h2_re = t;  
            t = dst[n - j - 1];  
  
            dst[j - 1] = h1_re + h2_re;  
            dst[n - j - 1] = h1_re - h2_re;  
            dst[j] = h1_im + h2_im;  
            dst[n - j] = h2_im - h1_im;  
        }  
  
        if (j <= n2) {  
            dst[n2 - 1] = t0*scale;  
            dst[n2] = -t*scale;  
        }  
    }  
  
    if (complex_output && ((n & 1) == 0 || n == 1)) {  
        dst[-1] = dst[0];  
        dst[0] = 0;  
        if (n > 1)  
            dst[n] = 0;  
    }  
}  
  
/* Inverse FFT of complex conjugate-symmetric vector 
   input vector format: 
    re[0], re[1], im[1], ... , re[n/2-1], im[n/2-1], re[n/2] OR 
    re(0), 0, re(1), im(1), ..., re(n/2-1), im((n+1)/2-1) [, re((n+1)/2), 0] */  
template<typename T>  
static void CCSIDFT_32f(const T* src, T* dst, int n, int nf, int* factors, const int* itab, const Complex<T>* wave, int tab_size, const void* spec, Complex<T>* buf, int flags, double _scale)  
{  
    int complex_input = (flags & DFT_COMPLEX_INPUT_OR_OUTPUT) != 0;  
    int j, k, n2 = (n + 1) >> 1;  
    T scale = (T)_scale;  
    T save_s1 = 0.;  
    T t0, t1, t2, t3, t;  
  
    assert(tab_size == n);  
  
    if (complex_input) {  
        assert(src != dst);  
        save_s1 = src[1];  
        ((T*)src)[1] = src[0];  
        src++;  
    }  
  
    if (n == 1) {  
        dst[0] = (T)(src[0] * scale);  
    } else if (n == 2) {  
        t = (src[0] + src[1])*scale;  
        dst[1] = (src[0] - src[1])*scale;  
        dst[0] = t;  
    } else if (n & 1) {  
        Complex<T>* _src = (Complex<T>*)(src - 1);  
        Complex<T>* _dst = (Complex<T>*)dst;  
  
        _dst[0].re = src[0];  
        _dst[0].im = 0;  
        for (j = 1; j < n2; j++) {  
            int k0 = itab[j], k1 = itab[n - j];  
            t0 = _src[j].re; t1 = _src[j].im;  
            _dst[k0].re = t0; _dst[k0].im = -t1;  
            _dst[k1].re = t0; _dst[k1].im = t1;  
        }  
  
        DFT_32f(_dst, _dst, n, nf, factors, itab, wave, tab_size, 0, buf, DFT_NO_PERMUTE, 1.);  
        dst[0] *= scale;  
        for (j = 1; j < n; j += 2) {  
            t0 = dst[j * 2] * scale;  
            t1 = dst[j * 2 + 2] * scale;  
            dst[j] = t0;  
            dst[j + 1] = t1;  
        }  
    } else {  
        int inplace = src == dst;  
        const Complex<T>* w = wave;  
  
        t = src[1];  
        t0 = (src[0] + src[n - 1]);  
        t1 = (src[n - 1] - src[0]);  
        dst[0] = t0;  
        dst[1] = t1;  
  
        for (j = 2, w++; j < n2; j += 2, w++) {  
            T h1_re, h1_im, h2_re, h2_im;  
  
            h1_re = (t + src[n - j - 1]);  
            h1_im = (src[j] - src[n - j]);  
  
            h2_re = (t - src[n - j - 1]);  
            h2_im = (src[j] + src[n - j]);  
  
            t = h2_re*w->re + h2_im*w->im;  
            h2_im = h2_im*w->re - h2_re*w->im;  
            h2_re = t;  
  
            t = src[j + 1];  
            t0 = h1_re - h2_im;  
            t1 = -h1_im - h2_re;  
            t2 = h1_re + h2_im;  
            t3 = h1_im - h2_re;  
  
            if (inplace) {  
                dst[j] = t0;  
                dst[j + 1] = t1;  
                dst[n - j] = t2;  
                dst[n - j + 1] = t3;  
            } else {  
                int j2 = j >> 1;  
                k = itab[j2];  
                dst[k] = t0;  
                dst[k + 1] = t1;  
                k = itab[n2 - j2];  
                dst[k] = t2;  
                dst[k + 1] = t3;  
            }  
        }  
  
        if (j <= n2) {  
            t0 = t * 2;  
            t1 = src[n2] * 2;  
  
            if (inplace) {  
                dst[n2] = t0;  
                dst[n2 + 1] = t1;  
            } else {  
                k = itab[n2];  
                dst[k * 2] = t0;  
                dst[k * 2 + 1] = t1;  
            }  
        }  
  
        factors[0] >>= 1;  
        DFT_32f((Complex<T>*)dst, (Complex<T>*)dst, n2,  
            nf - (factors[0] == 1),  
            factors + (factors[0] == 1), itab,  
            wave, tab_size, 0, buf,  
            inplace ? 0 : DFT_NO_PERMUTE, 1.);  
        factors[0] <<= 1;  
  
        for (j = 0; j < n; j += 2) {  
            t0 = dst[j] * scale;  
            t1 = dst[j + 1] * (-scale);  
            dst[j] = t0;  
            dst[j + 1] = t1;  
        }  
    }  
    if (complex_input)  
        ((T*)src)[0] = (T)save_s1;  
}  
  
template<typename _Tp, int chs>  
static void complementComplexOutput(Mat_<_Tp, chs>& dst, int len, int dft_dims)  
{  
    int i, n = dst.cols;  
    size_t elem_size = dst.elemSize1();  
    if (elem_size == sizeof(float)) {  
        float* p0 = (float*)dst.ptr();  
        size_t dstep = dst.step / sizeof(p0[0]);  
        for (i = 0; i < len; i++) {  
            float* p = p0 + dstep*i;  
            float* q = dft_dims == 1 || i == 0 || i * 2 == len ? p : p0 + dstep*(len - i);  
  
            for (int j = 1; j < (n + 1) / 2; j++) {  
                p[(n - j) * 2] = q[j * 2];  
                p[(n - j) * 2 + 1] = -q[j * 2 + 1];  
            }  
        }  
    } else {  
        double* p0 = (double*)dst.ptr();  
        size_t dstep = dst.step / sizeof(p0[0]);  
        for (i = 0; i < len; i++) {  
            double* p = p0 + dstep*i;  
            double* q = dft_dims == 1 || i == 0 || i * 2 == len ? p : p0 + dstep*(len - i);  
  
            for (int j = 1; j < (n + 1) / 2; j++) {  
                p[(n - j) * 2] = q[j * 2];  
                p[(n - j) * 2 + 1] = -q[j * 2 + 1];  
            }  
        }  
    }  
}  
  
template<typename dump>  
static void CopyColumn(const uchar* _src, size_t src_step, uchar* _dst, size_t dst_step, int len, size_t elem_size)  
{  
    int i, t0, t1;  
    const int* src = (const int*)_src;  
    int* dst = (int*)_dst;  
    src_step /= sizeof(src[0]);  
    dst_step /= sizeof(dst[0]);  
  
    if (elem_size == sizeof(int)) {  
        for (i = 0; i < len; i++, src += src_step, dst += dst_step)  
            dst[0] = src[0];  
    } else if (elem_size == sizeof(int) * 2) {  
        for (i = 0; i < len; i++, src += src_step, dst += dst_step) {  
            t0 = src[0]; t1 = src[1];  
            dst[0] = t0; dst[1] = t1;  
        }  
    } else if (elem_size == sizeof(int) * 4) {  
        for (i = 0; i < len; i++, src += src_step, dst += dst_step) {  
            t0 = src[0]; t1 = src[1];  
            dst[0] = t0; dst[1] = t1;  
            t0 = src[2]; t1 = src[3];  
            dst[2] = t0; dst[3] = t1;  
        }  
    }  
}  
  
template<typename dump>  
static void CopyFrom2Columns(const uchar* _src, size_t src_step, uchar* _dst0, uchar* _dst1, int len, size_t elem_size)  
{  
    int i, t0, t1;  
    const int* src = (const int*)_src;  
    int* dst0 = (int*)_dst0;  
    int* dst1 = (int*)_dst1;  
    src_step /= sizeof(src[0]);  
  
    if (elem_size == sizeof(int)) {  
        for (i = 0; i < len; i++, src += src_step) {  
            t0 = src[0]; t1 = src[1];  
            dst0[i] = t0; dst1[i] = t1;  
        }  
    } else if (elem_size == sizeof(int) * 2) {  
        for (i = 0; i < len * 2; i += 2, src += src_step) {  
            t0 = src[0]; t1 = src[1];  
            dst0[i] = t0; dst0[i + 1] = t1;  
            t0 = src[2]; t1 = src[3];  
            dst1[i] = t0; dst1[i + 1] = t1;  
        }  
    } else if (elem_size == sizeof(int) * 4) {  
        for (i = 0; i < len * 4; i += 4, src += src_step) {  
            t0 = src[0]; t1 = src[1];  
            dst0[i] = t0; dst0[i + 1] = t1;  
            t0 = src[2]; t1 = src[3];  
            dst0[i + 2] = t0; dst0[i + 3] = t1;  
            t0 = src[4]; t1 = src[5];  
            dst1[i] = t0; dst1[i + 1] = t1;  
            t0 = src[6]; t1 = src[7];  
            dst1[i + 2] = t0; dst1[i + 3] = t1;  
        }  
    }  
}  
  
template<typename dump>  
static void CopyTo2Columns(const uchar* _src0, const uchar* _src1, uchar* _dst, size_t dst_step, int len, size_t elem_size)  
{  
    int i, t0, t1;  
    const int* src0 = (const int*)_src0;  
    const int* src1 = (const int*)_src1;  
    int* dst = (int*)_dst;  
    dst_step /= sizeof(dst[0]);  
  
    if (elem_size == sizeof(int)) {  
        for (i = 0; i < len; i++, dst += dst_step) {  
            t0 = src0[i]; t1 = src1[i];  
            dst[0] = t0; dst[1] = t1;  
        }  
    } else if (elem_size == sizeof(int) * 2) {  
        for (i = 0; i < len * 2; i += 2, dst += dst_step) {  
            t0 = src0[i]; t1 = src0[i + 1];  
            dst[0] = t0; dst[1] = t1;  
            t0 = src1[i]; t1 = src1[i + 1];  
            dst[2] = t0; dst[3] = t1;  
        }  
    } else if (elem_size == sizeof(int) * 4) {  
        for (i = 0; i < len * 4; i += 4, dst += dst_step) {  
            t0 = src0[i]; t1 = src0[i + 1];  
            dst[0] = t0; dst[1] = t1;  
            t0 = src0[i + 2]; t1 = src0[i + 3];  
            dst[2] = t0; dst[3] = t1;  
            t0 = src1[i]; t1 = src1[i + 1];  
            dst[4] = t0; dst[5] = t1;  
            t0 = src1[i + 2]; t1 = src1[i + 3];  
            dst[6] = t0; dst[7] = t1;  
        }  
    }  
}  
  
template<typename dump>  
static void ExpandCCS(uchar* _ptr, int n, int elem_size)  
{  
    int i;  
    if (elem_size == (int)sizeof(float)) {  
        float* p = (float*)_ptr;  
        for (i = 1; i < (n + 1) / 2; i++) {  
            p[(n - i) * 2] = p[i * 2 - 1];  
            p[(n - i) * 2 + 1] = -p[i * 2];  
        }  
        if ((n & 1) == 0) {  
            p[n] = p[n - 1];  
            p[n + 1] = 0.f;  
            n--;  
        }  
        for (i = n - 1; i > 0; i--)  
            p[i + 1] = p[i];  
        p[1] = 0.f;  
    } else {  
        double* p = (double*)_ptr;  
        for (i = 1; i < (n + 1) / 2; i++) {  
            p[(n - i) * 2] = p[i * 2 - 1];  
            p[(n - i) * 2 + 1] = -p[i * 2];  
        }  
        if ((n & 1) == 0) {  
            p[n] = p[n - 1];  
            p[n + 1] = 0.f;  
            n--;  
        }  
        for (i = n - 1; i > 0; i--)  
            p[i + 1] = p[i];  
        p[1] = 0.f;  
    }  
}  
  
template<typename dump>  
static int DFTFactorize(int n, int* factors)  
{  
    int nf = 0, f, i, j;  
  
    if (n <= 5) {  
        factors[0] = n;  
        return 1;  
    }  
  
    f = (((n - 1) ^ n) + 1) >> 1;  
    if (f > 1) {  
        factors[nf++] = f;  
        n = f == n ? 1 : n / f;  
    }  
  
    for (f = 3; n > 1;) {  
        int d = n / f;  
        if (d*f == n) {  
            factors[nf++] = f;  
            n = d;  
        } else {  
            f += 2;  
            if (f*f > n)  
                break;  
        }  
    }  
  
    if (n > 1)  
        factors[nf++] = n;  
  
    f = (factors[0] & 1) == 0;  
    for (i = f; i < (nf + f) / 2; i++)  
        std::swap(factors[i], factors[nf - i - 1 + f]);  
  
    return nf;  
}  
  
template<typename dump>  
static void DFTInit(int n0, int nf, int* factors, int* itab, int elem_size, void* _wave, int inv_itab)  
{  
    int digits[34], radix[34];  
    int n = factors[0], m = 0;  
    int* itab0 = itab;  
    int i, j, k;  
    Complex<double> w, w1;  
    double t;  
  
    if (n0 <= 5) {  
        itab[0] = 0;  
        itab[n0 - 1] = n0 - 1;  
  
        if (n0 != 4) {  
            for (i = 1; i < n0 - 1; i++)  
                itab[i] = i;  
        } else {  
            itab[1] = 2;  
            itab[2] = 1;  
        }  
        if (n0 == 5) {  
            if (elem_size == sizeof(Complex<double>))  
                ((Complex<double>*)_wave)[0] = Complex<double>(1., 0.);  
            else  
                ((Complex<float>*)_wave)[0] = Complex<float>(1.f, 0.f);  
        }  
        if (n0 != 4)  
            return;  
        m = 2;  
    } else {  
        // radix[] is initialized from index 'nf' down to zero  
        assert(nf < 34);  
        radix[nf] = 1;  
        digits[nf] = 0;  
        for (i = 0; i < nf; i++) {  
            digits[i] = 0;  
            radix[nf - i - 1] = radix[nf - i] * factors[nf - i - 1];  
        }  
  
        if (inv_itab && factors[0] != factors[nf - 1])  
            itab = (int*)_wave;  
  
        if ((n & 1) == 0) {  
            int a = radix[1], na2 = n*a >> 1, na4 = na2 >> 1;  
            for (m = 0; (unsigned)(1 << m) < (unsigned)n; m++)  
                ;  
            if (n <= 2) {  
                itab[0] = 0;  
                itab[1] = na2;  
            } else if (n <= 256) {  
                int shift = 10 - m;  
                for (i = 0; i <= n - 4; i += 4) {  
                    j = (bitrevTab[i >> 2] >> shift)*a;  
                    itab[i] = j;  
                    itab[i + 1] = j + na2;  
                    itab[i + 2] = j + na4;  
                    itab[i + 3] = j + na2 + na4;  
                }  
            } else {  
                int shift = 34 - m;  
                for (i = 0; i < n; i += 4) {  
                    int i4 = i >> 2;  
                    j = BitRev(i4, shift)*a;  
                    itab[i] = j;  
                    itab[i + 1] = j + na2;  
                    itab[i + 2] = j + na4;  
                    itab[i + 3] = j + na2 + na4;  
                }  
            }  
  
            digits[1]++;  
  
            if (nf >= 2) {  
                for (i = n, j = radix[2]; i < n0;) {  
                    for (k = 0; k < n; k++)  
                        itab[i + k] = itab[k] + j;  
                    if ((i += n) >= n0)  
                        break;  
                    j += radix[2];  
                    for (k = 1; ++digits[k] >= factors[k]; k++) {  
                        digits[k] = 0;  
                        j += radix[k + 2] - radix[k];  
                    }  
                }  
            }  
        } else {  
            for (i = 0, j = 0;;) {  
                itab[i] = j;  
                if (++i >= n0)  
                    break;  
                j += radix[1];  
                for (k = 0; ++digits[k] >= factors[k]; k++) {  
                    digits[k] = 0;  
                    j += radix[k + 2] - radix[k];  
                }  
            }  
        }  
  
        if (itab != itab0) {  
            itab0[0] = 0;  
            for (i = n0 & 1; i < n0; i += 2) {  
                int k0 = itab[i];  
                int k1 = itab[i + 1];  
                itab0[k0] = i;  
                itab0[k1] = i + 1;  
            }  
        }  
    }  
  
    if ((n0 & (n0 - 1)) == 0) {  
        w.re = w1.re = DFTTab[m][0];  
        w.im = w1.im = -DFTTab[m][1];  
    } else {  
        t = -FBC_PI * 2 / n0;  
        w.im = w1.im = sin(t);  
        w.re = w1.re = std::sqrt(1. - w1.im*w1.im);  
    }  
    n = (n0 + 1) / 2;  
  
    if (elem_size == sizeof(Complex<double>)) {  
        Complex<double>* wave = (Complex<double>*)_wave;  
  
        wave[0].re = 1.;  
        wave[0].im = 0.;  
  
        if ((n0 & 1) == 0) {  
            wave[n].re = -1.;  
            wave[n].im = 0;  
        }  
  
        for (i = 1; i < n; i++) {  
            wave[i] = w;  
            wave[n0 - i].re = w.re;  
            wave[n0 - i].im = -w.im;  
  
            t = w.re*w1.re - w.im*w1.im;  
            w.im = w.re*w1.im + w.im*w1.re;  
            w.re = t;  
        }  
    } else {  
        Complex<float>* wave = (Complex<float>*)_wave;  
        assert(elem_size == sizeof(Complex<float>));  
  
        wave[0].re = 1.f;  
        wave[0].im = 0.f;  
  
        if ((n0 & 1) == 0) {  
            wave[n].re = -1.f;  
            wave[n].im = 0.f;  
        }  
  
        for (i = 1; i < n; i++) {  
            wave[i].re = (float)w.re;  
            wave[i].im = (float)w.im;  
            wave[n0 - i].re = (float)w.re;  
            wave[n0 - i].im = (float)-w.im;  
  
            t = w.re*w1.re - w.im*w1.im;  
            w.im = w.re*w1.im + w.im*w1.re;  
            w.re = t;  
        }  
    }  
}  
  
// Calculates the inverse Discrete Fourier Transform of a 1D or 2D array  
template<typename _Tp, int chs1, int chs2>  
int idft(const Mat_<_Tp, chs1>& src, Mat_<_Tp, chs2>& dst, int flags = 0, int nonzero_rows = 0)  
{  
    dft(src, dst, flags | DFT_INVERSE, nonzero_rows);  
}  
  
} // namespace fbc  
  
#endif // FBC_CV_DFT_HPP_ 

© 著作权归作者所有

共有 人打赏支持
Pulsar-V

Pulsar-V

粉丝 49
博文 96
码字总数 76428
作品 1
成都
后端工程师
【工具使用系列】关于 MATLAB 图像变换,你需要知道的事

如何进行图像变换(频域处理) 正交变换 傅里叶变换 二维连续傅里叶变换 二维离散傅里叶变换(DFT) 快速傅里叶变换(FFT) 离散余弦变换 一维离散余弦变换 二维离散余弦变换 离散正弦变换 ...

AllenMoore
01/26
8
0
【信号处理系列】关于信号变换,你需要知道的事

如何进行信号变换 什么是信号变换 离散时间傅里叶分析 离散时间傅里叶变换(DTFT) DTFT的特性 LTI系统的频域表达式 模拟信号的采样和重构 z 变换 双向z变换 z变换的重要特性 z反变换 z域的系...

AllenMoore
02/03
0
0
关于对傅里叶变换的一些理解

近日以来,由于学习图像处理,感觉其对傅里叶变换等内容要求较高,故重整旗鼓又过了一遍信号系统等章节,做了不少实验,有所感悟,特记录下,以便备忘! 首先,对于傅里叶变换,最为需要理解...

文剑Boy
2015/04/06
0
0
【工具使用系列】关于 MATLAB 电路与系统分析,你需要知道的事

如何进行电路与系统分析 信号的可视化及自变量变换及运算 连续时间信号的可视化 离散时间信号的可视化 连续时间信号的自变量变换及运算 离散时间信号的自变量变换及运算 LTI系统的时域数值分...

AllenMoore
02/19
0
0
Matlab图像处理入门教程(菜鸟级)

1,图像的读入和显示: (1)图像读入:A=imread('文件名',文件格式) [X,map]=imread('文件名') %X代表索引图像矩阵,map为颜色映射表 (2)图像显示:image(A); 2. 图像写回 A=imwrite('文件...

枫言风语
2012/09/22
0
0

没有更多内容

加载失败,请刷新页面

加载更多

linux 系统的运行级别

运行级别 运行级别 | 含义 0 关机 1 单用户模式,可以想象为windows 的安全模式,主要用于修复系统 2 不完全的命令模式,不含NFS服务 3 完全的命令行模式,就是标准的字符界面 4 系统保留 5 ...

Linux学习笔记
今天
2
0
学习设计模式——命令模式

任何模式的出现,都是为了解决一些特定的场景的耦合问题,以达到对修改封闭,对扩展开放的效果。命令模式也不例外: 命令模式是为了解决命令的请求者和命令的实现者之间的耦合关系。 解决了这...

江左煤郎
今天
3
0
字典树收集(非线程安全,后续做线程安全改进)

将500W个单词放进一个数据结构进行存储,然后进行快速比对,判断一个单词是不是这个500W单词之中的;来了一个单词前缀,给出500w个单词中有多少个单词是该前缀. 1、这个需求首先需要设计好数据结...

算法之名
昨天
15
0
GRASP设计模式

此文参考了这篇博客,建议读者阅读原文。 面向对象(Object-Oriented,OO)是当下软件开发的主流方法。在OO分析与设计中,我们首先从问题领域中抽象出领域模型,在领域模型中以适当的粒度归纳...

克虏伯
昨天
1
0
Coding and Paper Letter(四十)

资源整理。 1 Coding: 1.Tomislav Hengl撰写的非官方作者指南:Michael Gould•Wouter Gerritsma。 UnofficialGuide4Authors 2.R语言包rwrfhydro,社区贡献的工具箱,用于管理,分析和可视化...

胖胖雕
昨天
1
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部