问题描述
- 二维傅里叶变换转一维傅里叶变换
- 看了一个论文说二维傅里叶变换图谱经圆积分计算后可以得到一维傅里叶变换，求大神推导公式或者过程
解决方案
// ==============================================================================// 快速离散傅里叶变换和功率谱// 一维快速傅里叶变换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 << ""]= ""; vector<complex<double> >v = Data[i]; copy(v.begin() v.end() ostream_iterator<complex<double> >(cout "")); cout << endl; }}// ---------------------------------void PrintVCD(vector<vector<double> >Data) { cout << ""Data : "" << endl; for (unsigned i = 0; i < Data.size(); i++) { cout << ""Row["" << i << ""]= ""; vector<double>v = Data[i]; copy(v.begin() v.end() ostream_iterator<double>(cout "")); cout << endl; }}// ---------------------------------// --------------------------------------------------------------------------------#endif解决方案二：
傅里叶变换（转）
什么是傅里叶变换
小波变换与傅里叶变换
解决方案三：
DSP好像有具体的求解证明过程。
时间： 2024-04-25 02:25:17