Отправлено -=ВН=- 29 марта 2004 г. 13:02 В ответ на: Господа, Товарищи не дайте пропасть в бездне незнания. Хочется получить спектр из выборок... отправлено barmer 29 марта 2004 г. 10:25

может поможет. На 64 точки, прореживание во времени. Реально работает, но использовать не рекомендуется - медленный. Просто для понятия. Он с бл. пл. запятой. Возвращает эту самую бл. пл. запятую (EXP). К сугубо целочисленному виду может приведен ликвидацией вызовов expcalcul,nrmlzarr, вычисления EXP.

int DAT[128];
int COSARR[64];
int SINARR[64];
//int NSTAG=6;    //NSTAG=log2(length), ÷èñëî ñòóïåíåé fft
//int NSUBSTAG=1; //÷èñëî ïîäñòóïåíåé íà îäíîé ñòóïåíè. Íà ïåðâîé -1, íà 2-îé -2, 3 -4, è ò.ä.
//int NBATTER=32; //÷èñëî áàáî÷åê â ïîäñòóïåíè. Íà îäíîé è òîé æå ñòóïåíè â êàæäîé ïîäñòóïåíè îäíî è òî æå ÷èñëî áàáî÷åê.
union lonint
{
 long l;
 int i[2];
};
int bitrev32(int i)
{
 int j=((i&1)<<4) | ((i&2)<<2) | (i&4) | ((i&8)>>2) | ((i&16)>>4);
 return j;
}
int bitrev64(int i)
{
 int j=((i&1)<<5) | ((i&2)<<3) | ((i&4)<<1) | ((i&8)>>1) | ((i&16)>>3) | ((i&32)>>5);
 return j;
}
int expcalcul(void)
{
 int i,j,k,l;
 j=0;
 for(i=0;i<128;i++)  //search of max
  {
   k=DAT[i];
   if(k<0)
    {
     if(k==0x8000) k=32767;
     else k=-k;
    }
   if(k>j) j=k;
  }
// define exponent of maximum
 if(!j) return 13;
 k=0x4000;
 for(i=0;i<15;i++)
  {
   if((j&k)) {l=i;break;}
   k>>=1;
  }
 return (l-2);     //2 çàùèòíûõ áèòà äëÿ èñêëþ÷åíèÿ ïåðåïîëíåíè
}
void nrmlzarr(int shft)
{
 int i,j;
 if(!shft) return;
 if(shft<0)
  {
   j=-shft;
  for(i=0;i<128;i++) DAT[i]>>=j;
  }
 else for(i=0;i<128;i++) DAT[i]<<=shft;
}
int fft(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=6;
 NSUBSTAG=1;
 NBATTER=32;
 EXP=0;
 for(n=0;n<NSTAG;n++)
  {
   j=expcalcul();
   EXP-=j;
   nrmlzarr(j);
   mm=2*NBATTER;
   for(m=0;m<NSUBSTAG;m++)
    {
     k=bitrev32(m);
     RW=(long)COSARR[k];
     IW=(long)SINARR[k];
     for(l=0;l<NBATTER;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<64;i++)
   {
    j=bitrev64(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;
}
static double resul[64];
void main(void)
{
 int i,j,k,ex;
 double x,y,z,fi,v,w;
 double PI=3.14159265358979323846;
 fi=2.*PI/64.;
 for(i=0;i<64;i++)
  {
   z=fi*(double)i;
   x=32767.*cos(z);
   if(x>32767.) x=32767.;
   if(x<(-32767.)) x=-32767.;
   y=32767.*sin(z);
   if(y>32767.) y=32767.;
   if(y<(-32767.)) y=-32767.;
   COSARR[i]=(int)x;
   SINARR[i]=(int)y;
  }
 for(i=0;i<128;i++) DAT[i]=0.;
 for(i=0;i<8;i++) DAT[2*i]=32767.;
 i=fft();
 v=(double)(i-6);
 w=pow(2.,v);
 for(i=0;i<64;i++)
  {
   x=(double)DAT[2*i];
   y=(double)DAT[2*i+1];
   x/=32768.;
   y/=32768.;
   x*=v;
   y*=v;
   resul[i]=sqrt(x*x+y*y);
  } 
}

• Ответ: — barmer   (30.03.2004 05:33, 1145 байт)
  • Ответ: — -=ВН=-   (30.03.2004 11:05, 945 байт)
• Ответ: — barmer   (30.03.2004 05:31, пустое)