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