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.