Równanie różnicowe filtru dolnoprzepustowego pierwszego stopnia

Funkcja przejścia filtru dolnoprzepustowego 1 stopnia dana ogólnie

$H(s) = \displaystyle\frac{1}{s R C + 1}$

Korzystając z informacji w http://en.wikipedia.org/wiki/Low-pass_filter równanie różnicowe takiego filtra możemy zapisać jako:

$y(n) = b_0x(n) + b_1x(n-1) + b_2x(n-2) - a_1y(n-1) - a_2y(n-2)$

, gdzie:

$b_0 = \displaystyle\frac{T}{RC+T}$
$b_1 = 0$
$b_2 = 0$
$a_1 = b_0 - 1$
$a_2 = 0$

Jeśli skorzystamy tutaj z metody takiej jak przy filtrze drugiego stopnia to nasze wynikowe blepy będą się nieznacznie różnić. Maksymalny błąd względny będzie rzędu 0.001. Podana tutaj wyżej metoda jest pozbawiona tego błędu.

Istnieje wiele metod konwersji filtrów ciągłych w czasie na dyskretne. W zależności od charakterystyki funkcji H(s) metody te mogą dawać mniej lub bardziej dokładną postać H(z).

Kolejna uwaga: porównując wersję lsim do wersji bq widzimy, że powyższa metoda nie generuje takich samych wyników jak lsim. Takie same wyniki jak lsim generuje metoda taka jak dla drugiego rzędu. Jeśli uznać, że lsim jest najdokładniejszy, zaś metody różnicowe są mniej dokładne to metoda z filtru drugiego stopnia jest dokładniejsza. Czemu więc w WinUAE zastosowano metodę mniej dokładną.

Porównując charakterystykę częstotliwością sygnału lsim z bq wygenerowanym przez metodę taką jak dla filtra drugiego rzędu widzimy, że choć różnice w czasie są znaczne (0.001), różnice w częstotliwości są znacznie mniejsze (0.0001).

Kolejna uwaga: oryginalny kod w pythonie wykorzystuje tak zwane frequency prewarping, brak tego w kodzie C# powoduje błąd rzędu $10^{-5}$. Przy czym nie wiem czy ten błąd jest spowodowany głównie brakiem FQ PW. Tak więc wydaje się, że stosowanie tej techniki mija się tutaj z celem.

Window function - Kaiser w Python

Tutaj implementacja w Python z wykorzystaniem odwołania do scipy.special.iv, czyli Modified Bessel function of first order. Mam nadzieję, że w następnym poście podam kod na scipy.special.iv. Wtedy wszystko będzie już gotowe do przeniesienia w C#. By tam generować BLEP-y programowo jakie się nam wymarzy.

def kaiser(samples, beta):
    result = [0] * samples
    M = samples - 1
    for i in range(samples):
        num = special.iv(1, beta*math.sqrt(1-(2*i/M - 1)**2))
        den = special.iv(1, beta)
        result[i] = num / den
    return result

Transformata Z i równanie różnicowe filtru dolnoprzepustowego drugiego stopnia

Filtr dolnoprzepustowy pierwszego stopnia zostanie omówiony w innym poście. Nie możemy użyć poniższych wzorów podstawiając po prostu A=0. Ale za to podstawienie A=0 i B=0 daje poprawne wyniki.

Filtru dolnoprzepustowego drugiego stopnia:

$H(s) = \displaystyle\frac{1}{s^2 R_1 R_2 C_1 C_2 + s(R_1 C_1 + R_2 C_2) + 1}$

Ogólnie:

$H(s) = \displaystyle\frac{1}{s^2 A + s B + 1}$

Filtr pierwszego stopnia to filtr dla którego A = 0.

Korzystając z wzoru na związek pomiędzy transformatą Laplaca a Z:

$s =\displaystyle\frac{2}{T} \frac{(z-1)}{(z+1)}$

Po przekształceniach otrzymujemy:

$H(z) = \displaystyle\frac{T^2 + 2T^2z^{-1} + T^2 z^{-2}}{4A+2TB+T^2 + (2T^2-8A)z^{-1}+(4A-2TB+T^2)z^{-2}}$

Co możemy zapisać jako:

$H(z)=\displaystyle\frac{b_0+b_1z^{-1}+b_2z^{-2}} {1+a_1z^{-1}+a_2z^{-2} }$

, gdzie:

$a_0=4A+2TB+T^2$
$a_1=(2T^2-8A)/{a_0}$
$a_2=(4A-2TB+T^2)/{a_0}$
$b_0={T^2}/{a_0}$
$b_1={2T^2}/{a_0}$
$b_2={T^2}/{a_0}$

T to okres co jaki chcemy otrzymywać kolejne próbki, w naszym przypadku jest to częstotliwość z jakiej wielokrotnością custom chip jest w stanie zmieniać wyjście.

Filtr zapisany w takiej postaci to Digital biquad filter.

Takie funkcji przejścia odpowiada następujące równianie różnicowe:

$y(n) = b_0x(n) + b_1x(n-1) + b_2x(n-2) - a_1y(n-1) - a_2y(n-2)$

Wszystko to pozwala nam na generowania przebiegów BLEP wprost z równania różnicowego, którego parametry są bezpośrednio określone parametrami elektrycznymi filtra i częstotliwością z jaką chcemy ciągły sygnał próbkować.

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
}

Resampling

To kombinacja upsamplingu i downsamplingu. Pozwala nam ona zmienić próbkowanie sygnału o niecałkowitą skale. Ponieważ jest to kombinacja upsamplingu i downsamplingu dla pewnych współczynników skalowania może się okazać, że współczynnik upsamplingu może być bardzo duży. Takie przypadki w moim przykładzie staram się odrzucać. No i dla współczynników skalowania całkowitych albo ich odwrotności musimy tylko skorzystać z upsamplingu lub downsamplingu.

Funkcja:

function [ tout, fout ] = resampling( fin, fsin, fsout )

m = lcm(fsin, fsout);
u = m / fsin;
d = m / fsout;

if (u > 1000)
    error('up scale factor too big');
end

if (u > 1)
    [tout, fout] = upsampling(fin, fsin, fsin*u);
    if (d > 1)
        [tout, fout] = downsampling(fout, fsin*u, fsout);
    end
elseif (d > 1)
    [tout, fout] = downsampling(fin, fsin, fsout);
else
    tout = (1/fsin) * linspace(0, size(fin, 2), size(fin, 2)); 
    fout = fin;
end

end

Ograniczenie na błąd powinniśmy sobie sami ustawić, w zależności jakiej długości sygnały przetwarzamy i ile mamy pamięci, ... czasu.

Upsampling

Upsampling czyli zwiększenie częstotliwości próbkowania sygnału o zadany całkowity współczynik przy zachowaniu jego widma.

Funkcja:

function [ tout, fout ] = upsampling( fin, fsin, fsout )

scale = fsout / fsin;

if (fix(scale) ~= scale)
    error('scale should be integer');
end

if (scale <= 0)
    error('scale should be greater than zero');
end

fin = fin * scale;
fout = upsample(fin, scale);
tout =  linspace(0, (1 / fsout) * (size(fout, 2) - 1), size(fout, 2)); 

[~, k] = sinc_gen(fsout, fsin/2);
if (mod(size(k, 2), 2) == 0)
    [~, k] = sinc_gen(fsout, fsin/2, fix(size(k, 2) / 2) * 2 + 1);
end

fout = conv(fout, k, 'same');

fout = fout(1:1:end-scale+1);
tout = tout(1:1:end-scale+1);

end

Funkcja upsample uzupełnia fin zerami. Zera w ilości scale-1 są wstawiane między próbki i na koniec. My te końcowe zera usuwamy. Dzięki czemu ostatnia próbka sygnału wyjściowego odpowiada ostatniej próbce sygnału wejściowego.

Okno sinc musi mieć nieparzystą długość by uniknącć przesunięcia wyniku.

Sam splot to tzw. idealna rekonstrukcja, tutaj o tyle prostsza i szybsza, że współczynnik skalowania jest całkowity i sinc rekonstruujący może być wyliczony raz i używamy dla każdego punktu z osobna.

W stosunku do wbudowanej w Matlab funkcji interp błąd przez nas popełniany jest na tyle niewielki, że możemy powiedzieć, że wynika on z innej metody interpolacji.

Przykład:

%% input data

fsin = 6000;
fsout = 30000;
scale = fsout / fsin;

func = @(ti) sin(2*pi*650*ti-pi/4) + sin(2*pi*1500*ti+pi/2);
tend = 1/500*2;

% build input func.
Tin = 1/fsin;
tin = (0:Tin:tend);
fin = func(tin);

% reconstruct
[fout, tout] = upsampling(fin, fsin, fsout);

% calculate ideal reconstruction
forg = func(tout);

fout1 = interp(fin, scale);
fout1 = fout1(1:1:end-scale+1);

%% compare them

hold on;
plot(tout, forg, 'r');
plot(tin, fin, 'g');
plot(tout, fout, 'b');
plot(tout, fout1, 'm');
legend('ideal','in', 'out', 'interp');

%% compare fq

FIN = fftshift(abs(fft(fin)));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));
s = size(FIN, 2) / size(FOUT, 2);
FORG = s*fftshift(abs(fft(forg)));
FOUT = s*FOUT;
FOUT1 = s*FOUT1;

fqout = fsout/2*linspace(-1, 1, size(FORG, 2));
fqin = (1/Tin)/2*linspace(-1, 1, size(FIN, 2));

hold on;
plot(fqout, FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, FOUT, 'b');
plot(fqout, FOUT1, 'm');
legend('ideal','in', 'out', 'interp');

%% error

error_calc(fout, forg);
error_calc(fout1, forg);

Downsampling

Downsampling czyli zmniejszenie częstotliwości próbkowania sygnału o zadany całkowity współczynnik przy zachowaniu jego widma (poniżej częstotliwości do której zmniejszamy). Ponieważ tutaj rekonstrukcja odbywa się z częstotliwością mniejszą o całkowitą wartość w przeciwieństwie do rekonstrukcji wystarczy nam tylko jeden impuls sinc. Cały kod jest znacznie szybszy.

Funkcja:

function [ fout, tout ] = downsampling( fin, fsin, fsout )

scale = fsin / fsout;

if (fix(scale) ~= scale)
    error('scale should be integer');
end

if (scale <= 0)
    error('scale should be greater than zero');
end

% works only with odd window.
k = sinc_gen(fsin, fsout/2);
if (mod(size(k, 2), 2) == 0)
    k = sinc_gen(fsin, fsout/2, size(k, 2) + 1);
end

tout = linspace(0, (1/fsin) * (size(fin, 2)-1), ... 
                ceil(size(fin, 2)/scale));
fout = zeros(1, size(tout, 2));

for m=1:1:size(fout, 2)
    fout(m) = conv_point(fin, k, (m - 1) * scale + 1);
end

end

Przykład:

%% input data

% resampling
fsin = 30000;
fsout = 6000;
scale = fsin / fsout;

func = @(ti) sin(2*pi*650*ti-pi/4) + sin(2*pi*1500*ti+pi/2);
tend = 1/500*2.4; % czym dluzszy tym lepsze widmo

% build input func.
Tin = 1/fsin;
tin = (0:Tin:tend);
fin = func(tin);

% downsampling
[fout, tout] = downsampling(fin, fsin, fsout);

fout1 = decimate(fin, scale);

% calculate ideal reconstruction
forg = func(tout);

%% compare them

fig = figure;
hold on;
plot(tout, forg, 'r');
plot(tin, fin, 'g');
plot(tout, fout, 'b');
plot(tout, fout1, 'm');
legend('ideal','in', 'out', 'decimate');

%% compare fq

FIN = fftshift(abs(fft(fin)));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));
s = size(FIN, 2) / size(FOUT, 2);
FORG = s*fftshift(abs(fft(forg)));
FOUT = s*FOUT;
FOUT1 = s*FOUT1;

fqout = fsout/2*linspace(-1, 1, size(FORG, 2));
fqin = (1/Tin)/2*linspace(-1, 1, size(FIN, 2));

hold on;
plot(fqout, FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, FOUT, 'b');
plot(fqout, FOUT1, 'm');
legend('ideal','in', 'out', 'decimate');

%% error;

error_calc(fout, forg);
error_calc(fout1, forg);

Błąd mojej wersji jest porównywalny do tej wbudowanej w Matlab. Ilość próbek wejściowych jest równa wielokrotnością skalowania. Jak widać obie wersje rozmijają się, Matlab na początku, ja na końcu. Dla innych liczb próbek obie funkcję potrafią zwrócić bardzo porównywalne wyniki.

Idealna rekonstrukcja sygnału

Mamy sygnał fin, który został spróbkowany z częstotliwością fsin. Chcemy wygenerować sygnał, który został spróbkowany z częstotliwością fsout. Jeśli fsout > fsin w sygnale wyjściowym nie powinno być żadnej częstotliwości powyżej fsin. Jeśli fsout < fsin w sygnale wyjściowym powinny się znaleźć tylko częstotliwości mniejsze od fsout z sygnału wejściowego.

Wszelkie interpolacje (spline, bicubic, itp) w dziedzinie czasu odpadają, gdyż wprowadzimy tym częstotliwości do sygnału, których w nim nie było oryginalnie. Jedne sprawują się gorzej, inne lepiej, ale generalnie interpolacje odtworzą sygnał wizualnie niby dobrze, częstotliwościowo nie będzie to to samo. Ma to szczególne znaczenie kiedy pracujemy z sygnałem dźwiękowym. No dobra za wyjątkiem podanej niżej interpolacji. Rekonstrukcję możemy też przeprowadzić z dziedzinie częstotliwości.

Aby odtworzyć oryginalny sygnał potrzeba nam minimum dwie próbki na na okres, w praktyce powinno być ich więcej. Możemy próbkować dokładnie w miejscach gdzie sinusoida przechodzi przez zero. A nawet jeśli jesteśmy minimalnie minięci podczas próbkowania z zerem, pobierane przez nas wartości są małe i błąd przez to jest duży.

Teoria mówi, że sygnał oryginalny może zostać dokładnie zrekonstruowany ze swych próbek, jeśli próbki były pobierane w odstępach czasu T, zaś sygnał próbkowany nie zawiera częstotliwości większych niż połowa częstotliwości próbkowania. Czyli jeśli próbkujemy z częstotliwością 44100Hz, to sygnał nie powinien zawierać częstotliwości większych niż 22050Hz. W przeciwnym razie częstotliwości wyższe od limitu pojawią się w sygnale zrekonstruowanym zniekształcając go (aliasign). Twierdzenie Whittaker–Shannon'a mówi:

$x(t) = \sum_{n=-\infty}^{\infty} x[n] \cdot {\rm sinc}\left(\frac{t - nT}{T}\right)\,$

Inaczej sygnał oryginalny to splot impulsów z próbek z funkcją sinc o odpowiednim charakterze.

Teraz mając odtworzony sygnał ciągły, próbkujemy go ponownie z nową częstotliwością. Jeśli z zmniejszą to najpierw pozbywamy się z niego częstotliwości większych od połowy nowej częstotliwości próbkowania.

Całość robi poniższy kod. Jeśli mamy zrekonstruować sygnał z niższą częstotliwością to filtrowane są próbki wejściowe. Taki sygnał podlega rekonstrukcji zgodnie z wyżej podaną formułą tylko w interesujących nas czasach.

Kod funkcji:

function [ fout, tout ] = reconstruction( fin, fsin, fsout )

Tout = 1 / fsout;
Tin = 1 / fsin;

if (fsin == fsout)
    fout = fin;
    tout = Tout* linspace(0, size(fout, 2), size(fout, 2)); 
    return;
end

scale = fsout / fsin;

if (scale < 1)
    % for compatibility with downsampling and decimate
    fout = zeros(1, ceil(size(fin, 2)* scale));
else
    fout = zeros(1, fix(size(fin, 2) * scale));
end

tout = linspace(0, (1/fsin) * (size(fin, 2) - 1), size(fout, 2));            

% when downsampling cut frequencies below fsout/2 from fin 
% function.
if (fsin > fsout)
    k = sinc_gen(fsin, fsout/2);
    fin = conv(fin, k, 'same');
end

W = fix(160 / (fsout*2 / fsin)); % experimental
W = min(W, size(fin, 2));

tin = (0:1:size(fin, 2)-1) * Tin;

if (W > size(fin, 2))
    W = size(fin, 2);
end

for m=1:size(fout, 2)
    [~, min_index] = min(abs(tout(m) - tin));
    min_index = fix(min_index - W/2);
    if (min_index < 1)
        min_index = 1;
    end
    max_index = min_index + W - 1;
    if (max_index > size(fin, 2))
        max_index = size(fin, 2);
    end
    sc = sinc((tout(m) - (min_index-1:1:max_index-1)*Tin) / Tin);
    fout(m) = sum(fin(min_index:1:max_index) .* sc);
end

% for compatibility with upsampling
if (scale > 1)
    fout = fout(1:1:end-scale+1);
    tout = tout(1:1:end-scale+1);
end

Dla idealnej rekonstrukcji okno W powinno być równe długości funkcji, w praktyce może być znacznie krótsze. Jeśli mamy do czynienia z downsamplingiem to najpierw sygnał trzeba przefiltrować filtrem dolnoprzepustowym. Inaczej częstotliwości powyżej fsout, które były w fin odbiją się i pojawią się jako zakłócenia w sygnale. Z tak zakłóconym sygnałem już niż nie możemy zrobić - nie da się stwierdzić co jest oryginałem, a co zakłóceniem.

Rekonstrukcja taka jest dokładna, ale wyjątkowo czasochłonna, dla każdej próbki wyjściowej trzeba policzyć osobny impuls sinc.

Normalnie okno sinc powinno być równe wektorowi wejściowemu, tutaj jest krótsze, o długości powiedzmy eksperymentalnej dla której nie widać różnicy.

Jest też naniesionych parę poprawek na kompatybilność z innymi funkcjami dotyczącymi długości wektora wyjściowego.

Przykład:

%% input data

% resampling
fsin = 30000;
fsout = 7300;
Tin = 1/fsin;
scale = fsin / fsout;

func = @(ti) sin(2*pi*650*ti-pi/4) + sin(2*pi*1500*ti+pi/2);
tend = 1/500*2;

%build input func.
tin = (0:Tin:tend);
fin = func(tin);

% reconstruct
[fout, tout] = reconstruction(fin, fsin, fsout);

%resample
[fout1, tout1] = resampling(fin, fsin, fsout);

% calculate ideal reconstruction
forg = func(tout);

%% compare them

fig = figure;
hold on;
plot(tout, forg, 'r');
plot(0:Tin/20:tend, func(0:Tin/20:tend), 'g');
plot(tout, fout, 'b');
plot(tout1, fout1, 'black');
legend('ideal','in', 'out reconstr', 'out resampl');

%% compare fq

FORG = fftshift(abs(fft(forg)));
FIN = fftshift(abs(fft(func(0:Tin:tend))));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));

fqout = fsout/2*linspace(-1, 1, size(FORG, 2));
fqout1 = fsout/2*linspace(-1, 1, size(FOUT1, 2));
fqin = (1/Tin)/2*linspace(-1, 1, size(FIN, 2));

s = size(FIN, 2) / size(FOUT, 2);

hold on;
plot(fqout, s*FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, s*FOUT, 'b');
plot(fqout1, s*FOUT1, 'black');
legend('ideal','in', 'out reconstr', 'out resampl');

%% error

error_calc(fout, forg);
error_calc(fout1, forg);

Filtr FIR o minimalnym przesunięciu fazowym

Weźmy wszystkie filtry o takiej samej charakterystyce częstotliwościowej G, ale różniących się charakterystyką przesunięcia fazowe. Filtr o minimalnym przesunięciu fazowym to taki filtr, że w dowolnym okresie czasu $<0, t>$ energia sygnału wyjściowego jest większa lub równa dowolnemu innemu filtrowi z G. Czyli sygnał w takim filtrze najszybciej narasta i najszybciej maleje. Sygnał jest maksymalnie przesunięty ku zeru.

Dla filtrów FIR Matlab posiada funkcję do znajdowania MPF (minimum-phase filter) rceps. Niestety nie działa ona za dobrze. Rezultat nie wygląda tak jak powinien, a na dodatek dla niektórych długości wektora policzenie wyniku jest niemożliwe.

Moja wersja:


function [ minblip ] = fir_min_phase( impulse )

impulse1 = [impulse, zeros(1, size(impulse, 2)*7)];

y = ifft(log_safe(abs(fft(impulse1))));

%y(1) = y(1) * 1;
y(2:fix(end/2)) = y(2:fix(end/2)) * 2;
%y(fix(end/2+1) = y(fix(end/2+1) * 1;
y(fix(end/2)+2:end) = 0;

minblip = real(ifft(exp(fft(y))));

minblip = minblip(1:1:size(impulse, 2));

minblip = minblip / trapz(minblip);

end

function [ r ] = log_safe( v )

r = log(v);

for k=1:1:size(r, 2)
if (r(k) == -inf)
r(k) = -1e-100;
end
end

end


Kompletna magia dla mnie - matematycznie nie wiem jak to działa. Zmniejszając początkową rozdzielczość okna, albo źle dobierając okno +/-1 możemy otrzymać wersję Matlabową.

Porównanie impulsu sinc i jego wersji MPF.


fs = 7093790/2;
d1 = firwin(fs, 21000, 2047);
d2 = fir_min_phase(d1);
[~, d3] = rceps(d1);

hold on;
plot(d1, 'r');
plot(d2, 'g');
plot(d3, 'b');
legend('org', 'my', 'matlab');

D1 = abs(fftshift(fft(d1)));
D2 = abs(fftshift(fft(d2)));
D3 = abs(fftshift(fft(d3)));
fx = linspace(-fs/2, fs/2, 2047);

hold on;
plot(fx, D1, '*');
plot(fx, D2, 'g');
plot(fx, D3, 'o');
legend('org', 'my', 'matlab');


Widzimy, że we wszystkich 3 przypadkach charakterystyka częstotliwościowa filtrów jest taka sama.

Idealny filtr dolnoprzepustowy

Idealny filtr dolnoprzepustowy powinien po stronie częstotliwości być linią prostą o wartości dodatniej od f=0 do częstotliwości odcięcia i powyżej nie przyjmować wartość zero. Ponieważ w transformacie Fouriera odwrotnej chcąc po stronie czasu uzyskać wartości rzeczywiste prostokąt ten powinien być symetryczny względem f=0.

Najpierw zdefiniujmy sobie tzw. funkcję prostokątną:

$\text{rect}(t) = \begin{cases}0 & \text{if } \; |t| \;> \;^1/_2 \\ ^1/_2 & \text{if }\; |t| \;= \;^1/_2 \\1 & \text{if }\; |t| \;<\; ^1/_2 \\\end{cases}$

W naszym wypadku przyjmując za częstotliwość odcięcia B filtr po stronie częstotliwości powinien mieć postać $\displaystyle rect(\frac{f}{2B})$.

Jego transformatą po stronie czasu jest:

$\displaystyle f(t)=\mathcal{F}^{-1}\{ rect(\frac{f}{2B})\} = 2Bsinc(2Bt)$

Funkcja generująca taki filtr:

function [ k, t ] = sinc_gen( sampling_rate, cutfq, samples )

  T = 1 / sampling_rate; 

  if (nargin == 2)
    samples = 160 / (cutfq*2 / sampling_rate); % experimental
  end

  tout = T * (1:samples); 
  tout = tout - T*(samples+1)/2;
  B = 2*cutfq;
  k = sinc(B*tout);
  k = k / trapz(k);

  if (nargout == 2)
      t = tout;
  end

end

Jeśli nie podamy ilości próbek jest ona ustalana na najlepszą (według mnie) wartość, która przy filtrowaniu nie powoduje błędów.

Zmienna sampling_rate to częstotliwość próbkowania sygnału. sampling_rate/2 to największa częstotliwość w sygnale. Czyli np: jeśli sampling_rate=44100, to $f_{max}=22050$. Na samym końcu impuls jest normalizowany, tak by jego pole powierzchni było równe jeden, dzięki czemu splot z sygnałem zachowuje jego energię.

Zauważmy, że dla takiej samej wartości samples jeśli stosunek sampling_rate do cutoff_fq FIR będzie taki sam to odpowiedź impulsowa będzie taka sama (w sensie wartości w tablicy).

Przykładowy filtr sinc:

[~, w3] = sinc_gen(44100, 4000, 200);
tx = linspace(0, 200/44100, 200);
plot(tx, w3);
title('sampling rate=44100Hz, cutoff fq=4kHz, samples=200');

Filtr jest symetryczny i zmierza do zera jak czaszmierza do nieskończoności. Dla parzystej liczby sampling_rate dwie środkowe próbki mają takie same wartości. Dla nieparzystej największą wartość ma środkowa próbka (ten jest bardziej pożądany, gdyż nie przesuwa sygnału podczas splotu).

Wraz ze zwiększaniem sampling rate i utrzymywaniem wysokiej cutoff_fq kształt filtru będzie coraz bardziej przypominać impuls.

W dziedzinie częstotliwości filtr taki wygląda tak:

f3 = abs(fftshift(fft(w3)));
fx = linspace(-44100/2, 44100/2, 200);
plot(fx, f3);
title('sampling rate=44100Hz, cutoff fq=4kHz, samples=200');

Oscylacje na brzegach są spowodowane tym, że nas filtr jest gwałtownie się urywa na brzegach. Aby tego uniknąć na filtry nakłada się okno, jest ich wiele rodzajów, i nie zamierzam wnikać, który się do czego się nadaje.

w4 = w3 .* blackman(200)';
f3 = abs(fftshift(fft(w4)));
fx = linspace(-44100/2, 44100/2, 200);
plot(fx, f3);
title('sampling rate=44100Hz, cutoff fq=4kHz, samples=200');

Moja wyspecjalizowana funkcja generująca sinca:

function [ k, t ] = firwin( fs, fcut, len, prevent_shift )

k = sinc_gen(fs, fcut);

if (sum(size(len) - [1 1]) == 0)
    L = len;
else
    L = size(len, 2);
end

if (size(k, 2) > L)
    k = sinc_gen(fs, fcut, L);
end;

p = true;
if (nargin == 4)
    p = prevent_shift;
end

% prevent time shifting
if (p)
    if (mod(size(k, 2), 2) == 0)
        k = sinc_gen(fs, fcut, size(k, 2) - 1);
    end
end

% apply window
w = kaiser(size(k, 2), 9);
k = k .* w';

% normalize
k = k / trapz(k);

if (nargout == 2)
    t = linspace(0, (size(k, 2)-1)/fs, size(k, 2));
end

end

Korzysta ona z ustawionego na stałe okna (Kaiser, 9). Można także zapobiegać generacji sinca o parzystej liczbie próbek, splot z takim sincem w czasie przesuwa nam rezultat o jeden.

Porównanie do funkcji wbudowanej w Matlaba:

spectrum_cut = 50000;
Fs = 350000;
len = 100;

k1 = fir1(len-1, 2 * spectrum_cut / Fs, kaiser(len, 9));
k1 = k1 / trapz(k1);

k2 = firwin(Fs, spectrum_cut, len, false);

hold on;
plot(k1, '*');
plot(k2, 'r');

err = error_calc(k1, k2);

Funkcja pomocnicza:

function [ r ] = error_calc( v1, v2 )

res = sum(abs((v1 - v2) / v2) / size(v1, 2));

if (nargout == 0)
    display(sprintf('error: %f', res));
else
    r = res;
end

end

Splot

Poniżej podane są funkcje w Matlabie realizujące splot, działające jak ich odpowiedniki w Matlabie:

function [ fout ] = conv_full( fin, k )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  fout = zeros(1,s1+s2-1);
  
  for m=1:size(fout, 2)
    fin1 = max(m-s2+1, 1);
    fin2 = min(m, s1);
    k1 = max(m-s1+1, 1);
    k2 = min(m, s2);
    fout(m) = sum(fin(fin1:fin2) .* k(k2:-1:k1));
  end

end

function [ fout ] = conv_same( fin, k )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  L = fix(size(k,2)/2)+1;
  H = L + size(fin, 2) - 1;
  
  fout = zeros(1, H-L+1);
  
  for m=L:H
    fin1 = max(m-s2+1, 1);
    fin2 = min(m, s1);
    k1 = max(m-s1+1, 1);
    k2 = min(m, s2);
    fout(m-L+1) = sum(fin(fin1:fin2) .* k(k2:-1:k1));
  end

end

Funkcja conv_same przydaje się w przetwarzaniu sygnałów, wtedy z reguły u traktujemy jako sygnał zaś v jako odwrotną transformatę Fouriera jakiegoś filtru, którego chcemy użyć. Splotowi u i v w czasie odpowiada mnożenie po stronie częstotliwości. Wynik conv_same będzie takiego samego rozmiaru jak u.

Funkcja v może być także odpowiedzią impulsową funkcji przejścia H(s) dowolnego układu. Znając H(s) liczymy jaka jest odpowiedź układu na impuls (1 po stronie s), czyli jest to transformata odwrotna 1*H(s). Mówimy wtedy o odpowiedzi impulsowej filtru. Po stronie czasu mamy układ o funkcji przejścia H(s), podajemy na jego wejście impuls (po stronie częstotliwości widmo impulsu jest ciągłe od minus nieskończoności do plus nieskończoności i o stałej amplitudzie). To co otrzymujemy na wyjściu to odpowiedź w czasie przefiltrowanych wszystkich możliwych częstotliwości. Transformatą Fouriera tej odpowiedzi wyznacza nam charakterystykę częstotliwością filtru.

I jeszcze przykład funkcji realizującej splot dla danego punktu:

function [ r ] = conv_point( fin, k, point )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  if (point < 1)
    error('bad point index');
  end
  
  if (point > size(fin, 2))
    error('bad point index');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  L = fix(size(k,2)/2)+1;
  H = L + size(fin, 2) - 1;
  
  m = point + L - 1;
  
  fin1 = max(m-s2+1, 1);
  fin2 = min(m, s1);
  k1 = max(m-s1+1, 1);
  k2 = min(m, s2);
  r = sum(fin(fin1:fin2) .* k(k2:-1:k1));

end

Porównanie wersji same i full:

fsin1 = 500; % czestotliwosc pierwsze sinusoidy
fsin1_phase = pi/2; %-pi/4.14; % przesuniecie fazowe pierwszej sinusoidy
fsin2 = 2000; % czestotliwosc drugiej sinusoidy
fsin2_phase = -pi/4; %pi/2.11; % przesuniecie fazowe drugiej sinusoidy
fs = 24000; % czestotliwosc probkowania

% funkcja wejsciowa.
func_in = @(ti) sin(2*pi*fsin1*ti+fsin1_phase) + ...
sin(2*pi*fsin2*ti+fsin2_phase);

fcut = 1000; % czestotliwosc odciecia
tend = 3/fsin1; % dlugosc funkcji
sinc_samples = 200;

% przygotuj wektor czasu.
ts = 1/fs; % odstep czasu pomiedzy probkowaniami
tx = (0:ts:tend);

if (mod(size(tx, 2), 2) == 1)
tx = tx(1:1:end-1);
end

% wejscie
fin = func_in(tx)/15;

% filtr
w1 = sinc_gen(fs, fcut, sinc_samples);

fout1 = conv(fin, w1, 'same');
fout2 = conv(fin, w1, 'full');

figure(1);
hold on;
plot(fin, 'r');
plot(w1, 'g');
plot(fout1, 'b');
plot(fout2, 'c');
legend('wejscie', 'filtr', 'conv full', 'conv same');

Przykład:

fs = 30000; % czestotliwosc probkowania
fcut = 1000; % czestotliwosc odciecia
fsin1 = 500; % czestotliwosc pierwsze sinusoidy
fsin1_phase = pi/2; %-pi/4.14; % przesuniecie fazowe pierwszej sinusoidy
fsin2 = 2000; % czestotliwosc drugiej sinusoidy
fsin2_phase = -pi/4; %pi/2.11; % przesuniecie fazowe drugiej sinusoidy
tend = 1/fsin1*3.35; % dlugosc funkcji
sinc_samples = 1000;

% funkcja wejsciowa.
func_in = @(ti) sin(2*pi*fsin1*ti+fsin1_phase) + ...
sin(2*pi*fsin2*ti+fsin2_phase);

% funkcja wyjsciowa - idealne odciecie.
func_out = @(ti) sin(2*pi*fsin1*ti+fsin1_phase);

% przygotuj wektor czasu.
ts = 1/fs; % odstep czasu pomiedzy probkowaniami
tx = (0:ts:tend);

% oblicz wartosci obu funkcji w czasie.
fin = func_in(tx);
fout = func_out(tx);

%[~, fin] = sinc_gen(fs, (fsin1+fsin2)*1.12, size(tx, 2));
%[~, fout] = sinc_gen(fs, fsin1, size(tx, 2));

% przygotuj sinca
[~, w1] = sinc_gen(fs, fcut, sinc_samples);

% oba okna w czestotliwosci.
fw1 = fft(w1);
fwx = linspace(-fs/2, fs/2, size(fw1, 2));

% funkcja wejsciowa w czestotliwosci.
ffin = fftshift(abs(fft(fin)));
ffinx = linspace(-fs/2, fs/2, size(ffin, 2));

% przyciecie poprzez splot
fout2 = conv(fin, w1, 'same');

figure(1);
hold on;
plot(fwx, fftshift(abs(fw1)), 'r');
legend('sinc');
title('filtr dolnoprzepustowy w przestrzeni czestotliwosci');

figure(2);
hold on;
plot(fin, 'c');
plot(fout);
plot(fout2, 'r');
title('porowanie w czasie');
legend('wejscie', 'idealne wyjscie', 'przyciecie przez splot');

figure(3);
hold on;
plot(ffinx, fftshift(abs(fft(fin))), 'c');
plot(ffinx, fftshift(abs(fft(fout))));
plot(ffinx, fftshift(abs(fft(fout2))), 'r');
title('porowanie w czestotliwosci');
legend('wejscie', 'idealne wyjscie', 'przyciecie przez splot');

Po odkomentowaniu linii zastępujących nasze funkcje sincami, wyniki przedstawiają się tak:

Test naszych funkcji:

% przetestujmy poprawnosc dzialania naszej funkcji conv_cut.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'same');
x12 = conv(a, b, 'same');
x13 = conv(a, c, 'same');
x14 = conv(a, d, 'same');

x21 = conv(b, a, 'same');
x22 = conv(b, b, 'same');
x23 = conv(b, c, 'same');
x24 = conv(b, d, 'same');

x31 = conv(c, a, 'same');
x32 = conv(c, b, 'same');
x33 = conv(c, c, 'same');
x34 = conv(c, d, 'same');

x41 = conv(d, a, 'same');
x42 = conv(d, b, 'same');
x43 = conv(d, c, 'same');
x44 = conv(d, d, 'same');

y11 = conv_same(a, a);
y12 = conv_same(a, b);
y13 = conv_same(a, c);
y14 = conv_same(a, d);

y21 = conv_same(b, a);
y22 = conv_same(b, b);
y23 = conv_same(b, c);
y24 = conv_same(b, d);

y31 = conv_same(c, a);
y32 = conv_same(c, b);
y33 = conv_same(c, c);
y34 = conv_same(c, d);

y41 = conv_same(d, a);
y42 = conv_same(d, b);
y43 = conv_same(d, c);
y44 = conv_same(d, d);

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

% przetestujmy poprawnosc dzialania naszej funkcji conv_full.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'full');
x12 = conv(a, b, 'full');
x13 = conv(a, c, 'full');
x14 = conv(a, d, 'full');

x21 = conv(b, a, 'full');
x22 = conv(b, b, 'full');
x23 = conv(b, c, 'full');
x24 = conv(b, d, 'full');

x31 = conv(c, a, 'full');
x32 = conv(c, b, 'full');
x33 = conv(c, c, 'full');
x34 = conv(c, d, 'full');

x41 = conv(d, a, 'full');
x42 = conv(d, b, 'full');
x43 = conv(d, c, 'full');
x44 = conv(d, d, 'full');

y11 = conv_full(a, a);
y12 = conv_full(a, b);
y13 = conv_full(a, c);
y14 = conv_full(a, d);

y21 = conv_full(b, a);
y22 = conv_full(b, b);
y23 = conv_full(b, c);
y24 = conv_full(b, d);

y31 = conv_full(c, a);
y32 = conv_full(c, b);
y33 = conv_full(c, c);
y34 = conv_full(c, d);

y41 = conv_full(d, a);
y42 = conv_full(d, b);
y43 = conv_full(d, c);
y44 = conv_full(d, d);

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

%% przetestujmy poprawnosc dzialania naszej funkcji conv_point.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'same');
x12 = conv(a, b, 'same');
x13 = conv(a, c, 'same');
x14 = conv(a, d, 'same');

x21 = conv(b, a, 'same');
x22 = conv(b, b, 'same');
x23 = conv(b, c, 'same');
x24 = conv(b, d, 'same');

x31 = conv(c, a, 'same');
x32 = conv(c, b, 'same');
x33 = conv(c, c, 'same');
x34 = conv(c, d, 'same');

x41 = conv(d, a, 'same');
x42 = conv(d, b, 'same');
x43 = conv(d, c, 'same');
x44 = conv(d, d, 'same');

y11 = zeros(1, size(a, 2));
y12 = zeros(1, size(a, 2));
y13 = zeros(1, size(a, 2));
y14 = zeros(1, size(a, 2));
for i=1:1:size(a, 2)
y11(i) = conv_point(a, a, i);
y12(i) = conv_point(a, b, i);
y13(i) = conv_point(a, c, i);
y14(i) = conv_point(a, d, i);
end

y21 = zeros(1, size(b, 2));
y22 = zeros(1, size(b, 2));
y23 = zeros(1, size(b, 2));
y24 = zeros(1, size(b, 2));
for i=1:1:size(b, 2)
y21(i) = conv_point(b, a, i);
y22(i) = conv_point(b, b, i);
y23(i) = conv_point(b, c, i);
y24(i) = conv_point(b, d, i);
end

y31 = zeros(1, size(c, 2));
y32 = zeros(1, size(c, 2));
y33 = zeros(1, size(c, 2));
y34 = zeros(1, size(c, 2));
for i=1:1:size(c, 2)
y31(i) = conv_point(c, a, i);
y32(i) = conv_point(c, b, i);
y33(i) = conv_point(c, c, i);
y34(i) = conv_point(c, d, i);
end

y41 = zeros(1, size(d, 2));
y42 = zeros(1, size(d, 2));
y43 = zeros(1, size(d, 2));
y44 = zeros(1, size(d, 2));
for i=1:1:size(d, 2)
y41(i) = conv_point(d, a, i);
y42(i) = conv_point(d, b, i);
y43(i) = conv_point(d, c, i);
y44(i) = conv_point(d, d, i);
end

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

Funkcja pomocnicza:

function [ k, t ] = sinc_gen( sampling_rate, cutfq, samples )

T = 1 / sampling_rate;

if (nargin == 2)
samples = 160 / (cutfq*2 / sampling_rate); % experimental
end

tout = T * (1:samples);
tout = tout - T*(samples+1)/2;
B = 2*cutfq;
k = sinc(B*tout);
k = k / trapz(k);

if (nargout == 2)
t = tout;
end

end

Impuls, Splot, rodzaje filtrów, odpowiedź impulsowa

W poprzednich postach podałem przykład funkcji sinc i jej transformacji. W miarę jak będziemy brać coraz większy zakres czasu prostokąt po stronie częstotliwości będzie mieć coraz mniejsze oscylacje w punkcie skoku. W miarę jak nasza funkcja sinc będzie coraz węższa nasz prostokąt będzie coraz szerszy. Nasza funkcja sinc będzie coraz bardziej przypominać impuls. Innymi słowy transformatą impulsu jest nieskończone widmo częstotliwości o takiej samej amplitudzie. Impuls jest paczką sinusoid o wszystkich częstotliwościach i takiej samej amplitudzie.

Splot funkcji definiujemy jako:

$f(t)*h(t) = \int_{-\infty}^{\infty} f(\tau)h(t-\tau)d\tau$

Istnieje także jego dyskretny odpowiednik. W transformacie Fouriera splot ma bardzo ważną cechę. Mianowicie mnożeniu funkcji po stronie częstotliwości odpowiada splot po stronie czasu i na odwrót.

Co się stanie jeśli na jakiś filtr, czyli obiekt który odpowiednio filtruje częstotliwości sygnały wejściowego podamy impuls. Impuls to nieskończone widmo częstotliwości. Odpowiedz którą otrzymamy to suma odfiltrowanych sinusoid. Taka odpowiedź daje nam jednoznaczną informację o charakterystyce częstotliwościowej filtru.

W zależności od cech odpowiedzi możemy różnie klasyfikować filtry. Ze względu na energię sygnału wyjściowego: aktywne (dodają energii) i bierne (nie dodają energii, inaczej pasywne). Ze względu na czas po jakim następuje odpowiedź na sygnał wejściowy, czy jest on zawsze taki sam, czy od czegoś zależy np. od sygnału wcześniejszego. W szczególności mówimy, że filtr jest niezmienny w czasie (time-invariant) jeśli na każdą pojedyńczą sinusoidę odpowiada on w taki sam sposób i po takim samym czasie. Inną cechą filtrów jest ich linearność, bądź i nie. Filtr jest linearny jeśli dowolne rozłożenie sygnału wejściowego, w szczególności na pojedyncze sinusoidy i ich przefiltrowanie, a następnie złożenie sygnału wyjściowego równe jest odpowiedzi filtru na złożony sygnał wejściowy. Inaczej filtr filtruje każdą częstotliwość niezależnie od innych. Innym podziałem jest FIR (finite impulse response) i IIR (infinite impulse response). Odpowiedź filtru FIR zanika w skończonym czasie, zaś odpowiedź filtru IIR nie - asymptotycznie zmierza do zera z reguły oscylując. Filtr IIR możemy przyciąć do FIR zwłaszcza jeśli posługujemy się jego dyskretną postacią.

Wszystkie filtry o których będziemy mówić są pasywne, liniowe, time-invariant, FIR.

Robiąc z filtru IIR filtr FIR wprowadzamy błąd. Dobierając odpowiednio długość FIR możemy te zniekształcenia przenieść poniżej rozdzielczości sygnału. Lub powyżej częstotliwości poniżej których nie powinno ich być (jak w sygnale audio). Część błędu, czyli oscylacje na brzegach możemy zmniejszyć splatając odpowiedzieć impulsową z oknem. Funkcji okien jest bardzo wiele - po stronie częstotliwości mają one kształt okna, najprostszym jest RECT.

Jak wcześniej powiedzieliśmy mnożeniu po stronie częstotliwości odpowiada splot po stronie czasu. W tym przypadku splot funkcji i odpowiedzi impulsowej. Jest to bardzo ważna cecha pozwalająca nam bez użycia transformat filtrować sygnał.

Przesunięcie w czasie, a FFT

Przesunięcie w czasie nie wpływa na absolutne wartości amplitudy w FFT. Zmienia się tylko faza składowych sinusoid. Dla symetrycznego wykresu względem jego środka po stronie czasu część urojona amplitudy jest zerowa. Połówka symetrycznego wykresu będzie miała takie same absolutne amplitudy po stronie częstotliwości jak cały symetryczny wykres, różnica będzie w przesunięciach fazowych. Niesymetryczność po stronie FFT amplitud lub faz, powoduje, że po stronie czasu nasz wykres nie jest
rzeczywisty. Z uwagi na błędy zaokrąglenia rozsądnie jest jeśli spodziewamy się symetryczności dla bezpieczeństwa brać wartości bezwzględne
przy konwersji z czasu do częstotliwości, a przy konwersji odwrotnej część rzeczywistą wyniku.

Wpływ przesunięcia w czasie na FFT:

function [ k, t ] = sinc_gen( sampling_rate, cutfq, samples )

  T = 1 / sampling_rate; 

  if (nargin == 2)
    samples = 160 / (cutfq*2 / sampling_rate); % experimental
  end

  tout = T * (1:samples); 
  tout = tout - T*(samples+1)/2;
  B = 2*cutfq;
  k = sinc(B*tout);
  k = k / trapz(k);

  if (nargout == 2)
    t = tout;
  end

end

spectrum_cut = 20000;
fs = 3.5e6;

[k, t] = firwin(fs, spectrum_cut, 2047);

k1 = [zeros(1, 700), k];

K = fftshift(fft(k, 2047));
K1 = fftshift(fft(k1, 2747));
fx = linspace(-fs/2, fs/2, 2047);
fx1 = linspace(-fs/2, fs/2, 2747);

figure(1);
hold on;
plot(fx, abs(K), '.');
plot(fx1, abs(K1), 'r');
legend('org', 'shifted');

figure(2);
hold on;
plot(fx, unwrap(angle(K)), '.');
plot(fx1, unwrap(angle(K1)), 'r');
legend('org', 'shifted');

Niesymetryczność w czasie to niezerowa część urojona po stronie FFT:

figure(3);
hold on;
plot(fx, imag(K), '.');
plot(fx1, imag(K1), 'r');
legend('org', 'shifted');

Transformata Fouriera

Istnieje wiele transformat, które generalnie można powiedzieć przekształcają funkcję z jednej przestrzeni w inną. Z reguły istnieje także transformata odwrotna. Nie wszystkie funkcje da się także poddać takiej transformacji. Transformata Fouriera to takie przekształcenie, które możemy zinterpretować jako przekształcenie funkcji w dziedzinie czasie w funkcję w dziedzinie częstotliwości.

$f(\omega) = \int_{-\infty}^{\infty} f(x)\ e^{- 2\pi i x t}dt$

Transformata odwrotna:

$f(t) = \int_{-\infty}^{\infty} f(\omega)\ e^{2 \pi i x \omega}d\omega$

W wyniku transformacji funkcji o wartościach rzeczywistych zmiennej rzeczywistej transformata Fouriera jest symetryczna względem zera. Tylko taka symetria po stronie częstotliwości gwarantuje nam otrzymanie funkcji o wartościach rzeczywistych podczas transformacji odwrotnej. Ma to szczególne znaczenie podczas dyskretnej transformaty Fouriera kiedy niedokładności w symetrii są częste.

Po stronie częstotliwości widmo sygnału interpretujemy jako sumę sinusoid o danych częstotliwościach, przesunięciu fazowym i amplitudzie. Ponieważ we wzorach mamy całkę, a nie sumę, musiała by to być suma nieskończona. Przechodząc z całek na sumy przechodzimy na dyskretną transformatę Fouriera. Wynikiem transformaty jest funkcja liczby zespolonej - jej moduł to amplituda dla danej częstotliwości, jej przesunięcie fazowe to przesunięcie fazowe sinusoidy o danej częstotliwości.

Ponieważ z reguły jesteśmy zainteresowaniu amplitudą widma jej wykresy najczęściej podaje się po stronie częstotliwości. Przesunięcia fazowe z reguły pomija się. Zauważmy, że dla różnie poprzesuwanych składowych sinusoid widmo amplitud jest takie same. I drugiej strony przesuniecie funkcji w czasie nie zmienia jej widma. Właściwość ta jest wykorzystywana podczas kompresji obrazów i filmów. Szczególnym przypadkiem dobrania przesunięć fazowych jest funkcja o minimalnym przesunięciu fazowym (jest to taka funkcja w której w dowolnym czasie $<0,t>$ energia sygnału jest maksymalna).

W wyniku tego, że dyskretna transformata Fouriera operuje na dyskretnych wartościach wprowadza ona błędy. Wspomnianą wcześniej asymetrię funkcji rzeczywistej i nierzeczywistą postać jej transformacji odwrotnej. Innym rodzajem błędu jest to, że pojedyncza sinusoida może być reprezentowana poprzez rozmyty pik po stronie częstotliwości. Dlatego po retransformatą warto wziąć jej część rzeczystą (nie moduł).

W dalszej części będziemy mówili tylko o transformacji dyskretnej.

Oto przykład transformaty funkcji sinc:

Teoretycznie jej transformatą powinien być prostokąt. W rzeczywistości ponieważ wartości są dyskretne, a funkcja nie jest nieskończona, transformata nie jest idealnym prostokątem i wykazuje oscylacje na brzegach skoku. Są one tym mniejsze im nasz sinc jest dłuższy.

Dla transformaty odwrotnej: sinc po stronie częstotliwości będzie prostokątem po stronie czasu z takimi samymi oscylacjami na brzegach. Ponieważ nasz sinc jest skończony w czasie brakuje nam widma dla wyższych częstotliwości i one uniemożliwiają nam zreprodukowanie skoku i zlikwidowanie oscylacji - można powiedzieć ich wygaszenie.

Przykład jak może wyglądać transformacja złożenia dwóch sinusoid o f=650 i f=1500. Oczywiście w zależności jak dobierzemy krok próbkowania naszych ciągłych funkcji, ilości próbek czasu i ilości próbek częstotliwości możemy otrzymać wykres bardziej lub mniej przypominający piki.

Warto podkreślić, że dyskretna transformata jest bezwymiarowa podobnie jak wartości wejściowych próbek. Chcąc otrzymać wartości częstotliwości na osi musimy je policzyć.

Skończony sygnał w czasie nieokresowy wymaga do opisu nieskończonego widma częstotliwości. Podobnie jest z nieciągłościami w czasie. I drugą stronę: skończone widmo nie powoli nam na opisanie skończonego w czasie sygnału. Te dwa ograniczenia silnie wiążą się z DFFT, której liczba próbek jest skończona.