|
回复 #1 t12211017 的帖子
- #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
- fft( double* data, int nn, int isign)
- // numerical recipies routine "four1.c"
- {
- int n,mmax,m,j,istep,i;
- double wtemp,wr,wpr,wpi,wi,theta;
- double tempr,tempi;
- n=nn << 1;
- j=1;
- for (i=1;i<n;i+=2) {
- if (j > i) {
- SWAP(data[j],data);
- SWAP(data[j+1],data[i+1]);
- }
- m=n >> 1;
- while (m >= 2 && j > m) {
- j -= m;
- m >>= 1;
- }
- j += m;
- }
- mmax=2;
- while (n > mmax) {
- istep=2*mmax;
- theta=6.28318530717959/(isign*mmax);
- wtemp=sin(0.5*theta);
- wpr = -2.0*wtemp*wtemp;
- wpi=sin(theta);
- wr=1.0;
- wi=0.0;
- for (m=1;m<mmax;m+=2) {
- for (i=m;i<=n;i+=istep) {
- j=i+mmax;
- tempr=wr*data[j]-wi*data[j+1];
- tempi=wr*data[j+1]+wi*data[j];
- data[j]=data-tempr;
- data[j+1]=data[i+1]-tempi;
- data += tempr;
- data[i+1] += tempi;
- }
- wr=(wtemp=wr)*wpr-wi*wpi+wr;
- wi=wi*wpr+wtemp*wpi+wi;
- }
- mmax=istep;
- }
- }
- #undef SWAP
复制代码
[ 本帖最后由 风花雪月 于 2007-12-24 11:14 编辑 ] |
评分
-
1
查看全部评分
-
|