?? 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 + -