Filtracja sygnałów o nieznanych parametrach zaszumionych szumem białym

Spotykamy się w wielu przypadkach z potrzebą przefiltrowania sygnału który został poddany „wpływowi otoczenia”. W przypadku sygnałów natury elektrycznej spotykamy takie wpływy jak indukcyjność, przydźwięki, rozpraszanie się fal, etc. Są to dla nas składowe niepożądane w sygnale mierzonym. Staramy się te wpływy otoczenia wyeliminować separując przewody, używając filtrów ferrytowych. Co prawda polepsza to nasz sygnał, ale nie jest w stanie wyeliminować wszystkich jego defektów. Mierzony sygnał też jest narażony na niepoprawną pracę przetwornika pomiarowego który w stanie ustalonym pokazuje nam wartość zmieniającą się o jeden bit.

Powszechnie stosowanym rozwiązaniem jest użycie filtru dolnoprzepustowego gdyż w większości przypadków zakłócenia sygnału mają częstotliwość wyższa niż składowe sygnału bazowego. Takie rozwiązanie jest bardzo fajne jeżeli znamy częstotliwości składowe oraz są one względnie niskie.

Poważny problem zaczyna się w momencie gdy składowe sygnału nie są znane. Istnieje wtedy możliwość wytłumienia składowej z której składa się sygnał bazowy. Należało by zaprojektować taki filtr który sprawdzał by się w wielu przypadkach w których nie są znane informacje o sygnale. Problem wydaje się nie do rozwiązania, ale to błędne założenie.

W niniejszym artykule przedstawiono porównanie pracy filtrów dolno-przepustowych, filtrów adaptacyjnych i optymalnego filtru Wienera.

Sygnał bazowy

Załóżmy, że nie znamy wszystkich informacji o sygnale. Wiemy za to jakie jest jego próbkowanie (20ksps). Sygnał bazowy trwa 1 sekundę. Ponadto wiemy, że sygnał składa się z dziesięciu składowych harmonicznych które występują w zakresie. Wielomian opisujący sygnał bazowy można wyrazić poprzez.

%%  Sygnal bazowy
    %   Probkowanie 20ksps
    F_s = 20000;
    %   Czas trwania sygnalu 2 s
    t_final = 1;
    %   Wektor czasu
    t = 0 : ( 1 / F_s ) : t_final - ( 1 / F_s );
    t=t';
    L = length( t );

    n1 = 0;
    n2 = 0.05;

    Amplitude = abs( randn(10,1) * 10 ) + 1;
    Frequency = abs( randn(10,1) * 1000 ) + 10;

    disp('Based Amplitude:')
    Amplitude
    disp('Based Frequency:')
    Frequency
    Signal = 0;
    for i = 1 : length( Amplitude )
        Signal = Signal + ...
            Amplitude(i) * sin( 2 * pi * t * Frequency(i) );
    end

 

Zakłócenia są szumem białym czyli o widmie płaskim . Ponadto aby sygnał nie był możliwy do wykorzystania, amplituda szumu jest zbliżona do amplitudy składowych harmonicznych.

Sygnał wykorzystywany do badania filtrów jest sumą sygnału bazowego i sygnału zakłócenia v1.

v = 10*randn(L,1);
ar = [1, 1/2];
v1 = filter(1,ar,v*5); 
Signal_2 = Signal + v1;

 

Tak powstały sygnał bez jakiegokolwiek opisu statystycznego jest nierozpoznawalny jedynie opis statystyczny uzyskany na przykład przy użyciu szybkiej transformaty Fouriera.

Projektowanie filtrów dolnoprzepustowych

Podstawowym podejściem do filtracji tego typu sygnału jest zaprojektowanie filtra dolnoprzepustowego o wysokiej częstotliwości odcięcia () i paśmie przepustowym (). Wszystkie filtry użyte do badania posiadają zafalowanie w paśmie przepustowym rzędu i tłumienie na poziomie .

Filtry użyte do badań: filtr Butterwortha, Chebyszewa I rodzaju, Chebyszewa II rodzaju, Equiripple oraz Eliptyczny.

 

d = fdesign.lowpass('Fp,Fst,Ap,Ast',F_pass,F_stop,R_pass,R_stop,F_s);

mc_make_filter(d, 'butter', Signal, Signal_2, t)
mc_make_filter(d, 'cheby1', Signal, Signal_2, t)
mc_make_filter(d, 'cheby2', Signal, Signal_2, t)
mc_make_filter(d, 'equiripple', Signal, Signal_2, t)
mc_make_filter(d, 'ellip', Signal, Signal_2, t)
mc_make_filter(d, 'ifir', Signal, Signal_2, t)
mc_make_filter(d, 'kaiserwin', Signal, Signal_2, t)

 

 

%%  mc_make_filter
%   Mariusz Czajkowski
%   21.12.2014
%   www.uengineering.net
%   m.czajkowski@uengineering.net

function [hd] = mc_make_filter( d, FilterType, base, x, t )
    %   Tworzenie filtru
    hd = design(d,FilterType);
    
%%  Wynik filtracji 
    y = filter(hd, x);
    
    figure()
    subplot(2,1,1)
    plot(t, base, '-',...
        t, x, ':',...
        t, y, '--')
    title('Results of filtration')
    xlabel('Time [s]')
    ylabel('Amplitude')
    legend('Base signal', 'Noised signal', 'Filtered signal')
    xlim([0 0.05])
    
    subplot(2,1,2)
    plot(t, base, '-',...
        t, x, ':',...
        t, y, '--')
    title('Results of filtration')
    xlabel('Time [s]')
    ylabel('Amplitude')
    legend('Base signal', 'Noised signal', 'Filtered signal')
    xlim([0.1 0.105])
    
    n = 2^12;    
    figure()
    subplot(3,1,1)
    fft_analys(base, t, n);
    title('FFT analisys of base signal')
    xlim([0 5000])
    
    subplot(3,1,2)
    fft_analys(x, t, n);
    title('FFT analisys of noised signal')
    xlim([0 5000])
    
    subplot(3,1,3)
    fft_analys(y, t, n);
    title('FFT analisys of filtered signal')
    xlim([0 5000])
end

 

Filtr Butterworth’a

Filtr Chebyszewa I rodzaju

Chebyszewa II rodzaju

Filtr Eliptyczny

Filtr Equiripple

Z racji nie znania częstotliwości składowych wyższe składowe zostały stłumione przez filtry. Pokazuje to wadę używania takiego rozwiązania.

Filtry adaptacyjne

Do zaprojektowania filtrów adaptacyjnych użyto wzoru szumu bazowego o innych parametrach niż szum . Różnicę pomiędzy szumami przedstawiono na rysunku:

Do projektowania filtrów adaptacyjnych użyto algortymów LMS I NLMS.

%%  Konstruowanie filtru adaptacyjnego
LL = 7;                              % Ilosc wspolczynnikow filtru adaptacyjnego
hlms = adaptfilt.lms(LL);            % Filtr adaptacyjny FIR z uzyciem LMS
hnlms = adaptfilt.nlms(LL);          % Znormalizowany filtr najmniejszej średniej kwadratów

%%  Wybranie wielkosci kroku
[mumaxlms,mumaxmselms]   = maxstep(hlms,Signal_2)  % Maksymalny krok dla filtrow adaptacyjnych
[mumaxnlms,mumaxmsenlms] = maxstep(hnlms);  

%%  Ustawienei wielkosci kroku
hlms.StepSize  = mumaxmselms/30;    
hnlms.StepSize = mumaxmsenlms/20;   
                                  
%%  Projektowanie filtrow adaptacyjnych
[ylms,elms] = filter(hlms,v2,Signal_2);    
[ynlms,enlms] = filter(hnlms,v2,Signal_2);

 

 

Filtr optymalny

Filtr optymalny Wienera bazuje na splocie funkcji korelacji. Oznacza to, że nie bierze on pod uwagę przebiegu sygnału tylko momenty statystyczne, przez co wycięcie szumu białego jest o wiele łatwiejsze.

%%  Filtr optymalny
bw = firwiener(L-1,v2,Signal_2);           % Optymalny filtr FIR Wienera
yw = filter(bw,1,v2);               % Estymata sygnalu zaszumionego
ew = Signal_2 - yw;                        % Estymata rzeczywistego sygnalu

 

Podsumowanie

Analizując wykresy FFT dla poszczególnych filtrów można zauważyć, że optymalny filtr Wienera doskonale poradził sobie z zadaniem usuwając większość szumów oraz pozostawiając składowe harmoniczne sygnału.