📄 fretrans.cpp
字号:
#include "stdafx.h"
#include "GlobalApi.h"
#include "Cdib.h"
#include <math.h>
#include <direct.h>
#include <complex>
using namespace std;
// FOURBYTES就是用来计算离4最近的整倍数
#define FOURBYTES(bits) (((bits) + 31) / 32 * 4)
/**************************************************************************
* 文件名：FreTrans.cpp
*
* 正交变换API函数库：
*
* THREECROSS() - 将实对称矩阵化作三对角矩阵
* BSTQ() - 求三对角对称矩阵的特征值和特征向量
* FFT_1D() - 快速一维傅立叶变换
* IFFT_1D() - 快速一维傅立叶反变换
* IFFT_2D() - 快速二维傅立叶反变换
* DCT() - 离散余弦变换
* IDCT() - 离散余弦反变换
* WALSH() - 沃尔什－哈达玛变换
* IWALSH() - 沃尔什－哈达玛反变换
*
*
* DIBFFT_2D() - 图象的二维离散傅立叶快速变换
* DIBDFT_2D() - 图象的二维离散傅立叶变换
* DIBHOTELLING() - 图象的霍特林变换
* DIBDct() - 图像的离散余弦变换
* DIBWalsh() - 图像的沃尔什－哈达玛变换
*
************************************************************************
*/
/*************************************************************************
*
* \函数名称：
* FFT_1D()
*
* \输入参数:
* complex<double> * pCTData - 指向时域数据的指针，输入的需要变换的数据
* complex<double> * pCFData - 指向频域数据的指针，输出的经过变换的数据
* nLevel －傅立叶变换蝶形算法的级数，2的幂数，
*
* \返回值:
* 无
*
* \说明:
* 一维快速傅立叶变换。
*
*************************************************************************
*/
VOID WINAPI FFT_1D(complex<double> * pCTData, complex<double> * pCFData, int nLevel)
{
// 循环控制变量
int i;
int j;
int k;
double PI = 3.1415926;
// 傅立叶变换点数
int nCount =0 ;
// 计算傅立叶变换点数
nCount =(int)pow(2,nLevel) ;
// 某一级的长度
int nBtFlyLen;
nBtFlyLen = 0 ;
// 变换系数的角度 ＝2 * PI * i / nCount
double dAngle;
complex<double> *pCW ;
// 分配内存，存储傅立叶变化需要的系数表
pCW = new complex<double>[nCount / 2];
// 计算傅立叶变换的系数
for(i = 0; i < nCount / 2; i++)
{
dAngle = -2 * PI * i / nCount;
pCW[i] = complex<double> ( cos(dAngle), sin(dAngle) );
}
// 变换需要的工作空间
complex<double> *pCWork1,*pCWork2;
// 分配工作空间
pCWork1 = new complex<double>[nCount];
pCWork2 = new complex<double>[nCount];
// 临时变量
complex<double> *pCTmp;
// 初始化，写入数据
memcpy(pCWork1, pCTData, sizeof(complex<double>) * nCount);
// 临时变量
int nInter;
nInter = 0;
// 蝶形算法进行快速傅立叶变换
for(k = 0; k < nLevel; k++)
{
for(j = 0; j < (int)pow(2,k); j++)
{
//计算长度
nBtFlyLen = (int)pow( 2,(nLevel-k) );
//倒序重排，加权计算
for(i = 0; i < nBtFlyLen/2; i++)
{
nInter = j * nBtFlyLen;
pCWork2[i + nInter] =
pCWork1[i + nInter] + pCWork1[i + nInter + nBtFlyLen / 2];
pCWork2[i + nInter + nBtFlyLen / 2] =
(pCWork1[i + nInter] - pCWork1[i + nInter + nBtFlyLen / 2])
* pCW[(int)(i * pow(2,k))];
}
}
// 交换 pCWork1和pCWork2的数据
pCTmp = pCWork1 ;
pCWork1 = pCWork2 ;
pCWork2 = pCTmp ;
}
// 重新排序
for(j = 0; j < nCount; j++)
{
nInter = 0;
for(i = 0; i < nLevel; i++)
{
if ( j&(1<<i) )
{
nInter += 1<<(nLevel-i-1);
}
}
pCFData[j]=pCWork1[nInter];
}
// 释放内存空间
delete pCW;
delete pCWork1;
delete pCWork2;
pCW = NULL ;
pCWork1 = NULL ;
pCWork2 = NULL ;
}
/*************************************************************************
*
* \函数名称：
* IFFT_1D()
*
* \输入参数:
* complex<double> * pCTData - 指向时域数据的指针，输入的需要反变换的数据
* complex<double> * pCFData - 指向频域数据的指针，输出的经过反变换的数据
* nLevel －傅立叶变换蝶形算法的级数，2的幂数，
*
* \返回值:
* 无
*
* \说明:
* 一维快速傅立叶反变换。
*
************************************************************************
*/
VOID WINAPI IFFT_1D(complex<double> * pCFData, complex<double> * pCTData, int nLevel)
{
// 循环控制变量
int i;
// 傅立叶反变换点数
int nCount;
// 计算傅立叶变换点数
nCount = (int)pow(2,nLevel) ;
// 变换需要的工作空间
complex<double> *pCWork;
// 分配工作空间
pCWork = new complex<double>[nCount];
// 将需要反变换的数据写入工作空间pCWork
memcpy(pCWork, pCFData, sizeof(complex<double>) * nCount);
// 为了利用傅立叶正变换,可以把傅立叶频域的数据取共轭
// 然后直接利用正变换，输出结果就是傅立叶反变换结果的共轭
for(i = 0; i < nCount; i++)
{
pCWork[i] = complex<double> (pCWork[i].real(), -pCWork[i].imag());
}
// 调用快速傅立叶变换实现反变换，结果存储在pCTData中
FFT_1D(pCWork, pCTData, nLevel);
// 求时域点的共轭，求得最终结果
// 根据傅立叶变换原理，利用这样的方法求得的结果和实际的时域数据
// 相差一个系数nCount
for(i = 0; i < nCount; i++)
{
pCTData[i] =
complex<double> (pCTData[i].real() / nCount, -pCTData[i].imag() / nCount);
}
// 释放内存
delete pCWork;
pCWork = NULL;
}
/*************************************************************************
*
* \函数名称：
* FFT_2D()
*
* \输入参数:
* complex<double> * pCTData - 图像数据
* int nWidth - 数据宽度
* int nHeight - 数据高度
* complex<double> * pCFData - 傅立叶变换后的结果
*
* \返回值:
* 无
*
* \说明:
* 二维傅立叶变换。
*
************************************************************************
*/
VOID WINAPI DIBFFT_2D(complex<double> * pCTData, int nWidth, int nHeight, complex<double> * pCFData)
{
// 循环控制变量
int x;
int y;
// 临时变量
double dTmpOne;
double dTmpTwo;
// 进行傅立叶变换的宽度和高度，（2的整数次幂）
// 图像的宽度和高度不一定为2的整数次幂
int nTransWidth;
int nTransHeight;
// 计算进行傅立叶变换的宽度 （2的整数次幂）
dTmpOne = log(nWidth)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransWidth = (int) dTmpTwo ;
// 计算进行傅立叶变换的高度 （2的整数次幂）
dTmpOne = log(nHeight)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransHeight = (int) dTmpTwo ;
// x，y（行列）方向上的迭代次数
int nXLev;
int nYLev;
// 计算x，y（行列）方向上的迭代次数
nXLev = (int) ( log(nTransWidth)/log(2) + 0.5 );
nYLev = (int) ( log(nTransHeight)/log(2) + 0.5 );
for(y = 0; y < nTransHeight; y++)
{
// x方向进行快速傅立叶变换
FFT_1D(&pCTData[nTransWidth * y], &pCFData[nTransWidth * y], nXLev);
}
// pCFData中目前存储了pCTData经过行变换的结果
// 为了直接利用FFT_1D，需要把pCFData的二维数据转置，再一次利用FFT_1D进行
// 傅立叶行变换（实际上相当于对列进行傅立叶变换）
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTData[nTransHeight * x + y] = pCFData[nTransWidth * y + x];
}
}
for(x = 0; x < nTransWidth; x++)
{
// 对x方向进行快速傅立叶变换，实际上相当于对原来的图象数据进行列方向的
// 傅立叶变换
FFT_1D(&pCTData[x * nTransHeight], &pCFData[x * nTransHeight], nYLev);
}
// pCFData中目前存储了pCTData经过二维傅立叶变换的结果，但是为了方便列方向
// 的傅立叶变换，对其进行了转置，现在把结果转置回来
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTData[nTransWidth * y + x] = pCFData[nTransHeight * x + y];
}
}
memcpy(pCTData, pCFData, sizeof(complex<double>) * nTransHeight * nTransWidth );
}
/*************************************************************************
*
* \函数名称：
* IFFT_2D()
*
* \输入参数:
* complex<double> * pCFData - 频域数据
* complex<double> * pCTData - 时域数据
* int nWidth - 图像数据宽度
* int nHeight - 图像数据高度
*
* \返回值:
* 无
*
* \说明:
* 二维傅立叶反变换。
*
************************************************************************
*/
VOID WINAPI IFFT_2D(complex<double> * pCFData, complex<double> * pCTData, int nWidth, int nHeight)
{
// 循环控制变量
int x;
int y;
// 临时变量
double dTmpOne;
double dTmpTwo;
// 进行傅立叶变换的宽度和高度，（2的整数次幂）
// 图像的宽度和高度不一定为2的整数次幂
int nTransWidth;
int nTransHeight;
// 计算进行傅立叶变换的宽度 （2的整数次幂）
dTmpOne = log(nWidth)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransWidth = (int) dTmpTwo ;
// 计算进行傅立叶变换的高度 （2的整数次幂）
dTmpOne = log(nHeight)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransHeight = (int) dTmpTwo ;
// 分配工作需要的内存空间
complex<double> *pCWork= new complex<double>[nTransWidth * nTransHeight];
//临时变量
complex<double> *pCTmp ;
// 为了利用傅立叶正变换,可以把傅立叶频域的数据取共轭
// 然后直接利用正变换，输出结果就是傅立叶反变换结果的共轭
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTmp = &pCFData[nTransWidth * y + x] ;
pCWork[nTransWidth * y + x] = complex<double>( pCTmp->real() , -pCTmp->imag() );
}
}
// 调用傅立叶正变换
::DIBFFT_2D(pCWork, nWidth, nHeight, pCTData) ;
// 求时域点的共轭，求得最终结果
// 根据傅立叶变换原理，利用这样的方法求得的结果和实际的时域数据
// 相差一个系数
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTmp = &pCTData[nTransWidth * y + x] ;
pCTData[nTransWidth * y + x] =
complex<double>( pCTmp->real()/(nTransWidth*nTransHeight),
-pCTmp->imag()/(nTransWidth*nTransHeight) );
}
}
delete pCWork ;
pCWork = NULL ;
}
/*************************************************************************
*
* 函数名称：
* DCT()
*
* 参数:
* double * dpf - 指向时域值的指针
* double * dpF - 指向频域值的指针
* r －2的幂数
*
* 返回值:
* 无。
*
* 说明:
* 实现一维快速离散余弦变换。
*
***********************************************************************
*/
VOID WINAPI DCT(double *dpf, double *dpF, int r)
{
double PI = 3.1415926;
// 离散余弦变换点数
LONG lNum;
// 循环变量
int i;
// 中间变量
double dTemp;
complex<double> *comX;
// 离散余弦变换点数
lNum = 1<<r;
// 分配内存
comX = new complex<double>[lNum*2];
// 赋初值0
memset(comX, 0, sizeof(complex<double>) * lNum * 2);
// 将时域点转化为复数形式
for(i=0;i<lNum;i++)
{
comX[i] = complex<double> (dpf[i], 0);
}
// 调用快速付立叶变换
FFT_1D(comX,comX,r+1);
// 调整系数
dTemp = 1/sqrt(lNum);
// 求dpF[0]
dpF[0] = comX[0].real() * dTemp;
dTemp *= sqrt(2);
// 求dpF[u]
for(i = 1; i < lNum; i++)
{
dpF[i]=(comX[i].real() * cos(i*PI/(lNum*2))
+ comX[i].imag() * sin(i*PI/(lNum*2))) * dTemp;
}
// 释放内存
delete comX;
}
/*************************************************************************
*
* 函数名称：
* IDCT()
*
* 参数:
* double * dpF - 指向频域值的指针
* double * dpf - 指向时域值的指针
* r －2的幂数
*
* 返回值:
* 无。
*
* 说明:
* 实现一维快速离散余弦反变换。
*
************************************************************************/
VOID WINAPI IDCT(double *dpF, double *dpf, int r)
{
double PI = 3.1415926;
// 离散余弦反变换点数
LONG lNum;
// 计算离散余弦变换点数
lNum = 1<<r;
// 循环变量
int i;
// 中间变量
double dTemp, d0;
complex<double> *comX;
// 给复数变量分配内存
comX = new complex<double>[lNum*2];
// 赋初值0
memset(comX, 0, sizeof(complex<double>) * lNum * 2);
// 将频域变换后点写入数组comX
for(i=0;i<lNum;i++)
{
comX[i] = complex<double> (cos(i*PI/(lNum*2)) * dpF[i], sin(i*PI/(lNum*2) * dpF[i]));
}
// 调用快速付立叶反变换
IFFT_1D(comX,comX,r+1);
// 调整系数
dTemp = sqrt(2.0/lNum);
d0 = (sqrt(1.0/lNum) - dTemp) * dpF[0];
// 计算dpf(x)
for(i = 0; i < lNum; i++)
{
dpf[i] = d0 + 2 * lNum * comX[i].real()* dTemp ;
}
// 释放内存
delete comX;
}
/*************************************************************************
*
* 函数名称：
* WALSH()
*
* 参数:
* double * dpf - 指向时域值的指针
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -