Choisir le bon filtre pour les données de l'accéléromètre

Je suis assez nouveau sur DSP et j'ai fait des recherches sur les filtres possibles pour lisser les données de l'accéléromètre en python. Un exemple du type de données que je vais expérimenter peut être vu dans l'image suivante:

Essentiellement, je cherche des conseils pour lisser ces données pour éventuellement les convertir en vitesse et en déplacement. Je comprends que les accéléromètres des téléphones portables sont extrêmement bruyants.

Je ne pense pas pouvoir utiliser un filtre de Kalman pour le moment car je ne peux pas saisir l'appareil pour référencer le bruit produit par les données (j'ai lu qu'il est essentiel de placer l'appareil à plat et de trouver la quantité de bruit de ces lectures?)

La FFT a produit des résultats intéressants. L'une de mes tentatives a été de FFT le signal d'accélération, puis de rendre les basses fréquences pour avoir une valeur FFT absolue de 0. Ensuite, j'ai utilisé l'arithmétique oméga et la FFT inverse pour obtenir un tracé de la vitesse. Les résultats sont les suivants:

Est-ce une bonne façon de procéder? J'essaie de supprimer la nature bruyante globale du signal, mais des pics évidents, comme à environ 80 secondes, doivent être identifiés.

Je suis également fatigué d'utiliser un filtre passe-bas sur les données de l'accéléromètre d'origine, ce qui a fait un excellent travail de lissage, mais je ne sais pas vraiment où aller à partir d'ici. Tout conseil sur la route à suivre à partir d'ici serait vraiment utile!

EDIT: Un peu de 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

Donc, essentiellement, ive a effectué une FFT sur mes données d'accéléromètre, donnant Sz, filtré les hautes fréquences à l'aide d'un simple filtre de mur de briques (je sais que ce n'est pas idéal). Ensuite, j'utilise l'arithmétique oméga sur la FFT des données. Merci également à datageist d'avoir ajouté mes images dans mon message :)

Comme l'a souligné @JohnRobertson dans Bag of Tricks for Denoising Signals tout en maintenant des transitions nettes, Total Variaton (TV) débruitage est une autre bonne alternative si votre signal est constant par morceaux. Cela peut être le cas pour les données de l'accéléromètre, si votre signal continue de varier entre les différents plateaux.

Vous trouverez ci-dessous un code Matlab qui effectue le débruitage TV dans un tel signal. Le code est basé sur l'article An Augmented Lagrangian Method for Total Variation Video Restoration . Les paramètres et doivent être ajustés en fonction du niveau de bruit et des caractéristiques du signal.$\mu$$\rho$

Si est le signal bruyant et est le signal à estimer, la fonction à minimiser est , où est l'opérateur des différences finies.$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)];

Résultats:

Le problème est que votre bruit a un spectre plat. Si vous supposez un bruit blanc gaussien (qui s'avère être une bonne hypothèse), sa densité de spectre de puissance est constante. En gros, cela signifie que votre bruit contient toutes les fréquences. C'est pourquoi toute approche en fréquence, par exemple DFT ou filtres passe-bas, n'est pas bonne. Quelles seraient vos fréquences de coupure puisque votre bruit est sur tout le spectre?

Une réponse à cette question est le filtre de Wiener, qui nécessite la connaissance des statistiques de votre bruit et de votre signal souhaité. Fondamentalement, le signal bruyant (signal + bruit) est atténué sur les fréquences où le bruit devrait être plus important que votre signal, et il est amplifié là où votre signal devrait être plus grand que votre bruit.

Cependant, je suggérerais des approches plus modernes qui utilisent un traitement non linéaire, par exemple le débruitage en ondelettes. Ces méthodes donnent d'excellents résultats. Fondamentalement, le signal bruyant est d'abord décomposé en ondelettes, puis de petits coefficients sont mis à zéro. Cette approche fonctionne (et la DFT ne fonctionne pas) en raison de la nature multi-résolution des ondelettes. C'est-à-dire que le signal est traité séparément dans des bandes de fréquences définies par la transformée en ondelettes.

Dans MATLAB, tapez «menu d'onde» puis «SWT débruitage 1-D». Puis 'Fichier', 'Exemple d'analyse', 'Signaux bruyants', 'avec Haar au niveau 5, Blocs bruyants'. Cet exemple utilise l'ondelette Haar, qui devrait fonctionner correctement pour votre problème.

Je ne suis pas bon en Python, mais je pense que vous pouvez trouver des packages NumPy qui effectuent le débruitage des ondelettes Haar.

Selon la suggestion de Daniel Pipa, j'ai jeté un œil au débruitage en ondelettes et j'ai trouvé cet excellent article de Francisco Blanco-Silva.

Ici, j'ai modifié son code Python pour le traitement d'image pour travailler avec des données 2D (accéléromètre) plutôt que 3D (image).

Notez que le seuil est "compensé" pour le seuillage progressif dans l'exemple de Francisco. Considérez ceci et modifiez pour votre application.

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)

Où:

• wavelet- nom de chaîne de la forme en ondelettes à utiliser (voir pywt.wavelist(), par exemple 'haar')
• noise_sigma - écart type du bruit par rapport aux données
• data - tableau de valeurs à filtrer (par exemple, les données des axes x, y ou z)