Diskuze: IFFT pomocí FFTW

drob.cl

Člen

drob.cl:10.11.2014 19:26

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

coells

Tvůrce

coells:10.11.2014 20:56

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

drob.cl

Člen

drob.cl:11.11.2014 22:34

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

drob.cl

Člen

drob.cl:15.11.2014 18:17

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 2ⁿ?

Nahoru Odpovědět

15.11.2014 18:17

coells

Tvůrce

coells:15.11.2014 18:36

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

FFT se obvykle implementuje pouze pro velikosti 2ⁿ, občas pro velikosti f * 2ⁿ, 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.