Escolhendo o filtro correto para dados do acelerômetro

28

Eu sou bastante novo no DSP e fiz algumas pesquisas sobre possíveis filtros para suavizar dados do acelerômetro em python. Um exemplo do tipo de dados que estou enfrentando pode ser visto na imagem a seguir:

Dados do acelerômetro

Essencialmente, estou procurando conselhos para suavizar esses dados e eventualmente convertê-los em velocidade e deslocamento. Entendo que os acelerômetros dos telefones celulares são extremamente barulhentos.

Eu não acho que posso usar um filtro Kalman no momento porque não consigo me apossar do dispositivo para fazer referência ao ruído produzido pelos dados (eu li que é essencial colocar o dispositivo na horizontal e encontrar a quantidade de ruído dessas leituras?)

A FFT produziu alguns resultados interessantes. Uma das minhas tentativas foi realizar a FFT do sinal de aceleração e, em seguida, renderizar as frequências baixas para ter um valor absoluto da FFT igual a 0. Então usei a FFT aritmética e inversa ômega para obter um gráfico de velocidade. Os resultados foram os seguintes:

Sinal filtrado

Essa é uma boa maneira de fazer as coisas? Estou tentando remover a natureza geral barulhenta do sinal, mas picos óbvios, como em cerca de 80 segundos, precisam ser identificados.

Também me cansei de usar um filtro passa-baixo nos dados originais do acelerômetro, o que fez um ótimo trabalho em suavizá-lo, mas não tenho muita certeza de onde ir a partir daqui. Qualquer orientação sobre onde ir a partir daqui seria realmente útil!

EDIT: Um pouco de código:

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

Então, essencialmente, eu fiz uma FFT nos dados do meu acelerômetro, fornecendo Sz, filtrando altas frequências usando um simples filtro de parede de tijolos (eu sei que não é o ideal). Em seguida, use aritmética ômega na FFT dos dados. Também muito obrigado ao datageist por adicionar minhas imagens no meu post :)

Michael M
fonte
Bem-vindo ao DSP! A curva vermelha na sua segunda foto é uma versão "suavizada" dos dados originais (verdes)?
Phonon
A curva vermelha é (espero!) Uma curva de velocidade gerada a partir de fft seguida de filtragem, seguida por aritmética ômega (dividida por 2 * pi f j), seguida por inv. fft
Michael M
1
Talvez se você incluir uma expressão matemática mais precisa ou pseudocódigo para o que você fez, isso esclarecerá um pouco as coisas.
Phonon
Adicionado alguns agora, isso é a sensação geral do código ..
Michael M
1
Minha pergunta seria: o que você espera ver nos dados? Você não saberá se tem uma boa abordagem, a menos que saiba algo sobre o sinal subjacente que espera ver após a filtragem. Além disso, o código que você mostrou é confuso. Embora você não mostre a inicialização da fzmatriz, parece que você está aplicando um filtro passa-alto.
Jason R

Respostas:

13

Conforme apontado por @JohnRobertson no Pacote de truques para diminuir o ruído de sinais durante a manutenção de transições nítidas, o denoising de Variaton total (TV) é outra boa alternativa se o sinal for constante em partes. Pode ser o caso dos dados do acelerômetro, se o seu sinal continuar variando entre diferentes platôs.

Abaixo está um código do Matlab que executa o denoising da TV nesse sinal. O código é baseado no artigo Um método lagrangiano aumentado para restauração de vídeo com variação total . Os parâmetros e devem ser ajustados de acordo com o nível de ruído e as características do sinal.μρ

Se é o sinal ruidoso e é o sinal a ser avaliado, a função a ser minimizada é , onde é o operador das diferenças finitas.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)];

Resultados:

insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui

Daniel R. Pipa
fonte
Realmente gosto desta resposta, vou em frente e experimente. Desculpe por demorar tanto em responsder!
28712 Michael M
Excelente resposta. Obrigado pelos detalhes. Estou procurando a versão C deste código. Alguém aqui portou esse código do matlab para C que gostaria de compartilhar? Obrigado.
Rene Limberger,
O que significa constante peça a peça?
Tilaprimera
6

O problema é que seu ruído tem um espectro plano. Se você assume um ruído gaussiano branco (o que acaba sendo uma boa suposição), sua densidade de espectro de potência é constante. Grosso modo, isso significa que seu ruído contém todas as frequências. É por isso que qualquer abordagem de frequência, por exemplo, DFT ou filtros passa-baixo, não é boa. Quais seriam suas frequências de corte, já que seu ruído está em todo o espectro?

Uma resposta para essa pergunta é o filtro Wiener, que requer conhecimento das estatísticas do seu ruído e do sinal desejado. Basicamente, o sinal barulhento (sinal + ruído) é atenuado nas frequências em que se espera que o ruído seja maior do que o seu sinal e é amplificado onde seu sinal é maior do que o seu ruído.

No entanto, eu sugeriria abordagens mais modernas que usam processamento não linear, por exemplo, waenoet denoising. Esses métodos fornecem excelentes resultados. Basicamente, o sinal ruidoso é primeiro decomposto em wavelets e depois pequenos coeficientes são zerados. Essa abordagem funciona (e o DFT não) devido à natureza de multi-resolução das wavelets. Ou seja, o sinal é processado separadamente nas bandas de frequência definidas pela transformada wavelet.

No MATLAB, digite 'wavemenu' e depois 'SWT denoising 1-D'. Em seguida, 'Arquivo', 'Análise de exemplo', 'Sinais ruidosos', 'com Haar no nível 5, Blocos ruidosos'. Este exemplo usa a wavelet Haar, que deve funcionar bem para o seu problema.

Eu não sou bom em Python, mas acredito que você pode encontrar alguns pacotes NumPy que executam denoising Haar wavelet.

Daniel R. Pipa
fonte
1
Eu discordo de sua primeira declaração. Você está assumindo que o sinal de interesse cobre toda a largura de banda da sequência de entrada, o que é improvável. Ainda é possível obter uma relação sinal / ruído aprimorada usando filtragem linear neste caso, eliminando o ruído fora da banda. Se o sinal for altamente amostrado, você poderá obter uma grande melhoria com uma abordagem tão simples.
Jason R
É verdade, e isso é alcançado pelo filtro Wiener, quando você conhece as estatísticas do seu sinal e do seu ruído.
Daniel Daniel Pipa
Embora a teoria por trás do denoising wavelet seja complicada, a implementação é tão simples quanto a abordagem que você descreveu. Envolve apenas bancos de filtros e limiares.
Daniel Daniel Pipa
1
Estou pesquisando sobre isso agora, postarei meu progresso acima, graças a vocês e a Phonon por toda a sua ajuda até agora!
Michael M
@DanielPipa Não tenho acesso aos pacotes matlab em questão. Você pode fornecer um documento ou outra referência que descreva o método que corresponde ao seu código matlab.
9788 John Robertson
0

Conforme sugestão de Daniel Pipa, dei uma olhada no wavelet denoising e encontrei este excelente artigo de Francisco Blanco-Silva.

Aqui, modifiquei seu código Python para que o processamento de imagens funcione com dados 2D (acelerômetro) em vez de dados 3D (imagem).

Observe que o limite é "compensado" pelo limiar suave no exemplo de Francisco. Considere isso e modifique para seu aplicativo.

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)

Onde:

  • wavelet- nome da string da forma de wavelet a ser usada (veja pywt.wavelist(), por exemplo 'haar')
  • noise_sigma - desvio padrão de ruído dos dados
  • data - matriz de valores para filtrar (por exemplo, dados do eixo x, y ou z)
ryanjdillon
fonte