NOVINKA: Získej 40 hodin praktických dovedností s AI – ZDARMA ke každému akreditovanému kurzu!
Mezinárodní den IT společnosti je tady! Pouze nyní můžeš získat 90 % extra kreditů při nákupu od 1199 kreditů s promo kódem AJTACI90. Tak neváhej!

Diskuze: IFFT pomocí FFTW

V předchozím kvízu, Online test znalostí C++, jsme si ověřili nabyté zkušenosti z kurzu.

Jak se ti líbí článek?
Před uložením hodnocení, popiš prosím autorovi, co je špatněZnaků 0 z 50-500
Jak se ti kurz líbí?
Tvé hodnocení kurzuZnaků 0 z 50-500
Aktivity
Avatar
drob.cl
Člen
Avatar
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
Avatar
coells
Tvůrce
Avatar
Odpovídá na drob.cl
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
Avatar
drob.cl
Člen
Avatar
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
Avatar
drob.cl
Člen
Avatar
Odpovídá na coells
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 2n?

 
Nahoru Odpovědět
15.11.2014 18:17
Avatar
coells
Tvůrce
Avatar
Odpovídá na drob.cl
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 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.