Ты спрашивал про FFT

From

Sergey Gotsuljak (2:6055/9.24)

To

Andrey Lipsky

Date

2002-05-09T13:49:07Z

Area

RU.ALGORITHMS

AL> Ищy класснyю пpоцедypкy FFT для pеальных чисел. То есть на входе должен
AL> быть массив pеальных чисел, а на выходе - комплексных. Заpанее спасибо

Хочу предложить тебе маленькую библиотечку - сам использовал. Она работает как на реальных, так и на комплексных числах. Если знаешь язык C++ - разберешься за пару минут. Работает - проверено.


#define PI 3.14159265358979

/*
   Структура для описания сэмплов
*/

class TSample
{
   public:
      char Format[8];
      WORD MajorVersion;
      WORD MinorVersion;
      DWORD SamplingFrequency;
      WORD SourceBits;
      WORD StoreBits;
      WORD Channels;
      DWORD Number;
      float *Re;
      float *Im;
      char *Comment;
      TSample();
      TSample(int num);
      ~TSample();
      void Init();
};

TSample::TSample()
        :Re(NULL),
         Im(NULL),
         Comment(NULL)
{
   Format[0]='\0';
}

TSample::TSample(int num)
        :Re(NULL),
         Im(NULL),
         Comment(NULL)
{
   Number=num;
   Format[0]='\0';
   Re=new float[num];
   Im=new float[num];
}

TSample::~TSample()
{
   delete Format;
   delete Re;
   delete Im;
   delete Comment;
}

void TSample::Init()
{
   if(Re)
      for(DWORD i=0; i<Number; i++)
         Re[i]=0;
   if(Im)
      for(DWORD i=0; i<Number; i++)
         Im[i]=0;
}

/*
   Быстрое преобразование Фурье
*/

void FFT(TSample &in)
{
   int I=1;
   int N=in.Number;
   float c,s,t1,t2,t3,t4,u1,u2,u3;
   int i,j,p,l,L,M,M1,K;
   L=N;
   M=N/2;
   M1=N-1;
   while(L>=2)
   {
      l=L/2;
      u1=1.;
      u2=0.;
      t1=PI/(float)l;
      c=cos(t1);
      s=(-1)*I*sin(t1);
      for(j=0; j<l; j++)
      {
         for(i=j; i<N; i+=L)
         {
            p=i+l;
            t1=*(in.Re+i)+*(in.Re+p);
            t2=*(in.Im+i)+*(in.Im+p);
            t3=*(in.Re+i)-*(in.Re+p);
            t4=*(in.Im+i)-*(in.Im+p);
            *(in.Re+p)=t3*u1-t4*u2;
            *(in.Im+p)=t4*u1+t3*u2;
            *(in.Re+i)=t1;    *(in.Im+i)=t2;
         }
         u3=u1*c-u2*s;
         u2=u2*c+u1*s;  u1=u3;
      }
      L/=2;
   }
   j=0;
   for(i=0; i<M1; i++)
   {
      if(i>j)
      {
         t1=*(in.Re+j); t2=*(in.Im+j);
         *(in.Re+j)=*(in.Re+i); *(in.Im+j)=*(in.Im+i);
         *(in.Re+i)=t1; *(in.Im+i)=t2;
      }
      K=M;
      while(j>=K)
      {
         j-=K;K/=2;
      }
      j+=K;
   }
}

/*
   Обратное быстрое преобразование Фурье
*/

void IFFT(TSample &in)
{
   int I=-1;
   int N=in.Number;
   float c,s,t1,t2,t3,t4,u1,u2,u3;
   int i,j,p,l,L,M,M1,K;
   L=N;
   M=N/2;
   M1=N-1;
   while(L>=2)
   {
      l=L/2;
      u1=1.;
      u2=0.;
      t1=PI/(float)l;
      c=cos(t1);
      s=(-1)*I*sin(t1);
      for(j=0; j<l; j++)
      {
         for(i=j; i<N; i+=L)
         {
            p=i+l;
            t1=*(in.Re+i)+*(in.Re+p);
            t2=*(in.Im+i)+*(in.Im+p);
            t3=*(in.Re+i)-*(in.Re+p);
            t4=*(in.Im+i)-*(in.Im+p);
            *(in.Re+p)=t3*u1-t4*u2;
            *(in.Im+p)=t4*u1+t3*u2;
            *(in.Re+i)=t1;    *(in.Im+i)=t2;
         }
         u3=u1*c-u2*s;
         u2=u2*c+u1*s;  u1=u3;
      }
      L/=2;
   }
   j=0;
   for(i=0; i<M1; i++)
   {
      if(i>j)
      {
         t1=*(in.Re+j); t2=*(in.Im+j);
         *(in.Re+j)=*(in.Re+i); *(in.Im+j)=*(in.Im+i);
         *(in.Re+i)=t1; *(in.Im+i)=t2;
      }
      K=M;
      while(j>=K)
      {
         j-=K;K/=2;
      }                     
      j+=K;
   }
   for(DWORD i=0; i<in.Number; i++)
   {
      *(in.Re+i)=*(in.Re+i)/in.Number;
      *(in.Im+i)=*(in.Im+i)/in.Number;
   }
}

Всего хорошего!
--- FIPS/32 v0.99b W95/NT [M]
 * Origin: Жизнь не должна превращаться в рутину (2:6055/9.24)