二维傅里叶变换转一维傅里叶变换

Question

看了一个论文说二维傅里叶变换图谱经圆积分计算后可以得到一维傅里叶变换，求大神推导公式或者过程

Answer 1

 // ==============================================================================
// 快速离散傅里叶变换和功率谱
// 一维快速傅里叶变换FFT1和二维快速傅里叶变换FFT2
// 测试环境 C++ builder 2010
// MinGW (GCC) 4.5
// wuxuping 2010-08-10
// 一定要将文件保存为UTF-8的格式,这是MinGW的默认格式.
// ==============================================================================

#ifndef FFT_H
#define FFT_H
// --------------------
#include <vector>
#include<complex>
#include<bitset>
#include <iterator>
#include <algorithm>
#include <iostream>
#include <sstream>
using namespace std;

const double pai = 3.1415926535897932384626433832795028841971; // 圆周率

typedef vector<vector<complex<double> > >VVCD;
typedef vector<vector<double> >VVD;

// -------------------------------------------------
class FFT1 {
private:
    unsigned N; // 数据长度

    unsigned NFFT; // 使用快速傅里叶变化所计算的长度

    unsigned r; // 数据长度NFFT=2^r,且2^(r-1)<N<<NFFT

    vector<complex<double> >TS; // 时间域序列

    vector<complex<double> >FS; // 频率域序列FrequencySeries

    vector<double>PSD; // PSD功率谱密度(离散谱)

    vector<double>PS; // 返回PS=PSD*Frequency
    int sign;

    // sign > 0为正变换否则为逆变换
    // ------------------------------------
    void Set_NFFT_r(unsigned n);
    unsigned RevBit(unsigned num);
    void GetIndex(unsigned L, unsigned ALindex, unsigned &AL_1index0,
        unsigned &AL_1index1, unsigned &Windex);

    // ------------------------------------

    void doFFT1();

    // 执行傅里叶变换
    // ------------------------------------
public:

    // ------实序列构造函数-------------
    FFT1(vector<double>TimeSeries, int sign_ = 1) {
        N = TimeSeries.size();
        if (sign_ > 0) {
            sign = 1.0; // 正变换
        }
        else {
            sign = -1.0; // 逆变换
        }
        Set_NFFT_r(N);

        complex<double>W0(0.0, 0.0); // 旋转因子
        // --------------------------------
        for (unsigned i = 0; i < NFFT; i++) { // 赋初值
            if (i < N) {
                FS.push_back(W0);
                complex<double>Wi(TimeSeries[i], 0.0); // 旋转因子
                TS.push_back(Wi);
            }
            else {
                FS.push_back(W0); // 补零
                TS.push_back(W0); // 补零
            }

        }
        // --------------------------------

        doFFT1(); // 执行傅里叶变换
        // --------------------------------
    };

    // -------复序列构造函数-------------
    FFT1(vector<complex<double> >TimeSeries, int sign_ = 1) { // 赋初值
        N = TimeSeries.size();
        if (sign_ > 0) {
            sign = 1.0; // 正变换
        }
        else {
            sign = -1.0; // 逆变换
        }
        Set_NFFT_r(N);

        complex<double>W0(0.0, 0.0); // 旋转因子
        // --------------------------------
        for (unsigned i = 0; i < NFFT; i++) { // 赋初值
            if (i < N) {
                FS.push_back(W0);
                TS.push_back(TimeSeries[i]);
            }
            else {
                FS.push_back(W0); // 补零
                TS.push_back(W0); // 补零
            }

        }
        // --------------------------------

        doFFT1(); // 执行傅里叶变换
        // --------------------------------
    };

    // ----------成员函数------

    // ----------成员函数 -------------
    vector<complex<double> >GetFS() { // 返回频率域序列FrequencySeries
        return FS;
    };

    // ----------成员函数 -------------
    vector<double>GetFrequency() { // 返回频率
        vector<double>Frequency;
        // ----------------------
        for (unsigned int i = 0; i < FS.size(); i++) {
            Frequency.push_back(i);
        }
        return Frequency;
    };

    // ----------成员函数 -------------
    vector<double>GetPSD() { // 返回FS^2
        if (PSD.size() > 0) {
            PSD.clear();
        }
        // ----------------------
        for (unsigned int i = 0; i < FS.size(); i++) {
            PSD.push_back(real(FS[i] * conj(FS[i])));
        }
        return PSD;
    };

    vector<double>GetPS() { // 返回PS=PSD*Frequency
        if (PS.size() > 0) {
            PS.clear();
        }
        // ----------------------
        for (unsigned int i = 0; i < FS.size(); i++) {
            PS.push_back(i * 1.0 * real(FS[i] * conj(FS[i])));
        }
        return PS;
    };
    // ----------析构函数 -------------

    ~FFT1() {
    };
    // -----------------------------------
};

// ------设置NFFT和r的值--------
void FFT1::Set_NFFT_r(unsigned n) {
    bitset<32>bsn(n);
    r = 0;
    while ((n >> r) > 0) {
        r++;
    }

    if ((bsn.count() == 1) && (r > 0)) {
        r--;
    }
    bitset<32>bsNFFT(0);
    bsNFFT.set(r);
    NFFT = (unsigned)bsNFFT.to_ulong();
};

// ---- 比特倒序函数
unsigned FFT1::RevBit(unsigned num) {
    bitset<32>bs(num);
    bitset<32>bsR;
    // ---------逆序
    for (unsigned i = 0; i < r; i++) {
        bsR.set(i, bs[r - i - 1]);
    }
    return(unsigned)bsR.to_ulong();
};

// ------------------------------------------------------------------------------
// 计算公式中的索引AL(ALindex)=AL-1(AL_1index0)+AL-1(AL_1index1)W(Windex)
// ------------------------------------------------------------------------------
void FFT1::GetIndex(unsigned L, unsigned ALindex, unsigned &AL_1index0,
    unsigned &AL_1index1, unsigned &Windex) {
    // ALindex为AL()的索引
    // AL_1index0 = 0; // AL-1()的第一项索引
    // AL_1index1 = 0; // AL-1()的第二项索引
    // Windex = 0; // 相位角中的索引
    bitset<32>bs(ALindex);
    bitset<32>bs1(ALindex);
    bs1.set(r - L, 0);
    AL_1index0 = bs1.to_ulong();
    // -----------------------------------
    bitset<32>bs2(ALindex);
    bs2.set(r - L, 1);
    AL_1index1 = bs2.to_ulong();
    // -----------------------------------
    bitset<32>bs3; // 左L位
    for (unsigned i = 0; i < L; i++) {
        bs3.set(r - L + i, bs[r - i - 1]);
    }
    Windex = bs3.to_ulong();
}

// ------------------------------------
void FFT1::doFFT1() { // 一维快速Fourier变换
    unsigned AL_1index0 = 0; // AL-1()的第一项索引
    unsigned AL_1index1 = 0; // AL-1()的第二项索引
    unsigned Windex = 0; // 相位角中的索引
    double alpha; // 相位角
    // ----------------------------
    vector<complex<double> >X = TS; // X中间临时变量

    for (unsigned L = 1; L <= r; L++) { // 蝶形计算过程
        for (unsigned ALindex = 0; ALindex < NFFT; ALindex++) {
            // ALindex为AL()的索引
            GetIndex(L, ALindex, AL_1index0, AL_1index1, Windex);
            alpha = -2.0 * sign * pai * Windex / (1.0 * NFFT); // 旋转因子的相位角
            complex<double>W(cos(alpha), sin(alpha)); // 旋转因子
            FS[ALindex] = X[AL_1index0] + X[AL_1index1] * W;
        }
        X = FS; // 复制数组
    }
    //
    if (sign > 0) { // 为正变换
        for (unsigned i = 0; i < NFFT; i++) {
            FS[i] = X[RevBit(i)]; // 索引按二进制倒序
        }
    }
    else { // 为逆变换
        complex<double>WN(NFFT, 0.0); // 数据的个数
        for (unsigned i = 0; i < NFFT; i++) {
            FS[i] = X[RevBit(i)] / WN; // 索引按二进制倒序
        }
    }
}
// ==============================================================================
// 计算二维快速傅里叶变换FFT2
// ==============================================================================

// ----------------------------------------------
class FFT2 {

private:

    // ----------------------------------------------
    unsigned N1; // 数据行长度 0<=k1<N1  0<=j1<N1

    unsigned N2; // 数据列长度 0<=k2<N2  0<=j2<N2

    int sign; // sign > 0为正变换否则为逆变换

    VVCD Tk1k2; // 二维时间域序列,行k1列k2

    VVCD Yj1k2; // 中间序列对 Ak1k2每一行(k1)做傅里叶变换的结果

    VVCD Fj1j2; // 频率域序列FrequencySeries

    VVD PSD; // PSD功率谱密度(离散谱)

    VVD PS; // 返回PS=PSD*Frequency

    // ------------------------------------
    VVCD TranspositionVVCD(VVCD dv); // 矩阵转置
    // ------------------------------------

public:

    // ------实序列构造函数-------------
    FFT2(VVD TK1K2, int sign_ = 1) { // 赋初值

        N1 = TK1K2.size(); // 获取行数
        Tk1k2.resize(N1);
        for (unsigned k1 = 0; k1 < N1; k1++) {
            N2 = TK1K2[k1].size();
            Tk1k2[k1].resize(N2);
            for (unsigned k2 = 0; k2 < N2; k2++) {
                complex<double>W(TK1K2[k1][k2], 0.0);
                Tk1k2[k1][k2] = W;
            }

        }
        // -----------------------------------------
        for (unsigned k1 = 0; k1 < N1; k1++) {
            FFT1 fft1(Tk1k2[k1], sign_);
            Yj1k2.push_back(fft1.GetFS());
        }
        // -----------------------------------------
        Yj1k2 = TranspositionVVCD(Yj1k2); // 将矩阵转置
        // -----------------------------------------
        N2 = Yj1k2.size(); // 获取列数
        for (unsigned k2 = 0; k2 < N2; k2++) {
            FFT1 fft1(Yj1k2[k2], sign_);
            Fj1j2.push_back(fft1.GetFS());
        }
        // --------------
        Fj1j2 = TranspositionVVCD(Fj1j2); // 将矩阵转置
        // -----------------------
    };

    // -------复序列构造函数-------------
    FFT2(VVCD TK1K2, int sign_ = 1) { // 赋初值
        Tk1k2 = TK1K2;
        N1 = Tk1k2.size(); // 获取行数
        for (unsigned k1 = 0; k1 < N1; k1++) {
            FFT1 fft1(Tk1k2[k1], sign_);
            Yj1k2.push_back(fft1.GetFS());
        }
        // -----------------------------------------
        Yj1k2 = TranspositionVVCD(Yj1k2); // 将矩阵转置
        // -----------------------------------------
        N2 = Yj1k2.size(); // 获取列数
        for (unsigned k2 = 0; k2 < N2; k2++) {
            FFT1 fft1(Yj1k2[k2], sign_);
            Fj1j2.push_back(fft1.GetFS());
        }
        // --------------
        Fj1j2 = TranspositionVVCD(Fj1j2); // 将矩阵转置
        // -----------------------
    };

    // ----------成员函数 -------------
    VVCD GetFS() { // 返回频率域序列Fj1j2
        return Fj1j2;
    };

    // ----------成员函数 -------------
    VVD GetPSD() { // 返回FS^2
        PSD.resize(Fj1j2.size());
        for (unsigned i = 0; i < Fj1j2.size(); i++) {
            PSD[i].resize(Fj1j2[i].size());
        }
        // ----------------------
        for (unsigned i = 0; i < Fj1j2.size(); i++) {
            for (unsigned j = 0; j < Fj1j2[i].size(); j++) {
                PSD[i][j] = real(Fj1j2[i][j] * conj(Fj1j2[i][j]));
            }
        }
        return PSD;
    };

    // ----------成员函数 -------------
    VVD GetPS() { // 返回PS=PSD*Frequency
        PS.resize(Fj1j2.size());
        for (unsigned i = 0; i < Fj1j2.size(); i++) {
            PS[i].resize(Fj1j2[i].size());
        }
        // ----------------------
        for (unsigned i = 0; i < Fj1j2.size(); i++) {
            for (unsigned j = 0; j < Fj1j2[i].size(); j++) {
                PS[i][j] = (double)i * (double)j * real
                    (Fj1j2[i][j] * conj(Fj1j2[i][j]));
            }
        }
        return PS;
    };
    // ----------析构函数 -------------

    ~FFT2() {
    };
    // -----------------------------------
};

// ----------------------------------------
VVCD FFT2::TranspositionVVCD(VVCD dv) { // 将矩阵转置
    unsigned dvRow = dv.size();
    unsigned maxcol = 0;
    unsigned mincol = 99999;
    VVCD temp;
    if (dv.size() > 0) {
        // ------------------------------
        for (unsigned i = 0; i < dvRow; i++) {
            if (maxcol < dv[i].size()) {
                maxcol = dv[i].size();
            }
            if (mincol > dv[i].size()) {
                mincol = dv[i].size();
            }
        }
        // ------------------------------
        if (mincol == maxcol && mincol > 0) {
            temp.resize(mincol);
            for (unsigned i = 0; i < mincol; i++) {
                temp[i].resize(dvRow);
            }
            for (unsigned i = 0; i < dvRow; i++) {
                for (unsigned j = 0; j < mincol; j++) {
                    temp[j][i] = dv[i][j];
                }
            }
        }
    }
    return temp;
}

// ==============================================================================
// 此处定义的函数是为了控制台输出查看方便
// ==============================================================================
// VVCD==vector<vector<complex<double> > >
// VVD== vector<vector<double> >
// ---------------------------------
void PrintVVCD(vector<vector<complex<double> > >Data) {
    cout << "Data : " << endl;
    for (unsigned i = 0; i < Data.size(); i++) {
        cout << "Row[" << i << "]=\t";
        vector<complex<double> >v = Data[i];
        copy(v.begin(), v.end(), ostream_iterator<complex<double> >(cout,
                "\t"));
        cout << endl;
    }
}

// ---------------------------------
void PrintVCD(vector<vector<double> >Data) {
    cout << "Data : " << endl;
    for (unsigned i = 0; i < Data.size(); i++) {
        cout << "Row[" << i << "]=\t";
        vector<double>v = Data[i];
        copy(v.begin(), v.end(), ostream_iterator<double>(cout, "\t"));
        cout << endl;
    }
}
// ---------------------------------

// --------------------------------------------------------------------------------
#endif