По мотивам знаменитой Numerical Recipes
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr void four1(float data[], unsigned long nn); void realft(float data[], unsigned long n); void realft(float data[], unsigned long n) { unsigned long i,i1,i2,i3,i4,np3; float c1=0.5,c2,h1r,h1i,h2r,h2i; double wr,wi,wpr,wpi,wtemp,theta; theta=3.141592653589793/(double) (n>>1); c2 = -0.5; four1(data,n>>1); data--; wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0+wpr; wi=wpi; np3=n+3; for (i=2;i<=(n>>2);i++) { i4=1+(i3=np3-(i2=1+(i1=i+i-1))); h1r=c1*(data[i1]+data[i3]); h1i=c1*(data[i2]-data[i4]); h2r = -c2*(data[i2]+data[i4]); h2i=c2*(data[i1]-data[i3]); data[i1]=float(h1r+wr*h2r-wi*h2i); data[i2]=float(h1i+wr*h2i+wi*h2r); data[i3]=float(h1r-wr*h2r+wi*h2i); data[i4] = float(-h1i+wr*h2i+wi*h2r); wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; } data[1] = (h1r=data[1])+data[2]; data[2] = h1r-data[2]; } void four1(float data[], unsigned long nn) { unsigned long n,mmax,m,j,istep,i; double wtemp,wr,wpr,wpi,wi,theta; // двойная точность для тригонометрических вычислений float tempr,tempi; data--; n=nn << 1; j=1; for (i=1;i<n;i+=2) { // рекомбинация элементов в бит-реверсивном порядке if (j > i) { SWAP(data[j],data[i]); 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=mmax << 1; theta=6.28318530717959/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=float(wr*data[j]-wi*data[j+1]); tempi=float(wr*data[j+1]+wi*data[j]); data[j]=data[i]-tempr; data[j+1]=data[i+1]-tempi; data[i] += tempr; data[i+1] += tempi; } wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; } mmax=istep; } } #undef SWAP