Matlab 和 C++ （4），FFT（1）
有空整理了一下关于FFT的算法。目的就是希望独立于Matlab，也不想用人家的不知道的代码。FFT 现在最流行的就是蝶形算法，其理论基础就是，一个序列的FFT的变换，可以用原来序列的偶数列的变换，和奇数列的变换，相组合。
用下面的 Marlab 运算可以检验结果：
>>N = 1024;
>>A = rand( N, 1 ) + rand( N, 1 ) * i;
>>FA1 =fft(A); %get fft results
这里，我们得到了A的fft变换的结果。
下面我们用蝶形算法：
先得到序列的偶数和奇数项（注意，一定要从0开始算下标，因为历史原因，DFT在推算的时候，第一项用的是0，所以算为序列的偶数项）
>>A_Even = A( ( 1: N/2 ) * 2 - 1 ); % even start from the first number!
>>A_Odd = A( ( 1: N/2 ) * 2 );
>>fA_Even= fft( A_Even );
>>fA_Odd = fft( A_Odd );
到这里，我们有了奇数项和偶数项的变换，于是总的变换为：
>>N = 2*pi*( 0: N/2-1 ) /N;
>>FA2 = [ fA_Even + ( cos( N ) + i*sin( N ) )' .* fA_Odd ; fA_Even - ( cos( N ) + i*sin( N ) )' .* fA_Odd ];
怎么样，FA2 和 FA1 一样！
因为计算fft的时候，需要计算很多的sin，cos来得到 复指数值；一个长的序列如果用两个一半的序列代替，而相同的序列意味着这些sin，cos计算只要运算一次，所以，就节约了很多时间。
如果序列长度是2的指数次，这种对分可以一直进行，一直到最有只有两个数。
这就是蝶形算法。由 J.W. Cooley 和 J.W Tukey 在60年代中发明。
下面是我们的c++实现的一维的fft源程序，为了简单期间，我们规定的输入长度是2的指数次，并且，没有为此作检验。
文件名是myfft.cpp
编译后可以这样用：
>>N = 1024;
>>A = rand( N, 1 ) + rand( N, 1 ) * i;
>>FA =myfft(A);
-
//---------------------------------------------------------------------------
//--------------myfft.cpp
//---------------------------------------------------------------------------
#include <math.h>
#include "mex.h"
#define MY_LCC
//---------------------------------------------------------------------------
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
#define My_PI 3.141592653589793238462643383279502884
#define My_2_PI 6.283185307179586476925286766559005768
#define mSIN mysin_sf
#define mCOS mycos_sf
//#define mSIN sin
//#define mCOS cos
//---------------------------------------------------------------------------
#ifndef MY_LCC
double __declspec(naked) mysin_sf(double __a){
asm{
fld
fsin
ret
};
}
double __declspec(naked) mycos_sf(double __a){
asm{
fld
fcos
ret
};
}
#else
double mysin_sf(double __a)
{
//register double __result;
_asm("fldl 4(%esp)"); // a
_asm("fsin"); // sin(a)
// _asm("fstp (%__result)");
_asm("exit1:");
_asm("ret");
//return __result;// return a long double
}
double mycos_sf(double __a)
{
//register double __result;
_asm("fldl 4(%esp)"); // a
_asm("fcos"); // cos(a)
// _asm("fstp (%__result)");
_asm("exit2:");
_asm("ret");
//return __result;// return a long double
}
#endif
//---------------------------------------------------------------------------
void mySplit_C(double *pData, int nCount){
int i, j, m, n = (nCount << 1)-1;
double dTemp, *pDataI = pData + 1;
for (i=j=0;i<n;i+=2) {
if (j > i) {
dTemp = pData; pData = pData; pData = dTemp;
dTemp = pDataI; pDataI = pDataI; pDataI = dTemp;
}
m=nCount;
while (m >= 2 && j >= m) {
j -= m;
m >>= 1;
}
j += m;
}
}
//---------------------------------------------------------------------------
void myFFT_C(double *pData, int nCount){
int i, j, m, n, nMax, nStep;
double dTheta;
double dWC_r, dWC_i, dWK_r, dWK_i, dTempR, dTempI;
double *pDataI = pData + 1;
n = nCount << 1;
nMax=2;
while (n > nMax) {
nStep = nMax << 1;
dTheta= My_2_PI/nMax;
dWK_r = mCOS( dTheta );
dWK_i = mSIN( dTheta );
dWC_r = 1;
dWC_i = 0;
for(m=0; m<nMax; m+=2) { //Here are the two nested inner loops.
for (i=m; i<n; i+=nStep) {
j=i+nMax;
dTempR = dWC_r * pData - dWC_i* pDataI;
dTempI = dWC_r * pDataI + dWC_i* pData;
pData = pData - dTempR;
pDataI = pDataI - dTempI;
pData += dTempR;
pDataI += dTempI;
}
dWC_r = (dTempR = dWC_r) * dWK_r - dWC_i * dWK_i;
dWC_i = dWC_i * dWK_r + dTempR * dWK_i;
}
nMax=nStep;
}
}
//---------------------------------------------------------------------------
void fft_1(double *pData, int nCount, int nSign){
mySplit_C(pData, nCount);
myFFT_C(pData, nCount);
}
//---------------------------------------------------------------------------
void real_fft_1(double *pData, int nCount, int nSign){
int i, i1, i2, i3, i4, np3, ni;
double c1, c2, h1r, h1i, h2r, h2i;
double wr,wi,wpr,wpi,wtemp,theta;
theta= My_PI/(double)(nCount);
if (nSign == 1) {
fft_1(pData, nCount >>1, 1);
} else {
theta = -theta;
}
wtemp=sin(theta);
wpr = -2.0*wtemp*wtemp+1;
wpi=sin(theta + theta);
wr=wpr;
wi=wpi;
np3= nCount + 1; //+3
ni = nCount >> 2;
for (i=2; i<= ni; i++) { //for (i=2; i<= ni; i++) {
i4= 1+ ( i3 = np3 - ( i2= 1+ (i1= i+i-2)));
h1r=0.5*( (c1 = pData) + (c2 = pData) );
h2i=0.5*(c1 - c2);
h1i=0.5*((c1=pData) - (c2=pData));
h2r=0.5*(c1 + c2);
pData=h1r+ (c1= wr*h2r + wi*h2i );
pData=h1r- c1;
pData=h1i+ (c1 = wi*h2r - wr*h2i ) ;
pData=c1 - h1i;
//wr=(wtemp=wr)*wpr-wi*wpi+wr;
//wi=wi*wpr+wtemp*wpi+wi;
wr=(wtemp=wr)*wpr-wi*wpi;
wi=wi*wpr+wtemp*wpi;
}
if (nSign == 1) {
pData = (h1r=pData) + pData;
pData = h1r - pData;
} else {
pData=0.5* ((h1r=pData) + pData);
pData=0.5* (h1r-pData);
fft_1(pData, nCount >> 1, -1);
}
}
//---------------------------------------------------------------------------
#ifndef MY_LCC
//---------------------------------------------------------------------------
//#pragma argsused
//int WINAPI DllEntryPoint(HINSTANCE hinst, unsigned long reason, void* lpReserved)
//{
// return 1;
//}
//---------------------------------------------------------------------------
void _exportmexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
#else
voidmexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
#endif
{
int i, j, nArrayID;
double *xData, *xDataC, *pTemp, *pTempR, *pTempC;
double *pTempRp, *pTempCp;
bool bReal;
if (1!=nrhs ) mexErrMsgTxt("Using:= MYFFT(X); X is input data");
j = mxGetN(prhs) * mxGetM(prhs);
xData = (double *) mxGetData (prhs);
pTempC =(double *) mxGetPi (prhs);
nArrayID =mxGetClassID(prhs);
nLine=0;
if(nArrayID==mxDOUBLE_CLASS ){
pTemp = xDataC = (double *) malloc(j*2*sizeof(double));
if(pTempC==NULL){
bReal=true;
for(i=j;i>0;i--){ *pTemp++ = *xData++;}
real_fft_1(xDataC, j, 1);
}else{
bReal=false;
for(i=j;i>0;i--){ *pTemp++ = *xData++; *pTemp++ = *pTempC++;}
fft_1(xDataC, j, 1);
}
}else mexErrMsgTxt("X is unknow type !");
plhs = mxCreateDoubleMatrix(j , 1, mxCOMPLEX);
pTempR = (double *)mxGetPr(plhs);
pTempC = (double *)mxGetPi(plhs);
pTemp= xDataC;
if(bReal){
pTempRp = pTempR + j- 1;
pTempCp = pTempC + j- 1;
*pTempR = *pTemp++;
pTempR=*pTemp++;
*pTempC = pTempC = 0;
for(i=(j>>1)-1;i>0;i--){
(*pTempRp--) = (*(++pTempR)) = (*pTemp++);
(*pTempCp--) = -((*(++pTempC)) = (*pTemp++));
}
}else{
for(i=j;i>0;i--){
*pTempR++ = *pTemp++;
*pTempC++ = *pTemp++;
}
}
free(xDataC);
}
//---------------------------------------------------------------------------
这里有几点说明：
1）matlab用的是fftw.dll，计算fft用了exp(-i * 2*pi*N)，我们的程序因为原来的目的是处理2维的图像数据，为了修正数据排列，我们用exp(i * 2*pi*N)；所以，结果和matlab镜像对称。
2）速度在一维的时候，没有fftw.dll快，虽然算法和人家的是一样的，但是fftw用了多线程和汇编优化，这里我们只对sin cos进行了一定的优化，倒序算法和fft都还是c本身。以后我们会给一个多线程的版本，速度完全可以和fftw一比。
3）二维以上的算法由一维得到，因为我们的图像数据是unsigned int 的灰度。所以，用我们自己的程序比matlab fft2要快许多。（二维的fft算法以后会陆续给出） 假如用BC和VC，倒序算法可以这样优化：
void __declspec(naked) mySplit_A(double *pData, int nCount){//for intel single and AMD
asm{
push ebp
mov ebp,esp
add esp,-12
push ebx
push esi
xor ecx,ecx
mov eax,dword ptr
add eax,eax
dec eax
mov dword ptr ,eax
mov edx,dword ptr
add edx,8
mov dword ptr ,edx
mov edx,ecx
cmp ecx,dword ptr
jge short my@16
// ; EDX = j, ECX = i, EBX = @temp3
my@15:
cmp ecx,edx
jge short my@17
mov esi, dword ptr //pData
mov eax, dword ptr // eax = pData
mov ebx, dword ptr // ebx = pData
mov dword ptr , eax // pData = eax
mov dword ptr , ebx // pData = ebx
add esi, 4
mov eax, dword ptr // eax = pData
mov ebx, dword ptr // ebx = pData
mov dword ptr , eax // pData = eax
mov dword ptr , ebx // pData = ebx
mov esi, dword ptr //pDataI
mov eax, dword ptr // eax = pData
mov ebx, dword ptr // ebx = pData
mov dword ptr , eax // pData = eax
mov dword ptr , ebx // pData = ebx
add esi, 4
mov eax, dword ptr // eax = pData
mov ebx, dword ptr // ebx = pData
mov dword ptr , eax // pData = eax
mov dword ptr , ebx // pData = ebx
my@17:
mov eax,dword ptr
jmp short my@19
// ; EAX = m, EDX = j, ECX = i, EBX = @temp3
my@18:
sub edx,eax
sar eax,1
my@19:
cmp eax,2
jl short my@20
cmp eax,edx
jle short my@18
my@20:
add edx,eax
add ecx,2
// add dword ptr ,16
// add ebx,16
cmp ecx,dword ptr
jl short my@15
//?live16388@224: ;
my@16:
//@22:
pop esi
pop ebx
mov esp, ebp
pop ebp
ret
}
//---------------------------------------------------------------------------
void __declspec(naked) mySplit_A2(double *pData, int nCount){// for intel duro core & AMD
asm{
push ebp
mov ebp,esp
add esp,-12
push ebx
push esi
xor ecx,ecx
mov eax,dword ptr
add eax,eax
dec eax
mov dword ptr ,eax
mov edx,dword ptr
add edx,8
mov dword ptr ,edx
mov edx,ecx
cmp ecx,dword ptr
jge short my@216
// ; EDX = j, ECX = i, EBX = @temp3
my@215:
cmp ecx,edx
jge short my@217
mov esi, dword ptr //pData
fwait
fld // st0 = pData
fld // st0 = pData; st1 = pData
fstp // st0 -> pData ; st0 = old pData
fstp
mov esi, dword ptr //pDataI
fld // st0 = pData
fld // st0 = pData; st1 = pData
fstp // st0 -> pData ; st0 = old pData
fstp
my@217:
mov eax,dword ptr
jmp short my@219
// ; EAX = m, EDX = j, ECX = i, EBX = @temp3
my@218:
sub edx,eax
sar eax,1
my@219:
cmp eax,2
jl short my@220
cmp eax,edx
jle short my@218
my@220:
add edx,eax
add ecx,2
// add dword ptr ,16
// add ebx,16
cmp ecx,dword ptr
jl short my@215
//?live16388@224: ;
my@216:
//@22:
pop esi
pop ebx
mov esp, ebp
pop ebp
ret
}
}
//---------------------------------------------------------------------------
不过，假如单线程时，或者，假如对很大的数据块进行操作时，即使你的cpu是双核的，
你也要在调用汇编倒序函数之前 和 mov之间，插入适当的 fwait 指令和 memory lock函数。
要不然，还是让VC 和 BC的编译器自己优化，不要用汇编。
因为windows xp和vesta里面，对大的连续内存管理的很抠。
加上matlab自己没有动态内存管理，你不加入 lock & wait or fwait，一旦你的过程运行超过一个时间段，就会被windows的内存管理接管，从而花费很多时间，在无端的总线等待上，尤其是执行到mov指令时，你的指令也许就一个时钟周期，你的wait可以最差到4个周期。
所以，我们提倡用多线程，然后，将自己的权限提高，这样就可以得到更多的cpu资源了。
这里虽然贴了优化，但是不提倡用。
[ 本帖最后由 recbio 于 2007-8-2 18:09 编辑 ] 办 证:Q.⑥②⑥.⑥④④.③②◎.◆◆◆◆◆◆◆Θ
办 证：Q.626.644.320【】δキゲΘセツδ【】
办 证：Q.626.644.320【】δキゲΘセツδ【】
办 证：q.626.644.320〓〓★★〓〓〓★★〓〓〓
& z/ O, [+ u% r% s- F办 证Q.626.644.320★★★★★★★★★★★ 快速办理英语四. 8 O) ]" m$ l- G: n六N级证，雅思，公共英语证，
希望我们能在各种英语证方面对您有所帮助。如有需要请随时和我们联系。很多委托我司办理办
! \1 Q- r+ ^6 x3 d# X, S0 O$ [真的证件的客户都非常满意。我司在办各种英语证业内有着良好的口碑。在办理英语证行业有 # b6
希望我们能在各种英语证方面对您有所帮助。如有需要请随时和我们联系。很多委托我司办理办
! \1 Q- r+ ^6 x3 d# X, S0 O$ [真的证件的客户都非常满
意。我司在办各种英语证业内有着良好的口碑。在办理英语证行业有 # b6 b, p3 N8 E( ^1 F8 a! q8 T
着多年的历史。英语证是我司的主营业务之一。我们将竭诚为您服务。我们在英语证上有着丰富
) A, q'' n5 f( Y# D9 ?) M$ Y的经验。如有需要请随时和我们联系q.626.644.320另外办理各种学历文凭、业务及根据客
; q6 e& ?* U+ T; ?; J/ }2 f户样品及要求制作一切证联系q.626.644.320。 快速办理英语证，雅思，公共英语证，
: w% ^4 ]! y. n'' i* D希望我们能在各种英语证方面对您有所帮助。如有需要请随时和我们联系。很多委托我司办理办
$ @" v4 |5 r6 I7 `) g3 V# S真的证件的客户都非常满意。我司在办各种英语证业内有着良好的
" M* _. K! B6 B, u% _各种英语证业内有着良好的口碑。在办理英语证行业有着多年的历史。英语证是我司的主营业务 , B4 e# Q! U, K5 g1 z3 t, U) }, A
之一。我们将竭诚为您服务。我们在英语证上有着丰富的经验。如有需要请随时和我们联系QQ " l9 k5 O0 f9 ?! e$ B: ~
q.626.644.320。另外办理各种学历文凭、业务及根据客户样品及要求制作一切证联系QQ
8 M'' a& [- ^! |& [, q.626.644.320。 快速办理英语证，雅思，公共英语证，希望我们能在各种英语证方面对您 & r'' i+ N0 V1 p8 T6 v: `8 o
有所帮助。如有需要请随时和我们联系。很多委托我司办理办真的证件的客户都非常满意。我
6 P9 A) d; P4 q, Y3 } C司在办各种英语证业内有着良好的口碑。在办理英语证行业有着多年的历史。英语证是我司的主 " b; \$ ?3 p7 @, h0 y1 N
营业务之一。我们将竭诚为您服务。我们在英语证上有着丰富的经验快速办理英语四. 六N级证， % M" J. y# y* a; d+ j
雅思，公共英语证，希望我们能在各种英语证方面对您有所帮助。如有需要请随时和我们联系。
: g0 @4 w" J; d5 }! }0 f2 `很多委托我司办理办真的证件的客户都非常满意。我司在办各种英语证业内有着良好的口碑。 - b# g/ \2 R: [2 D" i9 U8 G8 W
在办理英语证行业有着多年的历史。英语证是我司的主营业务之一。我们将竭诚为您服务
页:
[1]