FFT és konvolúció

F.27.

A gyors Fourier transzformációt és a konvolúciót végző C függvények DLL-je

#include <windows.h>

#include <alloc.h>

#include <math.h>

#include <stdio.h>

#define PI 3.14159265359

typedef struct { double re; // valós rész

double im; // képzetes rész

} complex;   // komplex struktúra definiálása

typedef struct { unsigned long int n;    // a függvény elemek száma

unsigned long int n0;   // az első valódi elem indexe

double *fv;             // mutató a függvény elemekre

unsigned long int ndft; // a DFT alappontok száma

complex *dft;           // mutató a DFT értékekre

} fvtype;

// struktúra, amely tartalmazza a függvényt és annak transzformáltját

static complex *exp_tabl;      // exponenciális függvény táblázata

static complex *fftbuff;       // munkamező

static unsigned long int logn; // elemszám logaritmusa

static unsigned long int n=0;  // elemszám

// Komplex összeadás

void cx_add(complex *op1,complex *op2,complex *res)

{

res->re=op1->re+op2->re;

res->im=op1->im+op2->im;

return;

}

// Komplex kivonás

void cx_sub(complex *op1,complex *op2,complex *res)

{

res->re=op1->re-op2->re;

res->im=op1->im-op2->im;

return;

}

// Komplex szorzás

void cx_mul(complex *op1,complex *op2,complex *res)

{ double xres_re;

xres_re=op1->re*op2->re-op1->im*op2->im;

res->im=op1->re*op2->im+op1->im*op2->re;

res->re=xres_re;

return;

}

// Komplex szám másolása

void cx_copy(complex *op1,complex *res)

{

res->re=op1->re;

res->im=op1->im;

return;

}

// Komplex szám reciproka

void cx_rec(complex *op1,complex *res)

{ double s;

s=op1->re*op1->re+op1->im*op1->im;

if(s<1e-20) s=1e-20;

res->re=op1->re/s;

res->im=-op1->im/s;

return;

}

// Komplex osztás

void cx_div(complex *op1,complex *op2,complex *res)

{ complex s;

cx_rec(op2,&s);

cx_mul(op1,&s,res);

return;

}

// Komplex szám abszolút értéke

double cx_abs(complex *op1)

{ double s;

s=sqrt(op1->re*op1->re+op1->im*op1->im);

return(s);

}

// Komplex exponenciális függvény

void cx_exp(complex *op1,complex *res)

{ double s;

s=exp(op1->re);

res->re=s*cos(op1->im);

res->im=s*sin(op1->im);

return;

}

// Komplex szám feltöltése valós és képzetes résszel

void cx_fill(double real,double imag,complex *res)

{

res->re=real;

res->im=imag;

return;

}

// A függvény előállítja és tárolja az exp(-jx) táblázatot

// amit az n alappontú komplex FFT igényel.

// Az alappontok száma és az exp tábla kezdőcíme is sztatikus,

// globális változó. Az exp(-jx) táblának n/2-1 eleme van,

// x=2p/n-től x=(n/2-1)2p/n-ig.

// Ezt a függvényt közvetlenül nem kell meghívni,

// a refft_prepare gondoskodik róla.

void exp_tbl()

{ complex *p_exp;

unsigned long int i;

double s,s1;

p_exp=exp_tabl;

s=2.0*PI/(double)n;

s1=s;

for(i=1;i<n/4;i++)

{ p_exp->re=cos(s1);

p_exp->im=-sin(s1);

p_exp++;

s1+=s;

}

return;

}

// A függvény elvégzi a transzformációhoz szükséges előkészítő lépéseket

// feltölti a globális változókat, munkaterületet allokál,

// kiszámolja a szinusz és koszinusz táblázat értékeit.

// Paraméterként a transzformáció alappontjának kettes alapú

// logaritmusát kapja meg.

// Visszatérési értéke 0, ha minden rendben,

// 1, ha allokációs hiba lépett fel.

extern “C” __declspec(dllexport)

int WINAPI refft_prepare(unsigned long int log_n)

{

logn=log_n;

n=1<<logn;

fftbuff=(complex *)farcalloc(n/2,sizeof(complex));

if(fftbuff==NULL)

{ printf(“\nREFFT unable to alloc/1”);

return(1);

}

exp_tabl=(complex *)farcalloc(n/4-1,sizeof(complex));

if(exp_tabl==NULL)

{ printf(“\nREFFT unable to alloc/2”);

return(1);

}

exp_tbl();

return(0);

}

// A függvény felszabadítja a transzformáció azon munkaterületeit,

// amelyeket az előkészítés során nyitottunk meg

// refft_prepare függvénnyel.

extern “C” __declspec(dllexport)

void WINAPI refft_close()

{

farfree(exp_tabl);

farfree(fftbuff);

n=logn=0;

return;

}

// Egy egydimenziós, valós érétkű adattömb DFT transzformáltjának

// számítása az FFT algoritmus szerint.

// Az eredmény a bemeneti tömbbe íródik.

// A bemeneti és az eredmény tömb is természetes sorrendű.

// Az alappontok száma static int n, amely

// refft_prepare függvénnyel állítható be.

// Paraméterként egy double értékre mutató pointert vár,

// amelyet a transzformálandó vektor kezdetére kell irányítani.

// Visszatéréskor ugyanezen a helyen lesz az eredmény is.

extern “C” __declspec(dllexport)

void WINAPI refft(double *poi)

{ 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;

p=n/2;

poi1=poi;

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++;

}

p_from=fftbuff;

p_to=(complex *)poi;

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

{ p2=p;

p>>=1;

nper2p=n/p2;

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;

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;

}

}

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)

{ poi1=(double *)fftbuff;

poi2=poi;

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

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

}

s=1.0/(double)n;

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

return;

}

// Inverz DFT transzformáció FFT algoritmus alkalmazásával,

// feltételezve, hogy az eredmény elemei valósak.

// Hasonlóan működik refft-hez.

extern”C” __declspec(dllexport)

void WINAPI reifft(double *poi)

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

double *poi1,*poi2;

complex *p_exp;

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

complex cx,cxexp;

p_from=(complex *)poi;

p_to=fftbuff;

p=1;

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

{ p2=2*p;

nper2p=n/p2;

p2_to=p+(p1_to=p_to);

p1_from=p_from;

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

p2_to->re=p1_from->re-p1_from->im;

p1_to++;

p2_to++;

p1_from++;

}

p1_to+=p;

p2_to+=p;

p2_from=p_from+n/2-p;

if(n>=4*p)

{ p_exp=exp_tabl+p-1;

nper4p=nper2p/2;

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

{ cx_copy(p_exp,&cxexp);

cxexp.im=-cxexp.im;

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

{ p2_from->im=-p2_from->im;

cx_add(p1_from,p2_from,p1_to);

cx_sub(p1_from,p2_from,&cx);

cx_mul(&cx,&cxexp,p2_to);

p1_to++;

p2_to++;

p1_from++;

p2_from++;

}

p_exp+=p;

p1_to+=p;

p2_to+=p;

p2_from-=p2;

}

}

p2_to=p+(p1_to=p_to);

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

{ p1_to->im=2.0*p1_from->re;

p2_to->im=-2.0*p1_from->im;

p1_to++;

p2_to++;

p1_from++;

}

p1_to=p_to;

p_to=p_from;

p_from=p1_to;

p=p2;

}

p=n/2;

poi1=(double *)p_to;

poi2=poi1+p;

p1_from=p_from;

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

{ *poi1=p1_from->re+p1_from->im;

*poi2=p1_from->re-p1_from->im;

poi1++;

poi2++;

p1_from++;

}

if(logn%2)

{ poi1=(double *)fftbuff;

poi2=poi;

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

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

}

return;

}

// A függvény az üres állapotnak megfelelően tölti ki a fvtype struktúrát.

// Paraméterként egy fvtype struktúrára mutató pointert vár.

// Bármely fvtype struktúrával végzendő művelet előtt be kell hívni.

extern “C” __declspec(dllexport)

void WINAPI init_fvtype(fvtype *fvt)

{

fvt->n=0;

fvt->n0=0;

fvt->fv=NULL;

fvt->ndft=0;

fvt->dft=NULL;

return;

}

// Felszabadítja az operandus fvtype mutatói által mutatott tárterületeket.

void free_fvtype(fvtype *fvt)

{

if(fvt->fv!=NULL) farfree(fvt->fv);

if(fvt->dft!=NULL) farfree(fvt->dft);

fvt->fv=NULL;

fvt->dft=NULL;

return;

}

// Egy adatsornak megfelelő struktúra kitöltése.

// A szükséges hely allokálása és az adatsor átmásolása.

// Paraméterek:

// n: az adatsor adatainak száma,

// n0:index offset, az első adat valódi indexe,

// data: az fvtype struktúrába másolandó adatsor mutatója,

// fvt: a kitöltendő struktúra mutatója.

// Sikertelen allokáció esetén a függvény hibát jelez.

extern “C” __declspec(dllexport)

int WINAPI fill_fvtype(unsigned long int n,unsigned long int n0,double *data,fvtype *fvt)

{ unsigned long int i;

double *dp;

free_fvtype(fvt);

fvt->fv=(double*)farcalloc(n,sizeof(double));

if(fvt->fv==NULL) return(1);

dp=fvt->fv;

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

*(dp++)=*(data++);

fvt->n=n;

fvt->n0=n0;

return(0);

}

// Az fvtype struktúrában álló függvény DFT transzformálása

// ndft alappontra, ha még nem történt meg.

int dft_fvtype(fvtype *fvt,unsigned long int ndft)

{ unsigned long int i;

double *dp1,*dp2;

if(fvt->dft!=NULL && fvt->ndft==ndft) return(0);

if(fvt->fv==NULL || fvt->n==0) return(1);

if(fvt->n>ndft) return(3);

if(fvt->dft!=NULL) farfree(fvt->dft);

fvt->dft=(complex*)farcalloc(ndft,sizeof(double));

if(fvt->dft==NULL) return(2);

dp1=(double *)fvt->dft;

dp2=fvt->fv;

for(i=0;i<fvt->n;i++) *(dp1++)=*(dp2++);

for(i=fvt->n;i<ndft;i++) *(dp1++)=0.0;

refft((double *)fvt->dft);

fvt->ndft=ndft;

return(0);

}

// Két függvény konvolúcióját megvalósító rutin a DFT transzformáltak

// segítségével.

// Sikertelen allokáció esetén eredmény nélkül tér vissza a függvény.

extern “C” __declspec(dllexport)

int WINAPI convol_by_fft(fvtype *fvta,fvtype *fvtb,fvtype *fvtres)

{ long int i,na,nb,nres,ndft,logndft;

double s,*dbuff,*dp1,*dp2;

complex *cp1,*cp2,*cp3;

na=fvta->n;

nb=fvtb->n;

nres=na+nb-1;

ndft=1;

logndft=0;

while(ndft<nres) { ndft<<=1; logndft+=1; }

refft_prepare(logndft);

if(dft_fvtype(fvta,ndft)) goto alloc_err;

if(dft_fvtype(fvtb,ndft)) goto alloc_err;

free_fvtype(fvtres);

fvtres->dft=(complex*)farcalloc(ndft,sizeof(double));

if(fvtres->dft==NULL) goto alloc_err;

cp1=fvta->dft;

cp2=fvtb->dft;

cp3=fvtres->dft;

cp3->re=cp1->re*cp2->re;

cp3->im=cp1->im*cp2->im;

cp1++; cp2++; cp3++;

for(i=1;i<ndft/2;i++) cx_mul(cp1++,cp2++,cp3++);

fvtres->ndft=ndft;

dbuff=(double*)farcalloc(ndft,sizeof(double));

if(dbuff==NULL) goto alloc_err;

dp1=(double *)fvtres->dft;

dp2=dbuff;

s=(double)ndft;

for(i=0;i<ndft;i++) *(dp2++)=*(dp1++)*s;

reifft(dbuff);

fvtres->fv=(double*)farcalloc(nres,sizeof(double));

if(fvtres->fv==NULL) goto alloc_err2;

dp1=dbuff;

dp2=fvtres->fv;

for(i=0;i<nres;i++) *(dp2++)=*(dp1++);

fvtres->n=nres;

fvtres->n0=fvta->n0+fvtb->n0;

farfree(dbuff);

refft_close();

return(0);

alloc_err2:

farfree(dbuff);

alloc_err:

refft_close();

return(1);

}

Vélemény, hozzászólás?

Adatok megadása vagy bejelentkezés valamelyik ikonnal:

WordPress.com Logo

Hozzászólhat a WordPress.com felhasználói fiók használatával. Kilépés / Módosítás )

Twitter kép

Hozzászólhat a Twitter felhasználói fiók használatával. Kilépés / Módosítás )

Facebook kép

Hozzászólhat a Facebook felhasználói fiók használatával. Kilépés / Módosítás )

Google+ kép

Hozzászólhat a Google+ felhasználói fiók használatával. Kilépés / Módosítás )

Kapcsolódás: %s