По мотивам знаменитой 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
