|
ru.algorithms
- RU.ALGORITHMS ----------------------------------------------------------------
From : Sergey Gotsuljak 2:6055/9.24 09 May 2002 13:49:07
To : Andrey Lipsky
Subject : Ты спрашивал про FFT
--------------------------------------------------------------------------------
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)
Вернуться к списку тем, сортированных по: возрастание даты уменьшение даты тема автор
|
Ты спрашивал про FFT 
