Отправлено ВН 06 февраля 2003 г. 21:41 В ответ на: студенческая работа периода учебы отправлено net 06 февраля 2003 г. 17:46

Знакомый до боли текст. Кули-Тьюки.
Вот он же, тупо переписанный на C. Ну чуть-чуть измененный.
inv=0 - прямое, =1 - обратное.
len - длина, step - число ступеней, логарифм двоичный длины.
re - реальная, im-мнимая части вх.-вых. данных.
#include
void fft(double far *re,double far *im,int len,int step,int inv)
{
double PI=3.14159265358979323846;
double ur,ui,tr,ti,g,e;
int i,j,k,l,n,nv2,nm1,le,le1;
n=len;
nv2=n/2;
nm1=n-1;
j=1;
for(i=1;i<=nm1;i++)
{
if(i{
tr=re[j-1];
ti=im[j-1];
re[j-1]=re[i-1];
im[j-1]=im[i-1];
re[i-1]=tr;
im[i-1]=ti;
}
k=nv2;
while(k{
j-=k;
k/=2;
}
j+=k;
}
le=1;
for(l=1;l<=step;l++)
{
le1=le;
le*=2;
g=PI/((double)le1);
for(j=1;j<=le1;j++)
{
e=g*((double)(j-1));
ur=cos(e);
ui=sin(e);
if(!inv) ui=-ui;
for(i=j;i<=n;i+=le)
{
k=i+le1;
tr=re[k-1]*ur-im[k-1]*ui;
ti=re[k-1]*ui+im[k-1]*ur;
re[k-1]=re[i-1]-tr;
im[k-1]=im[i-1]-ti;
re[i-1]+=tr;
im[i-1]+=ti;
}
}
}
if(inv)
{
for(i=0;i{
re[i]/=n;
im[i]/=n;
}
}
}
А вот еще вариант. На 256 точек. Только в одну сторону, прямое т.е.
Обратное - знаки кое-где в бабочке поменять и результат на длину поделить.
float DAT[512];//DAT[2*i] -real, DAT[2*i+1] - imidg;
float COSARR[256];//COSARR[i]=cos(2*pi*i/256);
float SINARR[256];//SINARR[i]=sin(2*pi*i/256);
int bitrev128(int i)
{
int j=((i&1)<<6) | ((i&2)<<4) | ((i&4)<<2) | | ((i&16)>>2) | ((i&32)>>4) |((i&64)>>6);
return j;
}
int bitrev256(int i)
{
int j=((i&1)<<7) | ((i&2)<<5) | ((i&4)<<3) | ((i&8)<<1) | ((i&16)>>1) | ((i&32)>>3) | ((i&64)>>5) |((i&128)>>7);
return j;
}
void fft256(void)
{
float RA,IA,RB,IB,RW,IW;
float RE,IM;
float LI;
int NSTAG,NSUBSTAG,NBATTER;
int n,m,l,k,j,i,mm,ll;
NSTAG=8;
NSUBSTAG=1;
NBATTER=128;
for(n=0;n{
mm=2*NBATTER;
for(m=0;m{
k=bitrev128(m);
RW=COSARR[k];
IW=SINARR[k];
for(l=0;l{
k=4*m*NBATTER+2*l;
ll=k+mm;
RA=DAT[k];
IA=DAT[k+1];
RB=DAT[ll];
IB=DAT[ll+1];
RE=RB*RW+IB*IW;
IM=IB*RW-RB*IW;
DAT[k]=RA+RE;
DAT[k+1]=IA+IM;
DAT[ll]=RA-RE;
DAT[ll+1]=IA-IM;
}
}
NBATTER>>=1;
NSUBSTAG<<=1;
}
for(i=0;i<256;i++)
{
j=bitrev256(i);
if(i>=j) continue;
RE=DAT[2*i];
IM=DAT[2*i+1];
DAT[2*i]=DAT[2*j];
DAT[2*i+1]=DAT[2*j+1];
DAT[2*j]=RE;
DAT[2*j+1]=IM;
}
}
И последний, для целых чисел. С блочной плавающей запятой, она и возвращается. Тоже 256. DAT,COSARR,SINARR - теперь уже int, bitrev128,bitrev256 -те же. Ф-ии shiftcalcul,void nrmlzarr - для бл. ПЗ.
int shiftcalcul(void)
{
int i,j,k,l;
j=0;
for(i=0;i<512;i++)
{
k=DAT[i];
if(k<0)
{
if(k==0x8000) k=32767;
else k=-k;
}
if(k>j) j=k;
}
if(!j) return 13;
k=0x4000;
for(i=0;i<15;i++)
{
if((j&k)) {l=i;break;}
k>>=1;
}
return (l-2);
}
void nrmlzarr(int shft)
{
int i,j;
if(!shft) return;
if(shft<0)
{
j=-shft;
for(i=0;i<512;i++) DAT[i]>>=j;
}
else for(i=0;i<512;i++) DAT[i]<<=shft;
}
int fft256(void)
{
long RB,IB,RW,IW,LI;
int RE,IM,RA,IA;
int EXP,NSTAG,NSUBSTAG,NBATTER;
int n,m,l,k,j,i,ll,mm;
NSTAG=8;
NSUBSTAG=1;
NBATTER=128;
EXP=0;
for(n=0;n{
j=shiftcalcul();
EXP-=j;
nrmlzarr(j);
mm=2*NBATTER;
for(m=0;m{
k=bitrev128(m);
RW=(long)COSARR[k];
IW=(long)SINARR[k];
for(l=0;l{
k=4*m*NBATTER+2*l;
ll=k+mm;
RA=DAT[k];
IA=DAT[k+1];
RB=(long)DAT[ll];
IB=(long)DAT[ll+1];
LI=(RB*RW+IB*IW)<<1;
RE=(int)(LI>>16);
LI=(IB*RW-RB*IW)<<1;
IM=(int)(LI>>16);
DAT[k]=RA+RE;
DAT[k+1]=IA+IM;
DAT[ll]=RA-RE;
DAT[ll+1]=IA-IM;
}
}
NBATTER>>=1;
NSUBSTAG<<=1;
}
for(i=0;i<256;i++)
{
j=bitrev256(i);
if(i>=j) continue;
RE=DAT[2*i];
IM=DAT[2*i+1];
DAT[2*i]=DAT[2*j];
DAT[2*i+1]=DAT[2*j+1];
DAT[2*j]=RE;
DAT[2*j+1]=IM;
}
return EXP;
}

• Ответ: Да... — ВН   (06.02.2003 21:46, 84 байт)
  • Надо перед текстом (в отдельной строчке) в угловых скобках написать pre, а после - /pre — Xan   (07.02.2003 12:31, 54 байт)
    • Ответ: Спасибо, сейчас попробую. — ВН   (07.02.2003 16:01, 1340 байт)
      • Ответ: Лучше, но.... Еще раз. — ВН   (07.02.2003 16:15, 1503 байт)
        • Ответ: Еще — ВН   (07.02.2003 16:23, 1265 байт)
          • Ответ: И еще — ВН   (07.02.2003 16:29, 1375 байт)
            • Ответ: Ну хоть так. — ВН   (07.02.2003 16:31, пустое)