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

28

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:

Données de l'accéléromètre

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:

Signal filtré

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

Michael M
la source
Bienvenue chez DSP! La courbe rouge de votre deuxième image est-elle une version "lissée" des données originales (vertes)?
Phonon
La courbe rouge est (espérons-le!) Une courbe de vitesse générée à partir de fft suivie d'un filtrage, suivie d'une arithmétique oméga (divisée par 2 * pi f j), suivie par inv. fft
Michael M
1
Peut-être que si vous incluez une expression mathématique plus précise ou un pseudocode pour ce que vous avez fait, cela clarifierait un peu les choses.
Phonon
Ajouté maintenant, c'est la sensation générale du code ..
Michael M
1
Ma question serait: qu'espérez-vous voir dans les données? Vous ne saurez pas si vous avez une bonne approche à moins que vous ne sachiez quelque chose sur le signal sous-jacent que vous attendez après filtrage. De plus, le code que vous avez montré est déroutant. Bien que vous ne montriez pas l'initialisation du fztableau, il semble que vous appliquiez un filtre passe-haut à la place.
Jason R

Réponses:

13

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.μρ

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

Résultats:

entrez la description de l'image ici entrez la description de l'image ici entrez la description de l'image ici

Daniel R. Pipa
la source
Aimez vraiment cette réponse, allez-y et essayez-la. Désolé, il m'a fallu si longtemps pour répondre!
Michael M
Excellente réponse. Merci pour les détails. Je recherche la version C de ce code. Quelqu'un ici a-t-il porté ce code matlab vers C qu'il aimerait partager? Merci.
René Limberger du
Que signifie constante par morceaux?
tilaprimera
6

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.

Daniel R. Pipa
la source
1
Je ne suis pas d'accord avec votre première déclaration. Vous supposez que le signal d'intérêt couvre toute la bande passante de la séquence d'entrée, ce qui est peu probable. Il est toujours possible d'obtenir un meilleur rapport signal / bruit en utilisant un filtrage linéaire dans ce cas, en éliminant le bruit hors bande. Si le signal est fortement suréchantillonné, vous pouvez obtenir une grande amélioration avec une approche aussi simple.
Jason R
C'est vrai, et cela est réalisé par le filtre de Wiener, lorsque vous connaissez les statistiques de votre signal et de votre bruit.
Daniel R. Pipa
Bien que la théorie derrière le débruitage en ondelettes soit compliquée, la mise en œuvre est aussi simple que l'approche que vous avez décrite. Il implique uniquement des banques de filtres et des seuils.
Daniel R. Pipa
1
Je fais des recherches à ce sujet maintenant, publierai mes progrès ci-dessus, merci à vous deux et à Phonon pour toute votre aide jusqu'à présent!
Michael M
@DanielPipa Je n'ai pas accès aux packages matlab en question. Pouvez-vous fournir un document ou une autre référence décrivant la méthode qui correspond à votre code matlab.
John Robertson
0

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)
ryanjdillon
la source