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

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:

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:

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 :)

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.$\mu$$\rho$

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.$y$$x$$\mu\|{x-y}\|^2+\|{Dx}\|_1$$D$

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:

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.

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)