2011-12-22

Filtr Lanczos

Filtr Lanczos jest zbudowany na bazie funkcji sinc, która jest odpowiedzią impulsową idealnego filtru dolnoprzepustowego. Dowolny sygnał możemy odfiltrować takim filtrem poprzez splot.

Idealny filtr dolnoprzepustowy jest nieciągły w przestrzeni częstotliwości i rozciąga się nieskończenie w czasie. W grafice komputerowej, w filtrowaniu próbek obrazu podczas jego skalowania lub resamplingu musimy się posłużyć przyciętą w przestrzeni wersją sinc. Dla likwidacji oscylacji w jego widmie biorących się z nieciągłości powstałych w wyniku jego przycięcia musi zostać wymnożony z oknem. Jako okna używa się tego stworzonego przez niejakiego pana Lanczos - stąd nazwa filtru. Dlaczego takim oknem, a nie innym nie wiem. Okazało się także w trakcie eksperymentów, że najlepsze wyniki daje okno o wielkości 2 i 3.

Jedno o czym musimy pamiętać, że tworzenie filtrów to w pewnej mierze połączenie sztuki z nauką i dodatkowo to, że każdy ma inne gusta i inną percepcję. Filtr Lanczos jest prawie idealnym filtrem dolnoprzepustowym. Jednak wcale nie znaczy to najlepszym. Mówiąc o filtrowaniu obrazu chodzi nam o odcięcie wysokich częstotliwości, ale należy pamiętać, że oko nie postrzega widma (w przestrzeni, ale także samej fali elektromagnetycznej), poza tym receptory w oku nie są rozłożone równomiernie by eliminować możliwość aliasingu, sam mózg dokłada swoje w usuwaniu z postrzeganego obrazu niechcianych zakłóceń. Każdy ma inaczej rozłożone czopki w oku i inaczej utkane neurony i inny bagaż doświadczeń. Czopki u każdego reagują na tą samą częstotliwość fali trochę inaczej. Każdy ma inną optykę oka, inne zmętnienie soczewki. I pewnie można by tak wymieniać.

Oba filtry o szerokościach 2 i 3 mają dla pewnym wartości ujemne wartości, stąd istnieje możliwość, że wynikowa wartość splotu będzie ujemna. Czyli należy przewidzieć możliwość radzenia sobie z ujemnymi wartościami koloru. Przy samplowaniu Jitter dla niektórych metod resamplowania możemy otrzymać wartość koloru większą niż jeden.

Niezależnie od tego jak wygląda widmo filtru, jeśli filtr ma ujemne wybrzuszenia oznacza to, że będzie on miał właściwości wyostrzające.

Funkcja okna ma postać:

$\displaystyle W(x) = \begin{cases} \frac{sin \displaystyle\frac{x \pi}{\tau}}{\displaystyle\frac{x \pi}{\tau}} & \text{dla} \;\;\; -\tau \le x \le \tau \\ \\\\\\\\\\\\\\ \qquad 0 & \text{dla} \;\qquad |x|>\tau \end{cases}$

, gdzie $\tau$ to szerokość okna.

Po złożeniu z funkcją sinc otrzymujemy filtr:

$F(x) = \displaystyle \begin{cases} \frac{sin x \cdot sin \displaystyle\frac{x \pi}{\tau}}{x \cdot \displaystyle\frac{x \pi}{\tau}} & \text{dla } \;\;\; -\tau \le x \le \tau \\ \\\\\\\\\\\\\\\\\\\\ \qquad\qquad 0 & \text{dla } \;\qquad |x|>\tau \end{cases}$

Dla zera $\lim_{x \rightarrow 0}{\frac{sinx}{x}}\rightarrow 1$.

Dla $\tau=3$ wykresy filtru, okna, samej funkcji sinc mają postać:


Pole powierzchni filtra to 1, czyli dokonując z nim splotu nie zmieniamy energii sygnału.

Porównanie widma funkcji przyciętej funkcji sinc bez i z oknem Lanczos:


Filtr dwuwymiarowy $L(x,y)=L(x)*L(y)$ ma wygląd:


Jak widać nie powstaje on przez obrót filtru 1D. Objętość pod powierzchnią filtra 2D to jeden.

Przykładowy kod w Matlabie użytego do wygenerowania wykresów:

function [] = plot_lanczos()

    x = linspace(-4, 4, 1000);

    [y, ysinc, ywin] = lanczos(x, 3);

    trapz(x, y)

    figure(fig_index);
    fig_index = fig_index + 1;
    hold on;
    plot(x, ywin, 'r o');
    plot(x, y, 'b.');
    plot(x, ysinc, 'g');
    title('lanczos - spatial');
    legend('window', 'lanczos', 'sinc');

    x = linspace(-3, 3, 100);
    [y, ysinc, ~] = lanczos(x, 3);

    Y = fft(y, 1000);
    YSINC = fft(ysinc, 1000);

    figure(fig_index);
    fig_index = fig_index + 1;
    hold on;
    title('lanczos - freq');
    plot(fftshift(abs(Y)), 'r');
    plot(fftshift(abs(YSINC)));

end

function [y, ysinc, ywin] = lanczos(tau)

    y = zeros(1, size(x, 2));
    ysinc = zeros(1, size(x, 2));
    ywin = zeros(1, size(x, 2));

    for i=1:size(y, 2)
        xx = abs(x(i));
        if (xx == 0)
            ysinc(i) = 1;
            ywin(i) = 1;
        else
            xx = xx * pi;
            xxx = xx / tau;
            ysinc(i) = sin(xx) / xx;
            if (xx/pi <= tau)
                ywin(i) = sin(xxx) / xxx;
            end  
        end

        y(i) = ysinc(i) * ywin(i);
    end
end

function [] = plot_lanczos_3d()

    N = 51;
    tau = 3;

    x = linspace(-tau, tau, N);
    y = linspace(-tau, tau, N);

    lx = linspace(-tau, tau, N);
    ly = linspace(-tau, tau, N);

    for i=1:N 
        lx(i) = lanczos(lx(i), tau);
        ly(i) = lanczos(ly(i), tau);
    end

    [lx, ly] = meshgrid(lx, ly);

    lz = lx .* ly;

    surface(x, y, lz);

    trapz(y, trapz(x, lz))

    function [y] = lanczos(x, tau)

        xx = abs(x);
        if (xx == 0)
            y = 1;
        else
            xx = xx * pi;
            xxx = xx / tau;
            ysinc = sin(xx) / xx;
            ywin = 0;
            if (xx/pi <= tau)
                ywin = sin(xxx) / xxx;
            end  
            y = ysinc * ywin;
        end
    end
end

Przykładowy kod klasy w C# implementujące filtr:

public class LanczosFilter : Filter
{
    public double Tau = 3;

    public override double Ray
    {
        get
        {
            return Tau;
        }
    }

    public override double Evaluate(double a_value)
    {
        double v = Math.Abs(a_value);

        if (v > Tau)
            return 0;

        if (v < Constants.DOUBLE_PRECISION)
            return 1;

        v = v * Constants.PI;
        double p3 = v / Tau;
        return Math.Sin(v) * Math.Sin(p3) / (v * p3);
    }
}

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.

2011-12-19

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

Resampler LocalGrid

Budujemy dwuwymiarową (może też być jednowymiarowa) siatkę wag filtra. Środek filtru jest pozycji (0,0). Stąd dla filtru o promieniu 3 siatka musi mieć wymiary 7x7. Odległości z których pobieramy wartości filtra to -3, -2, -1, 0, 1, 2, 3. Zakładając, że wartości brzegowe w naszej siatce to zera możemy wymiary siatki zmniejszyć o 2.

Taką siatkę centrujemy w punkcie padania promienia i dla każdego piksela, na który nakłada się siatka dodajemy wartość koloru promienia ważoną wartością z siatki.

W przypadku supersamlingu wymiary siatki mnożymy przez pierwiastek z ilości promieni przewidywanych na piksel.

Bardzo ważne jest by uwzględnić wartości zwracane przez Sampler, tzn. zamiast ... 13.5, 14.5, 15.5, ... możemy dostać ... 13.49, 14.5, 15.51, ... (oczywiście z dużo mniejszym błędem). Tego rodzaju błędów nie da się uniknąć - zaokrąglanie liczb binarnych do postaci dziesiętnej.

W samym Resamplerze dodałem właściwość Subreslution, która dodatkowo zwiększa wymiary siatki. Jej zadaniem jest zwiększenie dokładności resamplowania, tak, by metoda ta w miarę zwiększania właściwości coraz bardziej przypominała resampler Exact, a coraz mniej Resizer. By kod mógł współdziałać zarówno z filtrem Box o promieniu 0.5 jak i z innymi filtrami o promieniach całkowitych Subreslution musi być nieparzyste. Inaczej wagi albo jednego albo drugiego filtru będą oscylować na granicy promienia w zależności od błędów numerycznych Samplera.

W analizie kodu najważniejsze jest dobre zrozumienie wszystkich $\pm\frac{1}{2}$, $\pm 1$, Floor, Ceiling, Round.

Kod:

public struct PixelInterpolator
{
    private ColorFloat m_numerator;
    private double m_denumerator;

    /// <summary>
    /// Note: a_ray_color in SRGB space.
    /// </summary>
    /// <param name="a_ray_color"></param>
    /// <param name="a_filter_weight"></param>
    public void Collect(ColorFloat a_ray_color, double a_filter_weight)
    {
        m_numerator += a_ray_color * a_filter_weight;
        m_denumerator += a_filter_weight;
    }

    /// <summary>
    /// Note: color in SRGB space.
    /// </summary>
    public ColorFloat Color
    {
        get
        {
            if (m_denumerator == 0)
                return ColorFloat.Black;

            return m_numerator / m_denumerator;
        }
    }

    public override string ToString()
    {
        return String.Format("num: {0}, denum: {1}", m_numerator, m_denumerator);
    }
}

public abstract class InterpolateResampler : FilterResampler
{
    protected PixelInterpolator[,] m_pixels;
    protected OverlayCorrector m_overlay_corrector;

    protected override ColorFloat ResamplePixel(int a_x, int a_y)
    {
        return Gamma.SRGBToLinear(m_pixels[a_x, a_y].Color);
    }

    internal override void RenderStart(RenderStartPhase a_phase)
    {
        if (a_phase == RenderStartPhase.PrepareObjectToRender)
        {
            m_pixels = new PixelInterpolator[Film.Width, Film.Height];
            m_overlay_corrector = OverlayCorrector.Create(
                OverlayMethod.Mirror, Film.Width, Film.Height, 
                Filter.Ray.Ceiling());
        }
    }

}

public class LocalGridResampler : InterpolateResampler
{
    private double[,] m_weights;

    [YAXNode]
    public int Subresolution = 3;

    private double m_scale;

    internal override void RenderStart(RenderStartPhase a_phase)
    {
        base.RenderStart(a_phase);

        if (a_phase == RenderStartPhase.PrepareObjectToRender)
        {
            var sampler = Film.Sampler as UniformSampler;

            if (sampler == null)
                throw new InvalidOperationException();
                    
            if (Subresolution % 2 == 0)
                throw new InvalidOperationException();

            m_scale = Subresolution * sampler.Subresolution;
            PrepareWeights();
        }
    }

    private void PrepareWeights()
    {
        List<List<PixelWeight>> horz_weights, vert_weights;
        ResizerResampler.PrecalculateWeights((int)m_scale, (int)m_scale, 1, 1, 
            out horz_weights, out vert_weights, Filter, false);
        double[] weights = horz_weights[0].Select(pw => pw.Weight).ToArray();

        Debug.Assert(weights.Length % 2 == 1);

        m_weights = new double[weights.Length, weights.Length];

        for (int x = 0; x < weights.Length; x++)
        {
            double wx = weights[x];

            for (int y = 0; y < weights.Length; y++)
            {
                double wy = weights[y];
                m_weights[x, y] = wx * wy;
            }
        }
    }

    private double GetWeight(double a_dx, double a_dy)
    {
        int x = (a_dx * m_scale).Round();
        x += (m_weights.GetLength(0) - 1) / 2;

        if (x < 0)
            return 0;
        else if (x >= m_weights.GetLength(0))
            return 0;

        int y = (a_dy * m_scale).Round();
        y += (m_weights.GetLength(1) - 1) / 2;

        if (y < 0)
            return 0;
        else if (y >= m_weights.GetLength(1))
            return 0;

        return m_weights[x, y];
    }

    internal override void SetRayColor(Vector2 a_pixel, ColorFloat a_color)
    {
        a_color = Gamma.LinearToSRGB(a_color);

        int start_x = (int)Math.Floor(a_pixel.X - Filter.Ray);
        int end_x = (int)Math.Ceiling(a_pixel.X + Filter.Ray) - 1;

        int start_y = (int)Math.Floor(a_pixel.Y - Filter.Ray);
        int end_y = (int)Math.Ceiling(a_pixel.Y + Filter.Ray) - 1;

        for (int y = start_y; y <= end_y; y++)
        {
            int cy = m_overlay_corrector.CorrectY(y);

            for (int x = start_x; x <= end_x; x++)
            {
                int cx = m_overlay_corrector.CorrectX(x);
                double w = GetWeight((x + 0.5) - a_pixel.X, (y + 0.5) - a_pixel.Y);
                m_pixels[cx, cy].Collect(a_color, w);
            }
        }

        #if DEBUG
        for (int y = start_y - 1; y <= end_y + 1; y++)
        {
            Debug.Assert(GetWeight((start_x + 0.5 - 1) - a_pixel.X, (y + 0.5) - 
                a_pixel.Y).IsAlmostEquals(0));
            Debug.Assert(GetWeight((end_x + 0.5 + 1) - a_pixel.X, (y + 0.5) - 
                a_pixel.Y).IsAlmostEquals(0));
        }
        for (int x = start_x; x <= end_x; x++)
        {
            Debug.Assert(GetWeight((x + 0.5) - a_pixel.X, (start_y + 0.5 - 1) - 
                a_pixel.Y).IsAlmostEquals(0));
            Debug.Assert(GetWeight((x + 0.5) - a_pixel.X, (end_y + 0.5 + 1) - 
                a_pixel.Y).IsAlmostEquals(0));
        }
        #endif
    }

    public override ResamplerType ResamplerType
    {
        get
        {
            return ResamplerType.LocalGrid;
        }
    }
}

Resampler Exact

Dla każdego promienia skanujemy piksele wokół (uwzględniamy promień filtra). Dla każdego piksela wagę liczymy dokładnie z filtra. Nie korzystamy z żadnych prekalkulowanych tabelek. Metoda ta powinna zapewniać największą dokładność.

W analizie kodu najważniejsze jest dobre zrozumienie wszystkich $\pm\frac{1}{2}$, $\pm 1$, Floor, Ceiling, Round.

Kod:
public struct ColorInterpolator
{
    private ColorFloat m_numerator;
    private double m_denumerator;

    public void Collect(ColorFloat a_ray_color, double a_filter_weight)
    {
        m_numerator += a_ray_color * a_filter_weight;
        m_denumerator += a_filter_weight;
    }

    public ColorFloat Color
    {
        get
        {
            if (m_denumerator == 0)
                return ColorFloat.Black;

            return m_numerator / m_denumerator;
        }
    }
}

public abstract class Resampler
{
    public Filter Filter = new BoxFilter();
    internal Film Film;

    public static Resampler Create(ResamplerType a_type)
    {
        if (a_type == ResamplerType.Exact)
            return new ExactResampler();
        else if (a_type == ResamplerType.LocalGrid)
            return new LocalGridResampler();
        else if (a_type == ResamplerType.Resizer)
            return new ResizerResampler();
        else
            throw new NotImplementedException();
    }

    public Rectangle GetAffectedRect(Rectangle a_src_rect)
    {
        int border = (int)Math.Ceiling(Filter.Ray) + 2;
        a_src_rect.Inflate(border, border);
        a_src_rect.Intersect(new Rectangle(0, 0, Film.Width, Film.Height));
        return a_src_rect;
    }

    public void Resample(Rectangle a_rect, ColorArrayBase a_dest)
    {
        for (int y = a_rect.Top; y < a_rect.Bottom; y++)
        {
            for (int x = a_rect.Left; x < a_rect.Right; x++)
            {
                a_dest.SetColor(x, y, ResamplePixel(x, y));
            }
        }
    }

    public abstract void RenderStart(RenderStartPhase a_phase);
    protected abstract ColorFloat ResamplePixel(int a_x, int a_y);
    public abstract void SetRayColor(Vector2 a_pixel, ColorFloat a_color);
}

public class ExactResampler : InterpolateResampler
{
    public override void SetRayColor(Vector2 a_pixel, ColorFloat a_color)
    {
        int start_x = (int)Math.Floor(a_pixel.X - Filter.Ray);
        int end_x = (int)Math.Ceiling(a_pixel.X + Filter.Ray) - 1;

        int start_y = (int)Math.Floor(a_pixel.Y - Filter.Ray);
        int end_y = (int)Math.Ceiling(a_pixel.Y + Filter.Ray) - 1;

        for (int y = start_y; y <= end_y; y++)
        {
            int cy = m_overlay_corrector.CorrectY(y);
            double wy = Filter.Evaluate((y + 0.5) - a_pixel.Y);

            for (int x = start_x; x <= end_x; x++)
            {
                int cx = m_overlay_corrector.CorrectX(x);
                double wx = Filter.Evaluate((x + 0.5) - a_pixel.X);

                m_rays[cx, cy].Collect(a_color, wy * wx);
            }
        }

        Debug.Assert(Filter.Evaluate(
            a_pixel.X - (start_x + 0.5 - 1)).IsAlmostEquals(0));
        Debug.Assert(Filter.Evaluate(
            a_pixel.X - (end_x + 0.5 + 1)).IsAlmostEquals(0));
        Debug.Assert(Filter.Evaluate(
            a_pixel.Y - (start_y + 0.5 - 1)).IsAlmostEquals(0));
        Debug.Assert(Filter.Evaluate(
            a_pixel.Y - (end_y + 0.5 + 1)).IsAlmostEquals(0));
    }
}

Resamplowanie Exact i LocalGrid - porównanie

Resampling LocalGrid jest nieznacznie szybszy. Dla filtru Lanczos który wymaga policzenia dwóch funkcji sinus i jest chyba najkosztowniejszym obliczeniowo filtrem jest to 10%.

Jeśli chodzi o jakość to już Resampler LocalGrid z minimalną Subresolution równą 2 jest prawie identyczny jak Exact. Dalsze zwiększanie Subresolution nie poprawia jakości generowanego obrazu.

2011-12-18

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');