Como aplico um filtro Chebishev?

7

Eu li um artigo sobre uma interface cérebro-computador. Neste artigo, os autores relataram "cada sinal foi filtrado com um filtro Chebishev Tipo I de passagem de banda de 8 ordens, cujas frequências de corte são de 0,1 e 10 Hz e foi dizimado de acordo com a alta frequência de corte". Eu tentei projetar esse filtro com scipy:

import scipy.signal as signal
signal.cheby1(8,0.05,[0.1,10.0],btype='band',analog=0,output='ba')

O resultado foi:

Warning: invalid value encountered in sqrt
(array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
    nan,  nan,  nan,  nan,  nan,  nan]), array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
    nan,  nan,  nan,  nan,  nan,  nan]))

Não tenho experiência no processamento de sinais, então, na verdade, não sei o que estou fazendo. Não sei se eles usaram um filtro IIR ou FIR ou se tenho que escalar as frequências de corte ou se estou usando a ondulação errada. Espero que você possa me ajudar.

alfa
fonte

Respostas:

5

O principal problema com o exemplo que você deu é que a função de design do filtro cheby1está retornando todos os NaNs, o que não será um filtro muito bom. O problema é como você está especificando as frequências de borda da faixa de passagem / banda de parada. Esta função específica é destinada a emular a cheby1função do MATLAB ; as frequências fornecidas devem ser normalizadas, de modo que um valor 1igual a metade da taxa de amostragem.

import scipy.signal as signal
fs = whatever_the_sample_rate_of_the_filter_input_is_going_to_be
signal.cheby1(8,0.05,[0.1/(fs/2),10.0/(fs/2)],btype='band',analog=0,output='ba')

Eu não tenho o SciPy à mão, mas isso deve pelo menos projetar corretamente o filtro que você deseja.

Jason R
fonte
Obrigado, vou tentar isso à noite. Minha taxa de amostragem é de 240 Hz. Você sabe como aplicar o filtro nos dados? Eu encontrei uma fórmula na wikipedia ( en.wikipedia.org/wiki/… ), mas esta me dá resultados oscilantes estranhos (números muito grandes ...). Então eu acho que implementei algo errado. Eu implementei isso em C ++, então não posso chamar a função scipy.signal.lfilter. Btw. o que mudará quando eu alterar a ondulação?
Alfa #
2
Este é um filtro numericamente desafiador, pois você possui pólos muito próximos ao círculo unitário. Você precisa dividir o filtro em seções de segunda ordem e aplicá-las sequencialmente.
Hilmar
E como eu dividiria o filtro em seções de segunda ordem? Eu acho que o Matlab tem uma função para isso, mas não está disponível no scipy.
Alfa #
@alfa: não scipy não tem isso, mas eu comecei a traduzir o código oitava para python: gist.github.com/endolith/4525003
endolith
2

Duas frequências de corte normalmente significam que é um filtro passa-banda com o passa-alto a 0,1 Hz e o passa-baixo a 10 Hz. O ponto de corte passa-baixo (que é a mais alta das duas frequências) determina quanto você pode obter uma amostra abaixo. Se o seu filtro passa-baixo fosse infinitamente íngreme, você poderia obter uma nova taxa de amostragem de 20 Hz (o dobro do ponto de corte). Como a inclinação é limitada, é necessário deixar uma faixa de proteção entre a frequência de corte e a nova frequência de Nyquist. Quanto você precisa depende da ordem do filtro e da quantidade de ruído de alias que você pode tolerar.

Neste exemplo específico, parece que eles tiveram uma amostragem reduzida por um fator de 12 ou mais ou menos, o que parece muito agressivo para mim e provavelmente resultará em muito ruído de alias.

Hilmar
fonte
Obrigado pelas explicações. Eu acho que a redução da amostragem por um fator de 12 faz sentido, porque há muito ruído com alta amplitude e frequência e o sinal no qual eles estavam interessados ​​tinha uma frequência muito baixa. Na verdade, eu queria reproduzir os resultados deles. Existem tutoriais práticos (com código) sobre esse tipo de filtro?
alfa