Hangkártyák programozása retro 19.

7.1.8. FFT valós vektor esetén

A (3.34) szerinti általános rész-transzformáltra áll, hogy

(7.40)

Ebből az f_i valós esetben egyenesen következik, hogy

(7.41)

Ahol * a konjugált képzést jelenti. Más szóval:

(7.42)

számsorozat redundáns, hiszen az r=N/2p fölötti rész konjugálással képezhető a kisebb indexűekből.

Nézzük, mennyi a transzformált nem redundáns része. Ez a következőkből áll:

valós érték (7.34) szerint egyértelmű,

valós érték,

komplex értékek.

Ha összeszámláljuk, kiderül, hogy ez éppen N darab adatot jelent (a komplex számot itt két adatnak vettük). Ez megfelel a várakozásnak, hiszen a transzformáció N darab adatból indul.

A valós bemenő vektorra egyszerűsített algoritmusnál tehát nem tároljuk a teljes rész-transzformáltat, csak a nem redundáns részt. Ennek megfelelően át kell rendeznünk a (7.39) szerinti rekurziós formulát. Ez a rendezés az alábbi eredményt adja, amelyet nem részletezek:

                                        (7.43)

 

(7.44)

 

(7.45)

Egy N/2 elemű komplex vektorban a 0 indexű elem nem tényleges komplex szám, hanem a valós illetőleg a képzetes rész helyén a transzformált két tisztán valós elemét, D₀ és D_N/2 –t tároljuk. A transzformáció eredményét ezen rend szerint kapjuk. Ez a tárolási séma a rész-transzformáltakra is vonatkozik.

A transzformációt végző C függvény [6] alapján a következő:

//A megfelelő paraméterátadási konvenció beállítása: WINAPI,

//és a függvény láthatóvá tétele a külvilág számára: __declspec(dllexport)

extern “C” __declspec(dllexport)

void WINAPI refft(double *poi) //transzformálandó vektor, mint paraméter

{ unsigned long int p,p2,q,r,i,nper2p,nper4p;

double s,*poi1,*poi2;

complex *p_exp;

complex *p_from,*p_to,*p1_from,*p1_to,*p2_from,*p2_to;

complex cx,cxexp;

//Első rekurziós lépés megvalósítása

p=n/2;

poi1=poi; //Munkaterületek kezdőcímeinek kijelölése

poi2=poi+p;

p1_to=fftbuff;

for(q=0;q<p;q++)

{ p1_to->re=*poi1+(*poi2);

p1_to->im=*poi1-(*poi2);

poi1++;

poi2++;

p1_to++;

}

//A további rekurziós lépések ciklusa

p_from=fftbuff;

p_to=(complex *)poi;

for(i=1;i<logn;i++)

{ p2=p;

p>>=1;

nper2p=n/p2;

//Az r=0 és r=n/2p eset

p2_from=p+(p1_from=p_from);

p1_to=p_to;

for(q=0;q<p;q++) { p1_to->re=p1_from->re+p2_from->re;

p1_to->im=p1_from->re-p2_from->re;

p1_from++;

p2_from++;

p1_to++;

}

p1_from+=p;

p2_from+=p;

p2_to=p_to+n/2-p;

//Általános eset

if(n>=4*p)

{ p_exp=exp_tabl+p-1;

nper4p=nper2p/2;

for(r=1;r<nper4p;r++)

{ cx_copy(p_exp,&cxexp);

for(q=0;q<p;q++)

{ cx_mul(p2_from,&cxexp,&cx);

cx_add(p1_from,&cx,p1_to);

cx_sub(p1_from,&cx,p2_to);

p2_to->im=-p2_to->im;

p1_from++;

p2_from++;

p1_to++;

p2_to++;

}

p_exp+=p;

p1_from+=p;

p2_from+=p;

p2_to-=p2;

}

}

//Az r=n/4p eset

p2_from=p+(p1_from=p_from);

for(q=0;q<p;q++)

{ p1_to->re=p1_from->im;

p1_to->im=-p2_from->im;

p1_from++;

p2_from++;

p1_to++;

}

p1_from=p_from;

p_from=p_to;

p_to=p1_from;

}

if(logn%2) //Visszamásolás a poi által mutatott adaterületre

{ poi1=(double *)fftbuff;

poi2=poi;

for(i=0;i<n;i++)

*(poi2++)=*(poi1++);

}

s=1.0/(double)n; //Az eredmény végigosztása n-nel

for(i=0;i<n;i++) *(poi++)*=s;

return;

}

Az inverz transzformáció hasonló meggondolások szerint kódolható, listája megtalálható a [6] alapján készült F.27-ben.

A fenti függvényt Delphi-ben a következőképpen kell deklarálni:

procedure ReFft(poi:pdouble); //pdouble=^double

stdcall; external ‘prjconvoldll.dll’ name ‘refft’;

ahol prjconvoldll.dll az állomány, illetve a dinamikusan szerkesztett könyvtár, ahonnét a kérdéses függvényt behívjuk, refft pedig a C függvény eredeti neve.

 

[6]  Székely Vladimír: Képkorrekció, hanganalízis, térszámítás PC-n (Gyors Fourier transzformációs módszerek), Computer Books, Budapest, 1994