Poprawa szybkości implementacji FFT

Question

Jestem początkującym programistą i obecnie próbuję pracować nad projektem wymagającym szybkiej transformacji Fouriera.

ja do tej pory udało się wdrożyć następujące:

Czy ktoś ma żadnych alternatyw i propozycje w celu poprawy szybkości programu, nie tracąc na dokładność.

short FFTMethod::FFTcalc(short int dir,long m,double *x,double *y) 
{ 
long n,i,i1,j,k,i2,l,l1,l2; 
double c1,c2,tx,ty,t1,t2,u1,u2,z; 

/* Calculate the number of points */ 
n = 1; 
for (i=0;i<m;i++) 
    n *= 2; 

/* Do the bit reversal */ 
i2 = n >> 1; 
j = 0; 
for (i=0;i<n-1;i++) { 
    if (i < j) { 
    tx = x[i]; 
    ty = y[i]; 
    x[i] = x[j]; 
    y[i] = y[j]; 
    x[j] = tx; 
    y[j] = ty; 
    } 
    k = i2; 
    while (k <= j) { 
    j -= k; 
    k >>= 1; 
    } 
    j += k; 
} 

/* Compute the FFT */ 
c1 = -1.0; 
c2 = 0.0; 
l2 = 1; 
for (l=0;l<m;l++) { 
    l1 = l2; 
    l2 <<= 1; 
    u1 = 1.0; 
    u2 = 0.0; 
    for (j=0;j<l1;j++) { 
    for (i=j;i<n;i+=l2) { 
     i1 = i + l1; 
     t1 = u1 * x[i1] - u2 * y[i1]; 
     t2 = u1 * y[i1] + u2 * x[i1]; 
     x[i1] = x[i] - t1; 
     y[i1] = y[i] - t2; 
     x[i] += t1; 
     y[i] += t2; 
    } 
    z = u1 * c1 - u2 * c2; 
    u2 = u1 * c2 + u2 * c1; 
    u1 = z; 
    } 
    c2 = sqrt((1.0 - c1)/2.0); 
    if (dir == 1) 
    c2 = -c2; 
    c1 = sqrt((1.0 + c1)/2.0); 
    } 

/* Scaling for forward transform */ 
if (dir == 1) { 
    for (i=0;i<n;i++) { 
     x[i] /= n; 
     y[i] /= n; 
    } 
} 


    return(1); 
} 

Answer 1

Niedawno znalazłem ten znakomity PDF na Construction of a high performance FFTs przez Eric Postpischil. Po samodzielnym opracowaniu kilku FFT wiem, jak trudno konkurować z bibliotekami komercyjnymi. Uwierz mi, że masz się dobrze, jeśli twój FFT jest tylko 4x wolniejszy niż Intel lub FFTW, a nie 40x! Możesz jednak rywalizować, a oto jak.

Podsumowując ten artykuł, autor stwierdza, że ​​FFT Radix2 są proste, ale nieefektywne, najbardziej wydajną konstrukcją jest FFT. Jeszcze bardziej wydajną metodą jest Radix8, jednak często nie pasuje do rejestrów na CPU, więc Radix4 jest preferowany.

FFT można budować etapami, więc aby obliczyć FFT 1024 point można wykonać 10 etapów Radix2 FFT (jako 2^10 - 1024) lub 5 etapów Radix4 FFT (4^5 = 1024) . Można nawet obliczyć FFT o 1024 punktach w etapach 8 * 4 * 4 * 4 * 2, jeśli tak zdecydujesz.Mniejsza liczba etapów oznacza mniejszą liczbę odczytów i zapisów w pamięci (wąskim gardłem dla wydajności FFT jest przepustowość pamięci), dlatego dynamiczne wybieranie radixów 4, 8 lub wyższych jest koniecznością. Etap Radix4 jest szczególnie wydajny, ponieważ wszystkie wagi wynoszą 1 + 0i, 0 + 1i, -1 + 0i, 0-1i, a kod motyla Radix4 można zapisać tak, aby zmieścił się całkowicie w pamięci podręcznej.

Po drugie, każdy etap w FFT nie jest taki sam. Pierwszy stopień wagi jest równy 1 + 0i. nie ma sensu obliczanie tego ciężaru, a nawet pomnożenie przez niego, ponieważ jest to pomnożenie złożone przez 1, więc pierwszy etap można wykonać bez ciężarków. Ostatni etap może być również traktowany inaczej i może być użyty do przeprowadzenia Decymacji w Czasie (odwrócenie bitowe). Dokument Erica Postpischila obejmuje to wszystko.

Wagi można wstępnie obliczyć i zapisać w tabeli. Obliczenia sin/cos zajmują około 100-150 cykli każdy na sprzęcie x86, więc ich wstępne obliczenia mogą zaoszczędzić 10-20% całkowitego czasu obliczeń, ponieważ dostęp do pamięci jest w tym przypadku szybszy niż obliczenia procesora. Używanie szybkich algorytmów do obliczania sincos w jednym przejściu jest szczególnie korzystne (zauważ, że cos jest równe sqrt (1.0 - sinus * sinus), lub przy użyciu wyszukiwań tabel, cos jest po prostu przesunięciem fazowym sinusa).

W końcu, gdy masz już super usprawnioną implementację FFT, możesz wykorzystać wektoryzację SIMD do obliczenia 4x operacji zmiennoprzecinkowych lub 2x podwójnych operacji zmiennoprzecinkowych na cykl w ramach procedury motylkowej, aby uzyskać kolejną poprawę prędkości o 100-300%. Biorąc wszystkie powyższe, będziesz miał całkiem sprytny i szybki FFT!

Aby przejść dalej, można przeprowadzić optymalizację w locie, zapewniając różne implementacje etapów FFT skierowanych do określonych architektur procesorów. Rozmiar pamięci podręcznej, liczba rejestrów, zestawy instrukcji SSE/SSE2/3/4 itp. Różnią się w zależności od maszyny, dlatego wybór jednego rozmiaru pasuje do wszystkich podejść jest często bity przez ukierunkowane procedury. W FFTW na przykład wiele mniejszych rozmiarów FFT to wysoce zoptymalizowane rozwijane (bez pętli) implementacje ukierunkowane na konkretną architekturę. Łącząc te mniejsze konstrukcje (takie jak procedury RadixN), możesz wybrać najszybszą i najlepszą procedurę dla danego zadania.

Answer 2

Chociaż nie mogę dać wskazówkę wydajności teraz, chciałbym dać kilka rad dla optymalizacji, który jest zbyt długo na komentarz:

1. Jeśli nie masz zrobić tak, napisz kilka testów poprawności dla twojego kodu właśnie teraz. Proste testy, takie jak "wykonaj FFT z tej tablicy i sprawdź, czy wyniki pasują do podanych przeze mnie wyników" wystarczają, ale przed zoptymalizowaniem kodu potrzebujesz solidnego i zautomatyzowanego testu jednostki, który potwierdza, że ​​zoptymalizowany kod jest poprawny.
2. Następnie wpisz swój kod, aby sprawdzić, gdzie znajduje się wąskie gardło. Chociaż podejrzewam, że najgłębsza pętla to for (i=j;i<n;i+=l2) {, widzenie jest lepsze niż wiara.

Answer 3

To wygląda na podstawową implementację FFT-2 ​​FFT prosto ze starego podręcznika. Istnieje wiele dziesiątek dziesięcioleci na temat optymalizacji FFT na różne sposoby, w zależności od wielu czynników. Na przykład, czy Twoje dane są mniejsze niż pamięć podręczna procesora?

Dodano: Na przykład, jeśli wektor danych plus tabela współczynników zmieści się w CPU dcache i/lub jeśli mnożniki będą znacznie wolniejsze niż dostęp do pamięci na CPU, wówczas wstępne obliczenie tabeli współczynników Twiddle może zmniejszyć całkowity cykl liczyć na wielokrotne użycie FFT. Ale jeśli nie, wstępne przetwarzanie może być wolniejsze. Reper. YMMV.

Answer 4

Istnieje kilka rzeczy mogę polecić stara:

1. Nie zamienić elementy napędowe, zamiast obliczania indeksu bit-odwrócone. Pozwoli to zaoszczędzić wiele odczytów i zapisów w pamięci.
2. Należy wstępnie obliczyć współczynniki, jeśli wykonujesz wiele FFT o tym samym rozmiarze. Pozwoli to zaoszczędzić niektóre obliczenia.
3. Użyj radix-4 FFT zamiast radix-2. Spowoduje to mniej iteracji w wewnętrznych pętlach.

Ostateczną odpowiedź można oczywiście znaleźć, profilując kod.