Diskuze: IFFT pomocí FFTW

C++ C a C++ IFFT pomocí FFTW

Avatar
drob.cl
Člen
Avatar
drob.cl:

Ahoj,
měl bych další dotaz k použití FFTW, viz pokračování tohoto:
http://www.itnetwork.cz/…45d041e281e7
Snažím se napsat Matlabovskou fci hilbert
(http://www.mathworks.com/…hilbert.html)
v C++ a narazil jsem na problém při výpočtu inverze FFT. Do tohoto kroku jsem měl výsledky stejné jako Matlab (výpočet fft, výpočet vektoru h -> modifikace vstupního vektoru pro ifft), ale u výpočtu inverze se mi již výsledky rozchází - pravděpodobně neumím správně použít nastavení nebo něco opomíjím - bohužel moje angličtina je velice chabá a co se týče DSP a td nemám moc žádné zkušenosti.

Na vstupu mam sice reálná data, ty si ale převedu na komplexní s imag. části 0. Vstup jednorozměrné double pole s 59 999 prvky.

Můj zdrojový kód:

fftw_complex* in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * size);
fftw_complex* out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * size);
fftw_complex* out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * size);

for (int i = 0; i < size; i++)
{
    in[i][0] = data[i];
    in[i][1] = 0;
}

fftw_plan plan = fftw_plan_dft_1d(size, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(plan);
out[0][1] = 0;
fftw_destroy_plan(plan);

// H
vector<double> h(size, 0);
if (2*floor(size/2) == size)
{
    // even
    h[0] = 1;
    h[size/2] = 1;
    for (int i = 1; i < size/2; i++)
    {
        h[i] = 2;
    }
}
else
{
    // odd
    h[0] = 1;
    for (int i = 1; i < (size+1)/2; i++)
    {
        h[i] = 2;
    }
}

// IFFT
 // modifikace vstupnich dat
for (int i = 0; i < size; i++)
{
    out[i][0] *= h[i];
    out[i][1] *= h[i];
}

plan = fftw_plan_dft_1d(size, out, out2, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_execute(plan);
fftw_destroy_plan(plan);

Dokázal by někdo poradit co může být špatně?

Editováno 10.11.2014 19:26
 
Odpovědět 10.11.2014 19:26
Avatar
coells
Redaktor
Avatar
Odpovídá na drob.cl
coells:

Začni s jednoduchými příklady, například velikost vstupního pole 4 nebo 8 a zkontroluj si správnost výpočtu.
Také se zamysli nad FFT o jiné velikosti, než je mocnina dvou.

No a když ani to nepomůže, stačí se zamyslet na funkcí hilbert(x)

y = hilbert(x) = IDFT(DFT(x) * h)
DFT(y) = DFT(x) * h

A vzhledem k tomu, že znáš x a díky matlabu si umíš spočítat také y, bude nalezení chyby triviální záležitostí.

 
Nahoru Odpovědět 10.11.2014 20:56
Avatar
drob.cl
Člen
Avatar
drob.cl:

Zatím jsem odhalil, že pro velikost 8 mi to dává nejpřesnější výsledky - jen jsou 10x větší. U jiných velikostí už to je jiné. Jinak vše až do IDFT sedí.

Zkusím ještě dál zkoumat, bohužel jsem teď nějak časově vytížen.

 
Nahoru Odpovědět 11.11.2014 22:34
Avatar
drob.cl
Člen
Avatar
Odpovídá na coells
drob.cl:

Tak chybka odhalena, výstupní data po druhém fft_execute(...) je ještě potřeba dělit velikostí vstupního pole. Co se týče té velikosti, doporučuješ dávat tomu pole o velikosti 2n?

 
Nahoru Odpovědět 15.11.2014 18:17
Avatar
coells
Redaktor
Avatar
Odpovídá na drob.cl
coells:

To jsem ti psal už minule, že DFT je ortogonální, ale není ortonormální.

FFT se obvykle implementuje pouze pro velikosti 2n, občas pro velikosti f * 2n, kde f = {3,5,7}.
Pokud budeš chtít rozšiřovat pole, pak se musíš zamyslet nad peridicitou dat a vhodným způsobem doplnit vstup.

 
Nahoru Odpovědět 15.11.2014 18:36
Děláme co je v našich silách, aby byly zdejší diskuze co nejkvalitnější. Proto do nich také mohou přispívat pouze registrovaní členové. Pro zapojení do diskuze se přihlas. Pokud ještě nemáš účet, zaregistruj se, je to zdarma.

Zobrazeno 5 zpráv z 5.