📄 fouriertransforms.cpp
字号:
/************************************************************
Copyright (C), 2008, LusterLightVision.
FileName: FourierTransforms.cpp
作 者: 赵敏
版 本: 1.00
日 期: 08/06/20
模块描述:
该模块为《付立叶变换》的代码实现。
主要函数及其功能:
FFT_1D() 一维付立叶正变换
IFFT_1D() 一维付立叶逆变换
FFT_2D() 二维付立叶变换
IFFT_2D() 二维傅立叶反变换
DIBFFT_2D 图像傅立叶变换
修改历史记录:
<编号> <作者> <时间> <版本> <描述>
1 赵敏 08/06/20 1.00 创建该模块，实现基本功能
***********************************************************/
#include "stdafx.h"
#include "LLVFourierTransforms.h"
#include "math.h"
#include <xmmintrin.h>
#include <emmintrin.h>
#define xmalloc(s) ::HeapAlloc(::GetProcessHeap(), HEAP_ZERO_MEMORY, (s)) //Dll中分配内存
#define xfree(p) ::HeapFree (::GetProcessHeap(), 0, (p)) //调用模块中释放内存
/*************************************************************************
*
* 函数名称：
* FFT_1D()
*
* 参数:
* complex<double> * pCTData - 指向时域数组的指针
* complex<double> * pCFData - 指向频域数组的指针
* nLevel －付立叶变换蝶形算法的级数，2的幂数
*
* 返回值:
* 无。
*
* 说明:
* 该函数用来实现快速付立叶变换。
*
************************************************************************/
LLVPRO_API void FFT_1D(complex<double> * pCTData, complex<double> * pCFData, int nLevel)
{
// 付立叶变换点数
long lcount;
// 循环变量
int i,j,k;
int p;
// 某一级长度
int nBtFlyLen = 0;
// 角度
double dAngle;
complex<double> *pCW,*pCWork1,*pCWork2,*pCTmp;
// 计算付立叶变换点数
lcount = 1 << nLevel;
// 分配运算所需存储器
pCW = new complex<double>[lcount / 2];
pCWork1 = new complex<double>[lcount];
pCWork2 = new complex<double>[lcount];
double dCof = PI * 2 / lcount;
// 计算加权系数
for(i = 0; i < lcount / 2; i++)
{
dAngle = -i * dCof;
pCW[i] = complex<double> (cos(dAngle), sin(dAngle));
}
// 将时域点写入pCWork1
memcpy(pCWork1, pCTData, sizeof(complex<double>) * lcount);
int nBtFlyLenHalf;
// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < nLevel; k++)
{
for(j = 0; j < 1 << k; j++)
{
// 计算长度
nBtFlyLen = 1 << (nLevel-k);
nBtFlyLenHalf = nBtFlyLen / 2;
for(i = 0; i < nBtFlyLenHalf; i++)
{
p = j * nBtFlyLen;
int nIndex1 = i + p;
int nIndex2 = nIndex1 +nBtFlyLenHalf;
pCWork2[nIndex1] = pCWork1[nIndex1] + pCWork1[nIndex2];
pCWork2[nIndex2] = (pCWork1[nIndex1] - pCWork1[nIndex2]) * pCW[i * (1<<k)];
}
}
// 交换pCWork1和pCWork2的数据
pCTmp = pCWork1;
pCWork1 = pCWork2;
pCWork2 = pCTmp;
}
// 变换后序列是码字倒序排列的，需要重新排序
for(j = 0; j < lcount; j++)
{
p = 0;
for(i = 0; i < nLevel; i++)
{
if (j&(1<<i))
{
p+=1<<(nLevel-i-1);
}
}
pCFData[j] = pCWork1[p];
}
// 释放内存
delete []pCW;
delete []pCWork1;
delete []pCWork2;
pCW = NULL ;
pCWork1 = NULL ;
pCWork2 = NULL ;
}
LLVPRO_API void FFT_1D_Opti(complex<double> * pCTData, complex<double> * pCFData, int nLevel)
{
// 付立叶变换点数
long lcount;
// 循环变量
int i,j,k;
// 某一级长度
int nBtFlyLen = 0;
// 角度
double dAngle;
complex<double> *pCW,*pCWork1,*pCWork2,*pCTmp;
// 计算付立叶变换点数
lcount = 1 << nLevel;
// 分配运算所需存储器
pCW = new complex<double>[lcount / 2];
pCWork1 = new complex<double>[lcount];
pCWork2 = new complex<double>[lcount];
double dCof = PI * 2 / lcount;
// 计算加权系数
for(i = 0; i < lcount / 2; i++)
{
dAngle = -i * dCof;
double dcosA = cos(dAngle);
double dsinA = sqrt(1- dcosA*dcosA);
pCW[i] = complex<double> (dcosA, -dsinA);
}
// 将时域点写入pCWork1
memcpy(pCWork1, pCTData, sizeof(complex<double>) * lcount);
int nBtFlyLenHalf,nIndex1,nIndex2,p;
__m128d md1,md2,md3,md4,md5,md6;
// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < nLevel; k++)
{
for(j = 0; j < 1 << k; j++)
{
// 计算长度
nBtFlyLenHalf = 1 << (nLevel-k-1);
nBtFlyLen = 1 << (nLevel-k);
p = j * nBtFlyLen;
for(i = 0; i < (nBtFlyLenHalf&~7) ; i+=8)
{
/*
nIndex1 = i + p;
nIndex2 = nIndex1 + nBtFlyLenHalf;
md1 = _mm_set_pd(pCWork1[nIndex1].imag(),pCWork1[nIndex1].real()); // 读取数据
md2 = _mm_set_pd(pCWork1[nIndex2].imag(),pCWork1[nIndex2].real()); // 读取数据
md3 = _mm_add_pd(md1,md2); // pCWork1[i + p] + pCWork1[i + p +nBtFlyLenHalf ];
md4 = _mm_sub_pd(md1,md2); // pCWork1[i + p] - pCWork1[i + p +nBtFlyLenHalf ];
int n = i <<k;
md1 = _mm_set_pd(pCW[n].real(),pCW[n].real());
md5 = _mm_mul_pd(md4,md1);
md1 = _mm_set_pd(pCW[n].imag(),pCW[n].imag());
md4 = _mm_shuffle_pd(md4,md4,1);
md4 = _mm_mul_pd(md4,md1);
md6 = _mm_add_pd(md5,md4);
md5 = _mm_sub_pd(md5,md4);
md6 = _mm_shuffle_pd(md5,md6,2);
double dWorks[2];
_mm_storeu_pd(dWorks,md3);
pCWork2[nIndex1] = complex<double>(dWorks[0], dWorks[1]);
_mm_storeu_pd(dWorks,md6);
pCWork2[nIndex2] = complex<double>(dWorks[0], dWorks[1]);
//*/
nIndex1 = i + p;
nIndex2 = i + p + nBtFlyLenHalf;
pCWork2[nIndex1] = pCWork1[nIndex1] + pCWork1[nIndex2];
pCWork2[nIndex2] = (pCWork1[nIndex1] - pCWork1[nIndex2]) * pCW[i <<k];
nIndex1 = i + p + 1;
nIndex2 = i + p + nBtFlyLenHalf + 1;
pCWork2[nIndex1] = pCWork1[nIndex1] + pCWork1[nIndex2];
pCWork2[nIndex2] = (pCWork1[nIndex1] - pCWork1[nIndex2]) * pCW[(i+1) <<k];
nIndex1 = i + p + 2;
nIndex2 = i + p + nBtFlyLenHalf + 2;
pCWork2[nIndex1] = pCWork1[nIndex1] + pCWork1[nIndex2];
pCWork2[nIndex2] = (pCWork1[nIndex1] - pCWork1[nIndex2]) * pCW[(i+2) <<k];
nIndex1 = i + p + 3;
nIndex2 = i + p + nBtFlyLenHalf + 3;
pCWork2[nIndex1] = pCWork1[nIndex1] + pCWork1[nIndex2];
pCWork2[nIndex2] = (pCWork1[nIndex1] - pCWork1[nIndex2]) * pCW[(i+3) <<k];
nIndex1 = i + p + 4;
nIndex2 = i + p + nBtFlyLenHalf + 4;
pCWork2[nIndex1] = pCWork1[nIndex1] + pCWork1[nIndex2];
pCWork2[nIndex2] = (pCWork1[nIndex1] - pCWork1[nIndex2]) * pCW[(i+4) <<k];
nIndex1 = i + p + 5;
nIndex2 = i + p + nBtFlyLenHalf + 5;
pCWork2[nIndex1] = pCWork1[nIndex1] + pCWork1[nIndex2];
pCWork2[nIndex2] = (pCWork1[nIndex1] - pCWork1[nIndex2]) * pCW[(i+5) <<k];
nIndex1 = i + p + 6;
nIndex2 = i + p + nBtFlyLenHalf + 6;
pCWork2[nIndex1] = pCWork1[nIndex1] + pCWork1[nIndex2];
pCWork2[nIndex2] = (pCWork1[nIndex1] - pCWork1[nIndex2]) * pCW[(i+6) <<k];
nIndex1 = i + p + 7;
nIndex2 = i + p + nBtFlyLenHalf + 7;
pCWork2[nIndex1] = pCWork1[nIndex1] + pCWork1[nIndex2];
pCWork2[nIndex2] = (pCWork1[nIndex1] - pCWork1[nIndex2]) * pCW[(i+7) <<k];
}
for(i = (nBtFlyLenHalf&~7); i < nBtFlyLenHalf ; i++)
{
nIndex1 = i + p;
nIndex2 = nIndex1 + nBtFlyLenHalf;
pCWork2[nIndex1] = pCWork1[nIndex1] + pCWork1[nIndex2];
pCWork2[nIndex2] = (pCWork1[nIndex1] - pCWork1[nIndex2]) * pCW[i <<k];
}
}
// 交换pCWork1和pCWork2的数据
pCTmp = pCWork1;
pCWork1 = pCWork2;
pCWork2 = pCTmp;
}
// 变换后序列是码字倒序排列的，需要重新排序
__m128i m0,m1,m2,m,m3,m4,m5,mp;
_declspec(align(16)) int nIndex[4] = {0,0,0,0} ;
m = _mm_set1_epi32(1);
m0 = _mm_setzero_si128();
for(j = 0; j < lcount; j++)
{
mp = _mm_set1_epi32(0);
m1 = _mm_set1_epi32(j);
for(i = 0; i < (nLevel&~3); i +=4)
{
// i
m2 = _mm_slli_epi32(m,i);
m2 = _mm_and_si128(m1,m2);
m2 = _mm_cmpgt_epi32(m2,m0);
m3 = _mm_slli_epi32(m,nLevel-i-1);
m2 = _mm_and_si128(m3,m2);
mp = _mm_add_epi32(mp,m2);
// i+1
m2 = _mm_slli_epi32(m,i+1);
m2 = _mm_and_si128(m1,m2);
m2 = _mm_cmpgt_epi32(m2,m0);
m3 = _mm_slli_epi32(m,nLevel-i-2);
m2 = _mm_and_si128(m3,m2);
mp = _mm_add_epi32(mp,m2);
// i+2
m2 = _mm_slli_epi32(m,i+2);
m2 = _mm_and_si128(m1,m2);
m2 = _mm_cmpgt_epi32(m2,m0);
m3 = _mm_slli_epi32(m,nLevel-i-3);
m2 = _mm_and_si128(m3,m2);
mp = _mm_add_epi32(mp,m2);
// i+3
m2 = _mm_slli_epi32(m,i+3);
m2 = _mm_and_si128(m1,m2);
m2 = _mm_cmpgt_epi32(m2,m0);
m3 = _mm_slli_epi32(m,nLevel-i-4);
m2 = _mm_and_si128(m3,m2);
mp = _mm_add_epi32(mp,m2);
}
_mm_store_si128((__m128i*)nIndex, mp);
p = nIndex[0];
for(i=nLevel & ~3; i<nLevel;i++)
{
if (j&(1<<i))
{
p+=1<<(nLevel-i-1);
}
}
pCFData[j] = pCWork1[p];
}
// 释放内存
delete []pCW;
delete []pCWork1;
delete []pCWork2;
pCW = NULL ;
pCWork1 = NULL ;
pCWork2 = NULL ;
}
/*************************************************************************
*
* 函数名称：
* IFFT_1D()
*
* 参数:
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -