void rfft(int n, double x[]) //實序列快速傅里葉變換 
{ 
int i,j,k,m,i1,i2,i3,i4,n1,n2,n4; 
double a,e,cc,ss,xt,t1,t2; 
for(j=1,i=1;i <16;i++) 
{ 
m=i; 
j= 2*j; 
if(j==n) break; 
} 
n1=n-1; 
for(j=0,i=0;i <n1;i++) 
{ 
if(i <j) 
{ 
xt=x[j]; 
x[j]=x[i]; 
x[i]=xt; 
} 
k=n/2; 
while (k <(j+1)) { 
j=j-k; 
k=k/2; 
} 
j=j+k; 
} 
for(i=0;i <n;i+=2) 
{ 
xt=x[i]; 
x[i]=xt+x[i+1]; 
x[i+1]=xt-x[i+1]; 
} 
n2=1; 
for(k=2;k <=m;k++) 
{ 
n4=n2; 
n2=2*n4; 
n1=2*n2; 
e=6.28318530718/n1; 
for(i=0;i <n;i+=n1) 
{ 
xt=x[i]; 
x[i]=xt+x[i+n2]; 
x[i+n2]=xt-x[i+n2]; 
x[i+n2+n4]=-x[i+n2+n4];// 
a=e; 
for(j=1;j <(n4-1);j++) 
{ 
i1=i+j; 
i2=i-j+n2; 
i3=i+j+n2; 
i4=i-j+n1; 
cc=cos(a); 
ss=sin(a); 
a=a+e; 
t1=cc*x[i3]+ss*x[i4]; 
t2=ss*x[i3]-cc*x[i4]; 
x[i4]=x[i2]-t2; 
x[i3]=-x[i2]-t2;// 
x[i2]=x[i1]-t1; 
x[i1]=x[i1]+t1; 
} 
} 
} 
}