Wybór właściwego filtra dla danych akcelerometru

28

Jestem dość nowy w DSP i przeprowadziłem badania dotyczące możliwych filtrów wygładzających dane akcelerometru w pythonie. Przykład rodzaju danych, których doświadczam, można zobaczyć na poniższym obrazie:

Dane akcelerometru

Zasadniczo szukam porady, jak wygładzić te dane, aby ostatecznie przekształcić je w prędkość i przemieszczenie. Rozumiem, że akcelerometry z telefonów komórkowych są bardzo głośne.

Nie sądzę, żebym mógł w tej chwili użyć filtra Kalmana, ponieważ nie mogę złapać urządzenia, aby odwoływać się do hałasu wytwarzanego przez dane (czytam, że konieczne jest umieszczenie urządzenia płaskiego i znalezienie poziomu hałasu z tych odczytów?)

FFT przyniosło kilka interesujących wyników. Jedną z moich prób było FFT sygnału przyspieszenia, a następnie renderowanie niskich częstotliwości w celu uzyskania bezwzględnej wartości FFT równej 0. Następnie zastosowałem arytmetykę omega i odwrotną FFT, aby uzyskać wykres prędkości. Wyniki były następujące:

Filtrowany sygnał

Czy to dobry sposób na załatwienie różnych spraw? Próbuję usunąć ogólną hałaśliwą naturę sygnału, ale należy zidentyfikować oczywiste piki, takie jak około 80 sekund.

Zmęczyłem się również używaniem filtra dolnoprzepustowego na oryginalnych danych akcelerometru, który wykonał świetną robotę, aby go wygładzić, ale tak naprawdę nie jestem pewien, dokąd się udać. Wszelkie wskazówki na temat tego, gdzie się udać, byłyby naprawdę pomocne!

EDYCJA: Trochę kodu:

for i in range(len(fz)): 
    testing = (abs(Sz[i]))/Nz

    if fz[i] < 0.05:
        Sz[i]=0

Velfreq = []
Velfreqa = array(Velfreq)
Velfreqa = Sz/(2*pi*fz*1j)
Veltimed = ifft(Velfreqa)
real = Veltimed.real

Zasadniczo więc przeprowadziłem FFT na danych mojego akcelerometru, dając Sz, odfiltrowałem wysokie częstotliwości za pomocą prostego filtru z cegły (wiem, że nie jest idealny). Następnie używam arytmetyki omega na FFT danych. Również bardzo dziękuję datageist za dodanie moich zdjęć do mojego posta :)

Michael M.
źródło
Witamy w DSP! Czy czerwona krzywa na drugim zdjęciu jest „wygładzoną” wersją oryginalnych (zielonych) danych?
Phonon
Czerwona krzywa jest (miejmy nadzieję!) Krzywą prędkości wygenerowaną z fft, po której następuje filtrowanie, a następnie arytmetyka omega (dzielona przez 2 * pi f j), po której następuje inv. fft
Michael M.
1
Być może, jeśli podasz bardziej precyzyjne wyrażenie matematyczne lub pseudokod tego, co zrobiłeś, nieco to wyjaśni.
Phonon
Dodano trochę teraz, to jest ogólny sposób działania kodu.
Michael M
1
Moje pytanie brzmiałoby: czego oczekujesz w danych? Nie będziesz wiedział, czy masz dobre podejście, chyba że wiesz coś o podstawowym sygnale, którego oczekujesz po filtrowaniu. Ponadto kod, który pokazałeś, jest mylący. Chociaż nie wyświetlasz inicjalizacji fztablicy, wydaje się, że zamiast tego stosujesz filtr górnoprzepustowy.
Jason R

Odpowiedzi:

13

Jak zauważył @JohnRobertson w Bag of Tricks for Denoising Signals Podczas utrzymywania ostrych przejść, denoising Total Variaton (TV) jest kolejną dobrą alternatywą, jeśli twój sygnał jest stały w kawałkach. Może tak być w przypadku danych akcelerometru, jeśli twój sygnał zmienia się w zależności od różnych płaskości.

Poniżej znajduje się kod Matlab, który wykonuje denoising TV w takim sygnale. Kod oparty jest na pracy An Augmented Lagrangian Method for Total Variation Video Restoration . Parametry i należy dostosować do poziomu hałasu i charakterystyki sygnału.μρ

Jeśli jest sygnałem zaszumionym, a jest sygnałem, który ma być oszacowany, funkcją do zminimalizowania jest , gdzie jest operatorem różnic skończonych.yxμxy2+Dx1D

function denoise()

f = [-1*ones(1000,1);3*ones(100,1);1*ones(500,1);-2*ones(800,1);0*ones(900,1)];
plot(f);
axis([1 length(f) -4 4]);
title('Original');
g = f + .25*randn(length(f),1);
figure;
plot(g,'r');
title('Noisy');
axis([1 length(f) -4 4]);
fc = denoisetv(g,.5);
figure;
plot(fc,'g');
title('De-noised');
axis([1 length(f) -4 4]);

function f = denoisetv(g,mu)
I = length(g);
u = zeros(I,1);
y = zeros(I,1);
rho = 10;

eigD = abs(fftn([-1;1],[I 1])).^2;
for k=1:100
    f = real(ifft(fft(mu*g+rho*Dt(u)-Dt(y))./(mu+rho*eigD)));
    v = D(f)+(1/rho)*y;
    u = max(abs(v)-1/rho,0).*sign(v);
    y = y - rho*(u-D(f));
end

function y = D(x)
y = [diff(x);x(1)-x(end)];

function y = Dt(x)
y = [x(end)-x(1);-diff(x)];

Wyniki:

wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj

Daniel R. Pipa
źródło
Naprawdę podoba mi się ta odpowiedź, idź i spróbuj. Przepraszam, że tyle czasu zajęło mi udzielenie odpowiedzi!
Michael M
Doskonała odpowiedź. Dziękuję za szczegóły. Szukam wersji C tego kodu. Czy ktoś tu przetransportował ten kod Matlab do C, którym chciałby się podzielić? Dzięki.
Rene Limberger,
Co oznacza stała kawałek?
tilaprimera
6

Problem polega na tym, że twój hałas ma płaskie spektrum. Jeśli przyjmiesz biały szum Gaussa (który okazuje się być dobrym założeniem), jego gęstość widma mocy jest stała. Z grubsza mówiąc, oznacza to, że twój hałas zawiera wszystkie częstotliwości. Dlatego każde podejście częstotliwościowe, np. DFT lub filtry dolnoprzepustowe, nie jest dobre. Jakie byłyby twoje częstotliwości odcięcia, ponieważ twój hałas jest w całym spektrum?

Jedną odpowiedzią na to pytanie jest filtr Wienera, który wymaga znajomości statystyki twojego hałasu i pożądanego sygnału. Zasadniczo zaszumiony sygnał (sygnał + szum) jest tłumiony w stosunku do częstotliwości, na których oczekuje się, że hałas będzie większy niż twój sygnał, i jest wzmacniany tam, gdzie oczekuje się, że twój sygnał będzie większy niż twój szum.

Sugerowałbym jednak bardziej nowoczesne podejścia wykorzystujące przetwarzanie nieliniowe, na przykład denoising falkowy. Te metody zapewniają doskonałe wyniki. Zasadniczo zaszumiony sygnał jest najpierw rozkładany na falki, a następnie małe współczynniki są zerowane. To podejście działa (a DFT nie działa) ze względu na charakter fal o wielu rozdzielczościach. Oznacza to, że sygnał jest przetwarzany osobno w pasmach częstotliwości określonych przez transformację falkową.

W MATLAB wpisz „wavemenu”, a następnie „SWT denoising 1-D”. Następnie „Plik”, „Przykładowa analiza”, „Zaszumione sygnały”, „z Haar na poziomie 5, zaszumione bloki”. W tym przykładzie użyto falki Haar, która powinna dobrze działać w przypadku Twojego problemu.

Nie jestem dobry w Pythonie, ale wierzę, że można znaleźć niektóre pakiety NumPy, które wykonują odszumianie falkowe Haar.

Daniel R. Pipa
źródło
1
Nie zgodziłbym się z twoim pierwszym oświadczeniem. Zakładasz, że sygnał zainteresowania obejmuje pełne pasmo sekwencji wejściowej, co jest mało prawdopodobne. W tym przypadku nadal można uzyskać lepszy stosunek sygnału do szumu za pomocą filtrowania liniowego, eliminując szum pozapasmowy. Jeśli sygnał jest bardzo nadpróbkowany, możesz uzyskać dużą poprawę dzięki tak prostemu podejściu.
Jason R
To prawda, a osiąga się to dzięki filtrowi Wienera, gdy znasz statystyki swojego sygnału i szumu.
Daniel R. Pipa,
Chociaż teoria denoisingu falkowego jest skomplikowana, wdrożenie jest tak proste, jak opisane podejście. Dotyczy tylko banków filtrów i progów.
Daniel R. Pipa,
1
Badam to teraz, opublikuję moje postępy powyżej, dzięki wam i Phononowi za całą waszą dotychczasową pomoc!
Michael M
@DanielPipa Nie mam dostępu do omawianych pakietów Matlab. Czy możesz podać dokument lub inną dokumentację opisującą metodę, która odpowiada Twojemu kodowi Matlab.
John Robertson,
0

Zgodnie z sugestią Daniela Pipy rzuciłem okiem na denerwowanie falek i znalazłem ten znakomity artykuł Francisco Blanco-Silvy.

Tutaj zmodyfikowałem jego kod Pythona do przetwarzania obrazu, aby działał z danymi 2D (akcelerometr), a nie danymi 3D (obraz).

Zauważ , że próg jest „nadrobiony” dla miękkiego progowania w przykładzie Francisco. Rozważ to i zmodyfikuj w swojej aplikacji.

def wavelet_denoise(data, wavelet, noise_sigma):
    '''Filter accelerometer data using wavelet denoising

    Modification of F. Blanco-Silva's code at: https://goo.gl/gOQwy5
    '''
    import numpy
    import scipy
    import pywt

    wavelet = pywt.Wavelet(wavelet)
    levels  = min(15, (numpy.floor(numpy.log2(data.shape[0]))).astype(int))

    # Francisco's code used wavedec2 for image data
    wavelet_coeffs = pywt.wavedec(data, wavelet, level=levels)
    threshold = noise_sigma*numpy.sqrt(2*numpy.log2(data.size))

    new_wavelet_coeffs = map(lambda x: pywt.threshold(x, threshold, mode='soft'),
                             wavelet_coeffs)

    return pywt.waverec(list(new_wavelet_coeffs), wavelet)

Gdzie:

  • wavelet- nazwa ciągu formy falkowej do użycia (patrz pywt.wavelist()np. 'haar')
  • noise_sigma - standardowe odchylenie hałasu od danych
  • data - tablica wartości do filtrowania (np. dane osi x, y lub z)
ryanjdillon
źródło