Matlab 和 C++ （4），FFT（1）

recbio · 发表于 2007-8-2 16:42

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？注册

x

有空整理了一下关于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 [esp+4]
fsin
ret
};
}
double __declspec(naked) mycos_sf(double __a){
asm{
fld [esp+4]
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[j]; pData[j] = pData[i]; pData[i] = dTemp;
dTemp = pDataI[j]; pDataI[j] = pDataI[i]; pDataI[i] = 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[j] - dWC_i* pDataI[j];
dTempI = dWC_r * pDataI[j] + dWC_i* pData[j];
pData[j] = pData[i] - dTempR;
pDataI[j] = pDataI[i] - dTempI;
pData[i] += dTempR;
pDataI[i] += 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[i1]) + (c2 = pData[i3]) );
h2i=0.5*(c1 - c2);
h1i=0.5*((c1=pData[i2]) - (c2=pData[i4]));
h2r=0.5*(c1 + c2);
pData[i1]=h1r + (c1= wr*h2r + wi*h2i );
pData[i3]=h1r - c1;
pData[i2]=h1i + (c1 = wi*h2r - wr*h2i ) ;
pData[i4]=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[0] = (h1r=pData[0]) + pData[1];
pData[1] = h1r - pData[1];
} else {
pData[0]=0.5* ((h1r=pData[0]) + pData[1]);
pData[1]=0.5* (h1r-pData[1]);
fft_1(pData, nCount >> 1, -1);
}
}
//---------------------------------------------------------------------------
#ifndef MY_LCC
//---------------------------------------------------------------------------
//#pragma argsused
//int WINAPI DllEntryPoint(HINSTANCE hinst, unsigned long reason, void* lpReserved)
//{
// return 1;
//}
//---------------------------------------------------------------------------
void _export mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
#else
void mexFunction(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[0]) * mxGetM(prhs[0]);
xData = (double *) mxGetData (prhs[0]);
pTempC =(double *) mxGetPi (prhs[0]);
nArrayID =mxGetClassID(prhs[0]);
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[0] = mxCreateDoubleMatrix(j , 1, mxCOMPLEX);
pTempR = (double *)mxGetPr(plhs[0]);
pTempC = (double *)mxGetPi(plhs[0]);
pTemp= xDataC;
if(bReal){
pTempRp = pTempR + j- 1;
pTempCp = pTempC + j- 1;
*pTempR = *pTemp++;
pTempR[j>>1]=*pTemp++;
*pTempC = pTempC[j>>1] = 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算法以后会陆续给出）

recbio · 发表于 2007-8-2 17:03

假如用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 [ebp+12]
add eax,eax
dec eax
mov dword ptr [ebp-4],eax
mov edx,dword ptr [ebp+8]
add edx,8
mov dword ptr [ebp-8],edx
mov edx,ecx
cmp ecx,dword ptr [ebp-4]
jge short my@16
// ; EDX = j, ECX = i, EBX = @temp3
my@15:
cmp ecx,edx
jge short my@17
mov esi, dword ptr [ebp+8] //pData
mov eax, dword ptr [esi + 8*edx] // eax = pData[j]
mov ebx, dword ptr [esi + 8*ecx] // ebx = pData[i]
mov dword ptr [esi + 8*ecx], eax // pData[i] = eax
mov dword ptr [esi + 8*edx], ebx // pData[j] = ebx
add esi, 4
mov eax, dword ptr [esi + 8*edx] // eax = pData[j]
mov ebx, dword ptr [esi + 8*ecx] // ebx = pData[i]
mov dword ptr [esi + 8*ecx], eax // pData[i] = eax
mov dword ptr [esi + 8*edx], ebx // pData[j] = ebx
mov esi, dword ptr [ebp-8] //pDataI
mov eax, dword ptr [esi + 8*edx] // eax = pData[j]
mov ebx, dword ptr [esi + 8*ecx] // ebx = pData[i]
mov dword ptr [esi + 8*ecx], eax // pData[i] = eax
mov dword ptr [esi + 8*edx], ebx // pData[j] = ebx
add esi, 4
mov eax, dword ptr [esi + 8*edx] // eax = pData[j]
mov ebx, dword ptr [esi + 8*ecx] // ebx = pData[i]
mov dword ptr [esi + 8*ecx], eax // pData[i] = eax
mov dword ptr [esi + 8*edx], ebx // pData[j] = ebx
my@17:
mov eax,dword ptr [ebp+12]
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 [ebp-12],16
// add ebx,16
cmp ecx,dword ptr [ebp-4]
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 [ebp+12]
add eax,eax
dec eax
mov dword ptr [ebp-4],eax
mov edx,dword ptr [ebp+8]
add edx,8
mov dword ptr [ebp-8],edx
mov edx,ecx
cmp ecx,dword ptr [ebp-4]
jge short my@216
// ; EDX = j, ECX = i, EBX = @temp3
my@215:
cmp ecx,edx
jge short my@217
mov esi, dword ptr [ebp+8] //pData
fwait
fld [esi + 8*edx] // st0 = pData[j]
fld [esi + 8*ecx] // st0 = pData[i]; st1 = pData[j]
fstp [esi + 8*edx] // st0 -> pData[j] ; st0 = old pData[j]
fstp [esi + 8*ecx]
mov esi, dword ptr [ebp-8] //pDataI
fld [esi + 8*edx] // st0 = pData[j]
fld [esi + 8*ecx] // st0 = pData[i]; st1 = pData[j]
fstp [esi + 8*edx] // st0 -> pData[j] ; st0 = old pData[j]
fstp [esi + 8*ecx]
my@217:
mov eax,dword ptr [ebp+12]
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 [ebp-12],16
// add ebx,16
cmp ecx,dword ptr [ebp-4]
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 编辑 ]

纤娴毅 · 发表于 2010-8-23 20:11

办证: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
在办理英语证行业有着多年的历史。英语证是我司的主营业务之一。我们将竭诚为您服务

账号		自动登录	找回密码
密码			注册

萍聚头条

[电子] Matlab 和 C++ （4），FFT（1）

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分

相关帖子