ru.algorithms

 
 - RU.ALGORITHMS ----------------------------------------------------------------
 From : Ilia Kantor                          2:5020/175.2   16 Dec 2002  21:27:27
 To : Nick Poroshin
 Subject : FFT - faq
 -------------------------------------------------------------------------------- 
 

 Mon Dec 16 2002 10:10, Nick Poroshin wrote to Evgeny Sharandin:
 
  NP> Можешь пpивести лучший на твой взгляд cdft/rdft в pамках 4-5 кб?
 
 Интересно, какие результаты покажет преобразование, данное ниже ?
 
 Прилагаемая реализация несколько упрощена в принципиальных местах, поэтому
 всегда будет уступать Ооуре. В реальных тестах, на больших количествах точек
 отставание где-то в 2 раза по времени и незначительное - по точности.
 
 Пример, тем не менее, учебный.. Потому - простой.
 ---------------- fft.cpp -------------------------
 
 #include <math.h>
 #include "general.h"
 #include "Complex.h"
 
 #define TRIG_VARS                  \
   ulong TLen,TNdx;int TDir;   \
   Complex PRoot,Root;
 
 #define INIT_TRIG(LENGTH,DIR)         \
   TNdx=0;TLen=(LENGTH);TDir=(DIR);     \
   PRoot.r=1.0;PRoot.i=0.0;            \
   Root.r=sin(CONST_PI/((LENGTH)*2.0));    \
   Root.r=-2.0*Root.r*Root.r;           \
   Root.i=sin(CONST_PI/(LENGTH))*(DIR);
 
 #define NEXT_TRIG_POW                     \
   if (((++TNdx)&15)==0) {                 \
    real Angle=(CONST_PI*(TNdx))/TLen;      \
    PRoot.r=sin(Angle*0.5);            \
    PRoot.r=1.0-2.0*PRoot.r*PRoot.r;   \
    PRoot.i=sin(Angle)*(TDir);         \
   } else {                     \
    Complex Temp;                       \
    Temp=PRoot;                        \
    PRoot = PRoot*Root;         \
    PRoot = PRoot+Temp;         \
   }
 inline ulong rev_next(ulong r, ulong n) { 
   do {
    n = n >> 1;
    r = r^n;
   } while ( (r&n) == 0);
   return r;
 }
 void FFTReOrder(Complex *Data, ulong Len) {
   Complex temp;
   if (Len <= 2) return;
   ulong r=0;
   for ( ulong x=1; x<Len; x++) {
    r = rev_next(r, Len);
    if (r>x) { temp=Data[x]; Data[x]=Data[r]; Data[r]=temp; } 
   }
 }
 void IFFT_T(Complex *Data, ulong Len, int Dir) {
   ulong Step, HalfStep;
   ulong b;
   TRIG_VARS;
 
   Step = 1;
   while (Step < Len) {
    HalfStep = Step;
    Step *= 2;
 
    INIT_TRIG(HalfStep,Dir);
 
    for (b = 0; b < HalfStep; b++) {
       ulong L,R;
       for (L=b; L<Len; L+=Step) {
        Complex TRight,TLeft;
        R=L+HalfStep;
        TLeft=Data[L];TRight=Data[R];
        TRight  = TRight * PRoot;
        Data[L] = TLeft + TRight;
        Data[R] = TLeft - TRight;
       }
       NEXT_TRIG_POW;
    }
 
   }
 }
 void FFT_T(Complex *Data, ulong Len, int Dir) {
   ulong k;
 
   TRIG_VARS;
 
   if (Len <= (CACHE_HALF/sizeof(Complex)) ) {
    IFFT_T(Data, Len,Dir);
    return;
   }
 
   Len /= 2;
 
   INIT_TRIG(Len, Dir);
 
   FFT_T(Data,    Len,Dir);
   FFT_T(Data+Len,Len,Dir);
 
   for (k=0; k<Len; k++) {
    Complex b,c;
    b=Data[k];
    c = Data[k+Len] * PRoot;
    Data[k] = b + c;
    Data[k+Len] = b - c;
    NEXT_TRIG_POW;
   }
 }
 void RealFFT(real *ddata, ulong Len, int Dir) {
   ulong i, j;
   Complex *Data=(Complex*)ddata;
   TRIG_VARS;
 
   Len /= 2;
 
   if (Dir > 0) {
    FFTReOrder(Data,Len);
    FFT_T(Data,Len,1);
   }
 
   INIT_TRIG(Len,Dir);
   NEXT_TRIG_POW;
 
   for (i = 1, j = Len - i; i < Len/2; i++, j--) {
        Complex p1,p2,t;
    t = conj(Data[j]);
    p1 = Data[i] + t;
    p2 = Data[i] - t;
    p2 = p2 * PRoot;
 
    t = Complex(-Dir*p2.i,Dir*p2.r);
 
    Data[i] = p1 + t;
    Data[j] = p1 - t;
    Data[j] = conj(Data[j]);
 
    Data[i] = Data[i]/2;
    Data[j] = Data[j]/2;
 
    NEXT_TRIG_POW;
   }
 
   {   real r,i;
    r=Data[0].r;i=Data[0].i;
    Data[0] = Complex(r+i,r-i);
   }
 
   if (Dir < 0) {
    Data[0] = Data[0]/2.0;
    FFTReOrder(Data,Len);
    FFT_T(Data,Len,-1);
   }
 }
 ------------------- general.h --------------
 #define CACHE_HALF 65536
 #define CONST_PI      3.1415926535897932384626433832
 typedef double real;
 typedef unsigned long ulong;
 typedef unsigned short ushort;
 ---------------- complex.h ----------------------------
 class Complex {
 public:
   real r, i; 
   Complex(void) { }
   Complex(real a, real b) { r=a; i=b; }
   inline const Complex operator+(const Complex &c) const {
    return Complex( r + c.r, i + c.i);
   }
   inline const Complex operator-(const Complex &c) const {
    return Complex( r - c.r, i - c.i);
   }   
   inline const Complex operator*(const Complex &c) const {
    return Complex( r*c.r - i*c.i, r*c.i + i*c.r);
   }
   inline const Complex operator/(const real &divisor) const {
    return Complex( r/divisor, i/divisor);
   }
 
 };
 
 inline const Complex conj(const Complex &c) {
   return Complex( c.r, -c.i);
 }
 
 --- ifmail v.2.15dev5
  * Origin: FidoNet Online - http://www.fido-online.com (2:5020/175.2)
 
 

Вернуться к списку тем, сортированных по: возрастание даты  уменьшение даты  тема  автор 

Тема:    Автор:    Дата:    FFT - faq   Nick Poroshin   10 Dec 2002 21:40:08   FFT - faq   Evgeny Sharandin   13 Dec 2002 03:06:00   FFT - faq   Ilia Kantor   14 Dec 2002 03:37:16   FFT - faq   Evgeny Sharandin   18 Dec 2002 02:27:00   FFT - faq   Nick Poroshin   16 Dec 2002 11:10:21   FFT - faq   Ilia Kantor   16 Dec 2002 21:27:27   FFT - faq   Nick Poroshin   18 Dec 2002 01:23:21   FFT - faq   Ilia Kantor   18 Dec 2002 00:59:27   FFT - faq   Nick Poroshin   18 Dec 2002 05:04:07   Re: FFT - faq   Dmitry Shuklin   18 Dec 2002 13:43:02   FFT - faq   Nick Poroshin   18 Dec 2002 23:56:18   FFT - faq   Ilia Kantor   14 Dec 2002 03:47:33   FFT - faq   Vladimir Vassilevsky   14 Dec 2002 06:55:03   FFT - faq   Ilia Kantor   14 Dec 2002 19:01:03   FFT - faq   Denis Artuhov   14 Dec 2002 17:09:55   Re: FFT - faq   Michael Ryazanov   17 Dec 2002 00:09:00   FFT - faq   Nick Poroshin   16 Dec 2002 11:23:51   FFT - faq   Ilia Kantor   16 Dec 2002 21:31:15   Re: FFT - faq   Michael Ryazanov   17 Dec 2002 00:16:00   FFT - faq   Evgeny Sharandin   18 Dec 2002 02:35:00