Szybka transformata Fouriera Radix-2

FFT czyli dyskretna transformata Fouriera w wersji szybkiej.

Dyskretną transformatę Fouriera definiujemy jako:

$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i 2 \pi \frac{k}{N} n}$

Definicja odwrotnej dyskretnej transformaty jest bardzo podobna:

$x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k \cdot e^{+i 2 \pi \frac{k}{N} n}$

W praktyce obie transformaty korzystają z tych samych procedur z różnie zadanymi parametrami. My będziemy tylko zajmować się FFT.

Wersja naiwna wprost z definicji:

public static Complex[] fft(Complex[] a_ar, int a_length = -1)
{
    if (a_length == -1)
        a_length = a_ar.Length;

    Complex[] result = new Complex[a_length];

    for (int k = 0; k < result.Length; k++)
    {
        for (int n = 0; n < a_ar.Length; n++)
        {
            double phase = -2 * MathExtensions.PI * k * n / a_ar.Length;
            result[k] += a_ar[n] * new Complex(Math.Cos(phase), Math.Sin(phase));
        }
    }

    return result;
}

Czas wykonania: 55s. Złożoność obliczeniowa to $O(n^2)$. Test to FFT w dwóch wymiarach.

Pierwszą rzeczą jaką możemy wykonać to zrezygnować z wyliczania za każdym razem wartości new Complex(Math.Cos(phase), Math.Sin(phase)). Zauważmy, że:

$coeff_{n+1}=e^{-i 2 \pi \frac{k}{N} (n+1)}=e^{-i 2 \pi \frac{k}{N} n}e^{-i 2 \pi \frac{k}{N} }=coeff_n \cdot e^{-i 2 \pi \frac{k}{N} }$

Implementacja w kodzie:

public static Complex[] fft(Complex[] a_ar, int a_length = -1)
{
    if (a_length == -1)
        a_length = a_ar.Length;

    Complex[] result = new Complex[a_length];

    for (int k = 0; k < result.Length; k++)
    {
        Complex step = Complex.FromPolar(-2 * MathExtensions.PI * k / a_ar.Length);
        Complex coeff = Complex.ONE;

        for (int n = 0; n < a_ar.Length; n++)
        {
            result[k] += a_ar[n] * coeff;
            coeff *= step;
        }
    }

    return result;
}

Czas wykonania: 16s. Całkiem niezłe przyspieszenie.

Czas na FFT w dwóch wymiarach. Jeśli zbiór danych potraktujemy jako tablicę dwuwymiarową, to FFT wykonujemy na rzędach, a później kolumnach, lub odwrotnie. Kod:

public static Complex[,] fft2(Complex[,] a)
{
    int w = a.GetLength(0);
    int h = a.GetLength(1);

    Complex[,] result = new Complex[w, h];

    for (int y = 0; y < h; y++)
    {
        Complex[] row = new Complex[w];

        for (int x = 0; x < w; x++)
            row[x] = a[x, y];

        Complex[] frow = fft(row);

        for (int x = 0; x < w; x++)
            result[x, y] = frow[x];
    }

    for (int x = 0; x < w; x++)
    {
        Complex[] col = new Complex[h];

        for (int y = 0; y < h; y++)
            col[y] = result[x, y];

        Complex[] fcol = fft(col);

        for (int y = 0; y < h; y++)
            result[x, y] = fcol[y];
    }

    return result;
}

Wynik zwracany przez fft(), fft2() warto przesunąć, tak by element o f=0 znajdował się w środku, a nie na początku tablicy.

public static T[,] fftshift<T>(T[,] a_ar)
{
    int deltax = (a_ar.GetLength(0) / 2.0).Floor();
    int deltay = (a_ar.GetLength(1) / 2.0).Floor();

    T[,] result = new T[a_ar.GetLength(0), a_ar.GetLength(1)];

    for (int x = 0; x < a_ar.GetLength(0); x++)
    {
        int cx = OverlapOverlayCorrector.Correct_XY(x + deltax, a_ar.GetLength(0));

        for (int y = 0; y < a_ar.GetLength(1); y++)
        {
            int cy = OverlapOverlayCorrector.Correct_XY(y + deltay, 
               a_ar.GetLength(1));
            result[cx, cy] = a_ar[x, y];
        }
    }

    return result;
}

Warto zwrócić uwagę, że dla parzystej długości tablicy, podwójne wywołanie fftshift() daje dane początkowe. Tak nie jest w przypadku tablicy o długości nieparzystej. Dla odwracalności definiujemy:

public static T[,] ifftshift<T>(T[,] a_ar)
{
    int deltax = (a_ar.GetLength(0) / 2.0).Ceiling();
    int deltay = (a_ar.GetLength(1) / 2.0).Ceiling();
    return generic_shift(a_ar, deltax, deltay);
}

Kompletny kod funkcji na przesuwanie:

public static T[] fftshift<T>(T[] a_ar)
{
    int delta = (a_ar.Length / 2.0).Floor();
    return generic_shift(a_ar, delta);
}

public static T[] ifftshift<T>(T[] a_ar)
{
    int delta = (a_ar.Length / 2.0).Ceiling();
    return generic_shift(a_ar, delta);
}

private static T[] generic_shift<T>(T[] a_ar, int delta)
{
    T[] result = new T[a_ar.Length];

    for (int i = 0; i < a_ar.Length; i++)
    {
        int ci = OverlapOverlayCorrector.Correct_XY(i + delta, a_ar.Length);
        result[ci] = a_ar[i];
    }

    return result;
}

public static T[,] fftshift<T>(T[,] a_ar)
{
    int deltax = (a_ar.GetLength(0) / 2.0).Floor();
    int deltay = (a_ar.GetLength(1) / 2.0).Floor();
    return generic_shift(a_ar, deltax, deltay);
}

public static T[,] ifftshift<T>(T[,] a_ar)
{
    int deltax = (a_ar.GetLength(0) / 2.0).Ceiling();
    int deltay = (a_ar.GetLength(1) / 2.0).Ceiling();
    return generic_shift(a_ar, deltax, deltay);
}

private static T[,] generic_shift<T>(T[,] a_ar, int deltax, int deltay)
{
    T[,] result = new T[a_ar.GetLength(0), a_ar.GetLength(1)];

    for (int x = 0; x < a_ar.GetLength(0); x++)
    {
        int cx = OverlapOverlayCorrector.Correct_XY(x + deltax, a_ar.GetLength(0));

        for (int y = 0; y < a_ar.GetLength(1); y++)
        {
            int cy = OverlapOverlayCorrector.Correct_XY(y + deltay, 
                a_ar.GetLength(1));
            result[cx, cy] = a_ar[x, y];
        }
    }

    return result;
}

Podana wcześniej wersja fft2() jest dla liczb zespolonych, tymczasem nasza tablica wejściowa to zbiór liczb rzeczywistych. Tak więc stwórzmy osobną wersję na to:

public static Complex[,] fft2(Complex[,] a)
{
    int w = a.GetLength(0);
    int h = a.GetLength(1);

    Complex[,] result = new Complex[w, h];

    for (int y = 0; y < h; y++)
    {
        Complex[] row = new Complex[w];

        for (int x = 0; x < w; x++)
            row[x] = a[x, y];

        Complex[] frow = fft(row);

        for (int x = 0; x < w; x++)
            result[x, y] = frow[x];
    }

    for (int x = 0; x < w; x++)
    {
        Complex[] col = new Complex[h];

        for (int y = 0; y < h; y++)
            col[y] = result[x, y];

        Complex[] fcol = fft(col);

        for (int y = 0; y < h; y++)
            result[x, y] = fcol[y];
    }

    return result;
}

Naiwna implementacja FFT dla liczb rzeczywistych:

public static Complex[] fft(double[] a_ar, int a_length = -1)
{
    Complex[] ar = new Complex[a_ar.Length];
    for (int i = 0; i < a_ar.Length; i++)
        ar[i] = new Complex(a_ar[i]);
    return fft(ar, a_length);
}

Dla takich funkcji jak wyżej jest zmierzony czas 15s podany wcześniej. Stwórzmy teraz wersję wyspecjalizowaną dla liczb rzeczywistych:

public static Complex[] fft(double[] a_ar, int a_length = -1)
{
    if (a_length == -1)
        a_length = a_ar.Length;

    Complex[] result = new Complex[a_length];

    for (int k = 0; k < result.Length; k++)
    {
        Complex step = Complex.FromPolar(-2 * MathExtensions.PI * k / a_ar.Length);
        Complex coeff = Complex.ONE;

        for (int n = 0; n < a_ar.Length; n++)
        {
            result[k] += a_ar[n] * coeff;
            coeff *= step;
        }
    }

    return result;
}

Jest ona identyczna z poprzednią. Czas wykonania: 14.3s.

Wersje 2D możemy zrównoleglić. Najłatwiej podział danych przeprowadzić po rzędach i kolumnach:

public static Complex[,] fft2(double[,] a)
{
    int w = a.GetLength(0);
    int h = a.GetLength(1);

    Complex[,] result = new Complex[w, h];

    Parallel.For(0, h, (y) =>
    {
        double[] row = new double[w];

        for (int x = 0; x < w; x++)
            row[x] = a[x, y];

        Complex[] frow = fft(row);

        for (int x = 0; x < w; x++)
            result[x, y] = frow[x];
    });

    Parallel.For(0, w, (x) =>
    {
        Complex[] col = new Complex[h];

        for (int y = 0; y < h; y++)
            col[y] = result[x, y];

        Complex[] fcol = fft(col);

        for (int y = 0; y < h; y++)
            result[x, y] = fcol[y];
    });

    return result;
}

Czas wykonania: 8s. Czyli skalowanie jest całkiem dobre (mam dwa rdzenie u siebie).

Dla wersji dwuwymierowej spóbujmy wstępnie policzyć nasze współczynniki coeff.

private static Complex[,] PrecalculateCoeffs(int a_in_length, int a_out_length)
{
    Complex[,] result = new Complex[a_out_length, a_in_length];

    for (int k = 0; k < a_out_length; k++)
    {
        Complex step = Complex.FromPolar(-2 * MathExtensions.PI * k / a_in_length);
        Complex coeff = Complex.ONE;

        for (int n = 0; n < a_in_length; n++)
        {
            result[k, n] = coeff;
            coeff *= step;
        }
    }

    return result;
}

Ważne jest znaczenie indeksów w tej tablicy.

Najpierw poprawmy nasze funkcje fft2 w następujący sposób, by wielkość tablicy wyjścia mogła się różnić od wejścia:

public static Complex[,] fft2(double[,] a_in, int a_width = -1, int a_height = -1)
{
    if (a_width == -1)
        a_width = a_in.GetLength(0);
    if (a_height == -1)
        a_height = a_in.GetLength(1);

    Complex[,] temp = new Complex[a_width, a_in.GetLength(1)];

    Parallel.For(0, a_in.GetLength(1), (y) =>
    {
        double[] row = new double[a_width];

        for (int x = 0; x < a_width; x++)
            row[x] = a_in[x, y];

        Complex[] frow = fft(row);

        for (int x = 0; x < a_width; x++)
            temp[x, y] = frow[x];
    });

    Complex[,] result2 = new Complex[a_width, a_height];

    Parallel.For(0, a_width, (x) =>
    {
        Complex[] col = new Complex[a_height];

        for (int y = 0; y < a_height; y++)
            col[y] = temp[x, y];

        Complex[] fcol = fft(col);

        for (int y = 0; y < a_height; y++)
            result2[x, y] = fcol[y];
    });

    return result2;
}

Wersja fft() korzystająca z prekalkulowanej tablicy:

private static Complex[] fft(Complex[] a_ar, Complex[,] a_coeffs)
{
    Complex[] result = new Complex[a_coeffs.GetLength(0)];

    for (int k = 0; k < result.Length; k++)
    {
        for (int n = 0; n < a_ar.Length; n++)
        {
            result[k] += a_ar[n] * a_coeffs[k, n];
        }
    }

    return result;
}

No i wreście nowa wersja fft2():

public static Complex[,] fft2(double[,] a_in, int a_width = -1, int a_height = -1)
{
    if (a_width == -1)
        a_width = a_in.GetLength(0);
    if (a_height == -1)
        a_height = a_in.GetLength(1);

    Complex[,] temp = new Complex[a_width, a_in.GetLength(1)];
    Complex[,] coeffs = PrecalculateCoeffs(a_in.GetLength(0), a_width);

    Parallel.For(0, a_in.GetLength(1), (y) =>
    {
        double[] row = new double[a_width];

        for (int x = 0; x < a_width; x++)
            row[x] = a_in[x, y];

        Complex[] frow = fft(row, coeffs);

        for (int x = 0; x < a_width; x++)
            temp[x, y] = frow[x];
    });

    coeffs = PrecalculateCoeffs(a_in.GetLength(1), a_height);

    Complex[,] result2 = new Complex[a_width, a_height];

    Parallel.For(0, a_width, (x) =>
    {
        Complex[] col = new Complex[a_height];

        for (int y = 0; y < a_height; y++)
            col[y] = temp[x, y];

        Complex[] fcol = fft(col, coeffs);

        for (int y = 0; y < a_height; y++)
            result2[x, y] = fcol[y];
    });

    return result2;
}

Czas wykonania: 6.42.

Po pierwsze z dużo wcześniej przeprowadzonych doświadczeń nie wynika, że tablicy Complex[][] są szybsze od Complex[,]. Po drugie zmiena prekalkulacji współczynników z rekurencyjnej na bezpośrednią nic nie zmienia w dokładności wyników. Przynajmniej w moim wypadku, kiedy FFT jest konwertowana na bitmapę odcieni szarości. Ciągłe kopiowania i tworzenie tablic w fft2() też nie stanowi problemu. W naszym przypadku to około 33ms.

Dalsze przyspieszenie to już implementacja tego co właściwie nazywamy FFT. Samo FFT ma złożonośc obliczeniową $O(n\log n)$, czyli nasz testowy przypadek powinien przyspieszyć do około 0.9s. Przełączając się na alternatywną wersję poskładaną z kodu z biblioteki AlgLib rzeczywiście tak się dzieje. Ich wersja, która stara się zachować złożoność $O(n\log n)$ niezależnie od wielkości danych wejściowych dokładnie tyle się wykonuje.

Istnieje wiele wersji FFT. Głównym problemem jest tutaj wielkość tablicy wejściowej. Samo FFT polega na rekurencyjnym podziale tablicy. Najlepiej jeśli wielkość tej tablicy to wielokrotność potęgi 2 (albo 4, 8, ...). Wtedy podział przebiega gładko, wszystko ładnie mieści się cache. Dla pozostałych przypków wielkości tablicy także istnieją metody gwarantujące zachowanie złożoności $O(n\log n)$, także dla liczb pierwszych. Istnieją także dedykowane algorytmy działające od razy w dwóch albo w większej ilości wymiarów. Ich złożoność obliczeniowa nie ulega zmianie.

My zaimplementujemy tylko algorytm wymagający by wielkość tablicy była rzędu $2^N$. Pozostałe zostawimy. Nie są one dla mnie istotne. Juź i tak implemtuje FFT z czystej ciekawości, gdyż na razie wszystko mi działa z AlgLib.

Algorytm ten nazywamy RADIX-2. Z definicji DFT:

$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{\textstyle -i 2 \pi \frac{k}{N} n}$

$X_k = \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N} (2m+k)}$

$X_k = \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N} 2mk} + \sum_{m=0}^{N/2-1} x_{2m+1} \cdot e^{\textstyle -i \frac{2 \pi}{N} (2m+1)k}$

$X_k = \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk} + e^{\textstyle -i \frac{2 \pi}{N} k} \sum_{m=0}^{N/2-1} x_{2m+1} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk}$

$X_k = E_k + e^{\textstyle -i \frac{2 \pi}{N} k}O_k$

$E_k$ to DFT parzystych wyrazów ciągu, $O_k$ to DFT nieparzystych wyrazów ciągów.

Rozbiliśmy więc jedną dużą DFT na dwie małe. Cały proces możemy powtarzać, aż długość DFT osiągnie 1. Jednak z punktu prędkości obliczeń może to nie być najlepszy moment na przerwanie rekursji.

Dodatkowo zauważmy, że:

$X_{\textstyle k+\frac{N}{2}} = E_{\textstyle k+\frac{N}{2}} + e^{\textstyle -i \frac{2 \pi}{N} (k+\frac{N}{2})}O_{\textstyle k+\frac{N}{2}}$

$\begin{align*} E_{k+N/2} &= \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} m(k+\frac{N}{2})} \\ &= \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk} \cdot e^{\textstyle -i \cdot {2 \pi} m} \\ &= \sum_{m=0}^{N/2-1} x_{2m} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk} \\ &= E_k \end{align*}$

$\begin{align*} O_{k+N/2} &= e^{\textstyle -i \frac{2 \pi}{N} (k+N/2)} \sum_{m=0}^{N/2-1} x_{2m+1} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} m(k+N/2))} \\ &= e^{\textstyle-i \cdot \pi } \cdot e^{\textstyle -i \frac{2 \pi}{N} k} \sum_{m=0}^{N/2-1} x_{2m+1} \cdot e^{\textstyle -i \cdot 2\pi} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk)} \\ &= -e^{\textstyle -i \frac{2 \pi}{N} k} \sum_{m=0}^{N/2-1} x_{2m+1} \cdot e^{\textstyle -i \frac{2 \pi}{N/2} mk)} \\ &= -e^{\textstyle -i \frac{2 \pi}{N} k} \cdot O_k \end{align*}$

$X_{\textstyle k+\frac{N}{2}} = E_k - e^{\textstyle -i \frac{2 \pi}{N} k} \cdot O_k$

Czyli tak naprawdę druga połówka DFT jest pewnym przekształceniem pierwszej i nie trzeba jej liczyć wprost.

Musiałem także zmienić rozmiar obrazka z 600x600 na 512x512. Nowy czas wykonania to 4.8s.

Wywaliłem także kod odpowiedzialny za możliwość tworzenia DFT o różnej długości. DFT jest długości tablicy wejściowej.

Korzystając z tych informacji możemy napisać rekurencyjną nową wersję fft():

public static Complex[] fft_radix2(Complex[] a_ar, int a_length, int a_from, int a_step, Complex[,] a_coeffs = null)
{
    if (a_length == 1)
        return new Complex[] { a_ar[a_from] };

    int len2 = a_length / 2;

    Complex[] odds = fft_radix2(a_ar, len2, a_from, a_step * 2, a_coeffs);
    Complex[] evens = fft_radix2(a_ar, len2, a_from + a_step, a_step * 2, a_coeffs);

    Complex[] result = new Complex[a_length];

    for (int k = 0; k < result.Length / 2; k++)
    {
        Complex even = evens[k];
        Complex odd = odds[k];

        Complex temp = Complex.FromPolar(
            -2 * MathExtensions.PI * k / a_length) * odd;
        result[k] = even + temp;
        result[k + len2] = even - temp;
    }

    return result;
}

Czas wykonania to 750ms. Zwiększamy obraz do 1024x1024. Czas wykonania to 2.8s.

Wprowadźmy rekurencyjne liczenie współczynnika temp:

public static Complex[] fft_radix2(
    Complex[] a_ar, int a_length, int a_from, int a_step, 
    Complex[,] a_coeffs = null)
{
    if (a_length == 1)
        return new Complex[] { a_ar[a_from] };

    int len2 = a_length / 2;

    Complex[] odds = fft_radix2(a_ar, len2, a_from, a_step * 2, a_coeffs);
    Complex[] evens = fft_radix2(a_ar, len2, a_from + a_step, a_step * 2, a_coeffs);

    Complex[] result = new Complex[a_length];

    Complex step = Complex.FromPolar(-2 * MathExtensions.PI / a_length);
    Complex coeff = Complex.ONE;

    for (int k = 0; k < result.Length / 2; k++)
    {
        Complex even = evens[k];
        Complex odd = odds[k];

        Complex temp = coeff * odd;
        result[k] = even + temp;
        result[k + len2] = even - temp;
        coeff *= step;
    }

    return result;
}

Czas wykonania to 2.3s.

Teraz spóbujmy ustalić rozmiar tablicy przy którym powinniśmy przerwać rekurencje. Po testach okazało się, że nic to nie daje. Zarówno jeśli wykorzystamy już istniejącą wersję fft(), lub napiszemy szybką wersję dla tablicy o wielkości 2. Tak więc zostaje jak jest.

W porównaniu do AlgLib moja wersja jest ciut szybsza. Na pewno można ją próbować przyspieszyć. Mnie to co jest zadawala i wystarcza. Nie potrzebuje szybciej ani bardziej uniwersalnie.

Pełny kod:

public static class FFT
{
    #region FFT

    public static Complex[] fft(Complex[] a_ar)
    {
        if (Bits.IsBase2(a_ar.Length))
            return fft_radix2(a_ar, a_ar.Length, 0, 1);
        else
            return fft_slow(a_ar);
    }

    public static Complex[] fft(double[] a_ar)
    {
        if (Bits.IsBase2(a_ar.Length))
            return fft_radix2(a_ar, a_ar.Length, 0, 1);
        else
            return fft_slow(a_ar);
    }

    private static Complex[] fft_radix2(Complex[] a_ar, int a_length, 
        int a_from, int a_step)
    {
        if (a_length == 1)
            return new Complex[] { a_ar[a_from] };

        int len2 = a_length / 2;

        Complex[] odds = fft_radix2(a_ar, len2, a_from, a_step * 2);
        Complex[] evens = fft_radix2(a_ar, len2, a_from + a_step, a_step * 2);

        Complex[] result = new Complex[a_length];

        Complex step = Complex.FromPolar(-2 * MathExtensions.PI / a_length);
        Complex coeff = Complex.ONE;

        for (int k = 0; k < result.Length / 2; k++)
        {
            Complex even = evens[k];

            Complex temp = coeff * odds[k];
            result[k] = even + temp;
            result[k + len2] = even - temp;
            coeff *= step;
        }

        return result;
    }

    private static Complex[] fft_radix2(double[] a_ar, int a_length, 
        int a_from, int a_step)
    {
        if (a_length == 1)
            return new Complex[] { new Complex(a_ar[a_from]) };

        int len2 = a_length / 2;

        Complex[] odds = fft_radix2(a_ar, len2, a_from, a_step * 2);
        Complex[] evens = fft_radix2(a_ar, len2, a_from + a_step, a_step * 2);

        Complex[] result = new Complex[a_length];

        Complex step = Complex.FromPolar(-2 * MathExtensions.PI / a_length);
        Complex coeff = Complex.ONE;

        for (int k = 0; k < result.Length/2; k++)
        {
            Complex even = evens[k];

            Complex temp = coeff * odds[k];
            result[k] = even + temp;
            result[k + len2] = even - temp;
            coeff *= step;
        }

        return result;
    }

    private static Complex[] fft_slow(Complex[] a_ar, Complex[,] a_coeffs = null)
    {
        if (a_coeffs == null)
        {
            Complex[] result = new Complex[a_ar.Length];

            for (int k = 0; k < result.Length; k++)
            {
                Complex step = Complex.FromPolar(
                    -2 * MathExtensions.PI * k / a_ar.Length);
                Complex coeff = Complex.ONE;

                for (int n = 0; n < a_ar.Length; n++)
                {
                    result[k] += a_ar[n] * coeff;
                    coeff *= step;
                }
            }

            return result;
        }
        else
        {
            Complex[] result = new Complex[a_ar.GetLength(0)];

            for (int k = 0; k < result.Length; k++)
            {
                for (int n = 0; n < a_ar.Length; n++)
                {
                    result[k] += a_ar[n] * a_coeffs[k, n];
                }
            }

            return result;
        }
    }

    private static Complex[] fft_slow(double[] a_ar, Complex[,] a_coeffs = null)
    {
        if (a_coeffs == null)
        {
            Complex[] result = new Complex[a_ar.Length];

            for (int k = 0; k < result.Length; k++)
            {
                Complex step = Complex.FromPolar(
                    -2 * MathExtensions.PI * k / a_ar.Length);
                Complex coeff = Complex.ONE;

                for (int n = 0; n < a_ar.Length; n++)
                {
                    result[k] += a_ar[n] * coeff;
                    coeff *= step;
                }
            }

            return result;
        }
        else
        {
            Complex[] result = new Complex[a_ar.Length];

            for (int k = 0; k < result.Length; k++)
            {
                for (int n = 0; n < a_ar.Length; n++)
                {
                    result[k] += a_ar[n] * a_coeffs[k, n];
                }
            }

            return result;
        }
    }

    #endregion

    #region FFT2

    public static Complex[,] fft2(double[,] a_in)
    {
        int w = a_in.GetLength(0);
        int h = a_in.GetLength(1);

        Complex[,] temp = new Complex[w, a_in.GetLength(1)];
        Complex[,] coeffs = null;

        if (Bits.IsBase2(w))
            PrecalculateCoeffs(a_in.GetLength(0), w);

        Parallel.For(0, a_in.GetLength(1), (y) =>
        {
            double[] row = new double[w];

            for (int x = 0; x < w; x++)
                row[x] = a_in[x, y];

            Complex[] frow;
            if (Bits.IsBase2(row.Length))
                frow = fft_radix2(row, row.Length, 0, 1);
            else
                frow = fft_slow(row, coeffs);

            for (int x = 0; x < w; x++)
                temp[x, y] = frow[x];

        });

        if (Bits.IsBase2(h))
            coeffs = PrecalculateCoeffs(a_in.GetLength(1), h);

        Complex[,] result2 = new Complex[w, h];

        Parallel.For(0, w, (x) =>
        {
            Complex[] col = new Complex[h];

            for (int y = 0; y < h; y++)
                col[y] = temp[x, y];

            Complex[] fcol;
            if (Bits.IsBase2(col.Length))
                fcol = fft_radix2(col, col.Length, 0, 1);
            else
                fcol = fft_slow(col, coeffs);

            for (int y = 0; y < h; y++)
                result2[x, y] = fcol[y];
        });

        return result2;
    }

    public static Complex[,] fft2(Complex[,] a_in)
    {
        int w = a_in.GetLength(0);
        int h = a_in.GetLength(1);

        Complex[,] temp = new Complex[w, a_in.GetLength(1)];
        Complex[,] coeffs = null;

        if (Bits.IsBase2(w))
            PrecalculateCoeffs(a_in.GetLength(0), w);

        Parallel.For(0, a_in.GetLength(1), (y) =>
        {
            Complex[] row = new Complex[w];

            for (int x = 0; x < w; x++)
                row[x] = a_in[x, y];

            Complex[] frow;
            if (Bits.IsBase2(row.Length))
                frow = fft_radix2(row, row.Length, 0, 1);
            else
                frow = fft_slow(row, coeffs);

            for (int x = 0; x < w; x++)
                temp[x, y] = frow[x];

        });

        if (Bits.IsBase2(h))
            coeffs = PrecalculateCoeffs(a_in.GetLength(1), h);

        Complex[,] result2 = new Complex[w, h];

        Parallel.For(0, w, (x) =>
        {
            Complex[] col = new Complex[h];

            for (int y = 0; y < h; y++)
                col[y] = temp[x, y];

            Complex[] fcol;
            if (Bits.IsBase2(col.Length))
                fcol = fft_radix2(col, col.Length, 0, 1);
            else
                fcol = fft_slow(col, coeffs);

            for (int y = 0; y < h; y++)
                result2[x, y] = fcol[y];
        });

        return result2;
    }

    private static Complex[,] PrecalculateCoeffs(int a_in_length, int a_out_length)
    {
        Complex[,] result = new Complex[a_out_length, a_in_length];

        for (int k = 0; k < a_out_length; k++)
        {
            Complex step = Complex.FromPolar(
                -2 * MathExtensions.PI * k / a_in_length);
            Complex coeff = Complex.ONE;

            for (int n = 0; n < a_in_length; n++)
            {
                result[k, n] = coeff;
                coeff *= step;
            }
        }

        return result;
    }

    #endregion

    #region FFT SHIFT

    public static T[] fftshift<T>(T[] a_ar)
    {
        int delta = (a_ar.Length / 2.0).Floor();
        return generic_shift(a_ar, delta);
    }

    public static T[] ifftshift<T>(T[] a_ar)
    {
        int delta = (a_ar.Length / 2.0).Ceiling();
        return generic_shift(a_ar, delta);
    }

    private static T[] generic_shift<T>(T[] a_ar, int delta)
    {
        T[] result = new T[a_ar.Length];

        for (int i = 0; i < a_ar.Length; i++)
        {
            int ci = OverlapOverlayCorrector.Correct_XY(i + delta, a_ar.Length);
            result[ci] = a_ar[i];
        }

        return result;
    }

    public static T[,] fftshift<T>(T[,] a_ar)
    {
        int deltax = (a_ar.GetLength(0) / 2.0).Floor();
        int deltay = (a_ar.GetLength(1) / 2.0).Floor();
        return generic_shift(a_ar, deltax, deltay);
    }

    public static T[,] ifftshift<T>(T[,] a_ar)
    {
        int deltax = (a_ar.GetLength(0) / 2.0).Ceiling();
        int deltay = (a_ar.GetLength(1) / 2.0).Ceiling();
        return generic_shift(a_ar, deltax, deltay);
    }

    private static T[,] generic_shift<T>(T[,] a_ar, int deltax, int deltay)
    {
        T[,] result = new T[a_ar.GetLength(0), a_ar.GetLength(1)];

        for (int x = 0; x < a_ar.GetLength(0); x++)
        {
            int cx = OverlapOverlayCorrector.Correct_XY(
                x + deltax, a_ar.GetLength(0));

            for (int y = 0; y < a_ar.GetLength(1); y++)
            {
                int cy = OverlapOverlayCorrector.Correct_XY(
                    y + deltay, a_ar.GetLength(1));
                result[cx, cy] = a_ar[x, y];
            }
        }

        return result;
    }
    #endregion
}