Auswahl des richtigen Filters für Beschleunigungsmesserdaten

Ich bin ziemlich neu in DSP und habe einige Nachforschungen über mögliche Filter zum Glätten von Beschleunigungsmesserdaten in Python angestellt. Ein Beispiel für den Datentyp, den ich erlebe, ist in der folgenden Abbildung zu sehen:

Im Wesentlichen suche ich einen Rat, um diese Daten zu glätten und sie schließlich in Geschwindigkeit und Verschiebung umzuwandeln. Ich verstehe, dass Beschleunigungsmesser von Mobiltelefonen extrem laut sind.

Ich glaube nicht, dass ich im Moment einen Kalman-Filter verwenden kann, weil ich das Gerät nicht in die Hände kriegen kann, um auf das von den Daten erzeugte Rauschen zu verweisen.

FFT hat einige interessante Ergebnisse hervorgebracht. Einer meiner Versuche war, das Beschleunigungssignal zu feilen und dann die tiefen Frequenzen auf einen absoluten FFT-Wert von 0 zu bringen. Dann verwendete ich Omega-Arithmetik und inverse FFT, um eine Darstellung der Geschwindigkeit zu erhalten. Die Ergebnisse waren wie folgt:

Ist dies ein guter Weg, um Dinge zu erledigen? Ich versuche, das allgemeine Rauschen des Signals zu beseitigen, aber es müssen offensichtliche Spitzen wie etwa bei 80 Sekunden identifiziert werden.

Ich habe es auch satt, einen Tiefpassfilter für die ursprünglichen Beschleunigungsmesserdaten zu verwenden, was eine hervorragende Glättung bewirkt hat, aber ich bin mir nicht sicher, wohin ich von hier aus gehen soll. Jede Anleitung, wohin Sie von hier aus gehen, wäre wirklich hilfreich!

EDIT: Ein bisschen Code:

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

Im Grunde genommen habe ich eine FFT mit meinen Beschleunigungsmesserdaten durchgeführt, die Sz ergab. Dabei wurden die hohen Frequenzen mit einem einfachen Ziegelwandfilter herausgefiltert (ich weiß, dass dies nicht ideal ist). Dann verwende ich Omega-Arithmetik für die FFT der Daten. Vielen Dank auch an datageist für das Hinzufügen meiner Bilder in meinen Post :)

Wie @JohnRobertson in Bag of Tricks zum Entrauschen von Signalen unter Beibehaltung scharfer Übergänge hervorhob, ist das Entrauschen mit Total Variaton (TV) eine weitere gute Alternative, wenn Ihr Signal stückweise konstant ist. Dies kann bei den Beschleunigungsmesserdaten der Fall sein, wenn Ihr Signal zwischen verschiedenen Plateaus variiert.

$\mu$$\rho$

$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)];

Ergebnisse:

Das Problem ist, dass Ihr Rauschen ein flaches Spektrum hat. Wenn Sie weißes Gaußsches Rauschen annehmen (was sich als gute Annahme herausstellt), ist dessen Leistungsspektrumsdichte konstant. Grob gesagt bedeutet dies, dass Ihr Rauschen alle Frequenzen enthält. Aus diesem Grund ist ein Frequenzansatz, z. B. DFT- oder Tiefpassfilter, nicht gut. Was wären Ihre Grenzfrequenzen, da Ihr Rauschen über das gesamte Spektrum verteilt ist?

Eine Antwort auf diese Frage ist der Wiener-Filter, der die Kenntnis der Statistik Ihres Rauschens und Ihres gewünschten Signals erfordert. Grundsätzlich wird das verrauschte Signal (Signal + Rauschen) über den Frequenzen gedämpft, bei denen das Rauschen höher als Ihr Signal ist, und es wird verstärkt, bei denen erwartet wird, dass Ihr Signal höher als Ihr Rauschen ist.

Ich würde jedoch modernere Ansätze vorschlagen, die nichtlineare Verarbeitung verwenden, beispielsweise Wavelet-Denoising. Diese Methoden liefern hervorragende Ergebnisse. Grundsätzlich wird das verrauschte Signal zuerst in Wavelets zerlegt und dann werden kleine Koeffizienten auf Null gesetzt. Dieser Ansatz funktioniert (und DFT nicht) aufgrund der Mehrfachauflösung von Wavelets. Das heißt, das Signal wird separat in Frequenzbändern verarbeitet, die durch die Wavelet-Transformation definiert sind.

Geben Sie in MATLAB 'wavemenu' und dann 'SWT denoising 1-D' ein. Dann 'File', 'Example Analysis', 'Noisy signals', 'mit Haar auf Level 5, Noisy blocks'. In diesem Beispiel wird Haar Wavelet verwendet, was für Ihr Problem gut funktionieren sollte.

Ich bin nicht gut in Python, aber ich glaube, Sie können einige NumPy-Pakete finden, die Haar-Wavelet-Denoising durchführen.

Auf Anregung von Daniel Pipa habe ich mir das Wavelet-Denoising angeschaut und diesen hervorragenden Artikel von Francisco Blanco-Silva gefunden.

Hier habe ich seinen Python-Code für die Bildverarbeitung so geändert, dass er mit 2D-Daten (Beschleunigungsmesser) und nicht mit 3D-Daten (Bilddaten) arbeitet.

Beachten Sie , dass der Schwellenwert für das Soft-Thresholding in Franciscos Beispiel "aufgeholt" wird. Berücksichtigen Sie dies und passen Sie es an Ihre Anwendung an.

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)

Woher:

• wavelet- Stringname der zu verwendenden Wavelet-Form (siehe pywt.wavelist()zB 'haar')
• noise_sigma - Standardabweichung des Rauschens von Daten
• data - Array mit zu filternden Werten (z. B. X-, Y- oder Z-Achsendaten)