Como escrever filtro lowpass para sinal amostrado em Python?

16

Eu tenho algum sinal que amostrou cada 1 ns (1e-9 segundos) e tem, digamos, 1e4 pontos. Preciso filtrar altas frequências desse sinal. Digamos que eu precise filtrar frequências maiores que 10 MHz. Eu quero que para frequências menores que o sinal de frequência de corte seja passado inalterado. Isso significa que o ganho do filtro será 1 para frequências inferiores à frequência de corte. Gostaria de poder especificar a ordem dos filtros. Quero dizer, o filtro de primeira ordem tem uma inclinação de 20 db / década (queda de energia) após a frequência de corte, o filtro de segunda ordem tem uma inclinação de 40 db / dec após a frequência de corte e assim por diante. Alto desempenho do código é importante.

Alex
fonte

Respostas:

19

A resposta de frequência para o filtro projetado usando a função de manteiga é:

Resposta do filtro Butterworth

Mas não há razão para limitar o filtro a um design de filtro monotônico constante. Se você deseja uma atenuação mais alta na faixa de parada e na faixa de transição mais íngreme, existem outras opções. Para obter mais informações sobre a especificação de um filtro usando iirdesing, consulte isso . Como mostra os gráficos de resposta em frequência para o design da manteiga , a frequência de corte (ponto de -3dB) está longe do objetivo. Isso pode ser atenuado pela amostragem reduzida antes da filtragem (as funções de design terão um tempo difícil com um filtro tão estreito, 2% da largura de banda). Vamos analisar a taxa de amostragem original com o limite especificado.

import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

from scipy.signal import fir_filter_design as ffd
from scipy.signal import filter_design as ifd

# setup some of the required parameters
Fs = 1e9           # sample-rate defined in the question, down-sampled

# remez (fir) design arguements
Fpass = 10e6       # passband edge
Fstop = 11.1e6     # stopband edge, transition band 100kHz
Wp = Fpass/(Fs)    # pass normalized frequency
Ws = Fstop/(Fs)    # stop normalized frequency

# iirdesign agruements
Wip = (Fpass)/(Fs/2)
Wis = (Fstop+1e6)/(Fs/2)
Rp = 1             # passband ripple
As = 42            # stopband attenuation

# Create a FIR filter, the remez function takes a list of 
# "bands" and the amplitude for each band.
taps = 4096
br = ffd.remez(taps, [0, Wp, Ws, .5], [1,0], maxiter=10000) 

# The iirdesign takes passband, stopband, passband ripple, 
# and stop attenuation.
bc, ac = ifd.iirdesign(Wip, Wis, Rp, As, ftype='ellip')  
bb, ab = ifd.iirdesign(Wip, Wis, Rp, As, ftype='cheby2') 

Filtros de taxa de amostra original

Como mencionado, como estamos tentando filtrar uma porcentagem tão pequena da largura de banda, o filtro não terá um corte agudo. Nesse caso, filtro lowpass, podemos reduzir a largura de banda para obter um filtro com melhor aparência. A função de reamostragem python / scipy.signal pode ser usada para reduzir a largura de banda.

Observe que a função de reamostragem realizará a filtragem para evitar aliases. A pré-filtragem também pode ser realizada (para reduzir o aliasing) e, nesse caso, poderíamos simplesmente re-amostrar em 100 e pronto , mas a pergunta foi feita sobre a criação de filtros. Neste exemplo, reduziremos a amostra em 25 e criaremos um novo filtro

R = 25;            # how much to down sample by
Fsr = Fs/25.       # down-sampled sample rate
xs = signal.resample(x, len(x)/25.)

Se atualizarmos os parâmetros de design para o filtro FIR, a nova resposta será.

# Down sampled version, create new filter and plot spectrum
R = 25.             # how much to down sample by
Fsr = Fs/R          # down-sampled sample rate
Fstop = 11.1e6      # modified stopband
Wp = Fpass/(Fsr)    # pass normalized frequency
Ws = Fstop/(Fsr)    # stop normalized frequency
taps = 256
br = ffd.remez(taps, [0, Wp, Ws, .5], [1,0], maxiter=10000) 

Resposta de filtro com redução de amostra

O filtro que opera nos dados de amostragem reduzida tem uma resposta melhor. Outro benefício do uso de um filtro FIR é que você terá resposta de fase linear.

Christopher Felton
fonte
1
Obrigado. Como você cria um gráfico do espectro de sinal?
1212 Alex
Muito obrigado por uma excelente resposta! Gostaria de saber se você poderia explicar como aplicar um filtro FIR com base nos coeficientes calculados usando o Remez? Estou tendo problemas para entender o que filtfiltdeseja para o aparâmetro.
23413 ali13m
Depois de ter os coeficientes de um projeto do filtro, ( b para FIR b e um para IIR), você pode utilizar um par de funções diferentes para realizar a filtragem: lfilter , convolve , filtfilt . Normalmente, todas essas funções operam de maneira semelhante: y = filtfilt (b, a, x) Se você tiver um filtro FIR, basta definir a = 1 , x é o sinal de entrada, b é o coeficiente FIR. Esta postagem também pode ajudar.
Christopher Felton
5

Isto funciona?

from __future__ import division
from scipy.signal import butter, lfilter

fs = 1E9 # 1 ns -> 1 GHz
cutoff = 10E6 # 10 MHz
B, A = butter(1, cutoff / (fs / 2), btype='low') # 1st order Butterworth low-pass
filtered_signal = lfilter(B, A, signal, axis=0)

Você está certo, porém, a documentação não é muito completa. Parece que butteré um invólucro para iirfilter, que está melhor documentado :

N: int A ordem do filtro. Wn: array_like Uma sequência escalar ou comprimento-2 que fornece as frequências críticas.

Porém, a maioria dessas coisas é clonada do matlab, para que você possa ver também a documentação deles :

a frequência de corte normalizada Wn deve ser um número entre 0 e 1, em que 1 corresponde à frequência de Nyquist, π radianos por amostra.

Atualizar:

Eu adicionei documentação para essas funções. :) O Github facilita isso.

endólito
fonte
1

Não sabe ao certo qual é o seu aplicativo, mas você pode conferir o Gnuradio: http://gnuradio.org/doc/doxygen/classgr__firdes.html

Os blocos de processamento de sinal são escritos em C ++ (embora os gráficos de fluxo do Gnuradio estejam em Python), mas você disse que o alto desempenho é importante.

wrapperapps
fonte
1

Estou tendo bons resultados com este filtro FIR. Observa que aplica o filtro duas vezes, avançando e retrocedendo, de modo a compensar o deslocamento do sinal (a filtfiltfunção não funcionou, não sei por quê):

def firfilt(interval, freq, sampling_rate):
    nfreq = freq/(0.5*sampling_rate)
    taps =  sampling_rate + 1
    a = 1
    b = scipy.signal.firwin(taps, cutoff=nfreq)
    firstpass = scipy.signal.lfilter(b, a, interval)
    secondpass = scipy.signal.lfilter(b, a, firstpass[::-1])[::-1]
    return secondpass

Um ótimo recurso para filtrar o design e o uso, de onde peguei esse código e de onde exemplos de filtros passa-banda e passe-hi podem ser obtidos, é ISTO .

heltonbiker
fonte
Não acredito que exista muito benefício em filtrar para frente e para trás um filtro FIR. Um filtro IIR pode se beneficiar de avançar / reverter (filtfilt) porque você pode obter fase linear de um filtro de fase não linear por meio de filtragem reversa.
Christopher Felton
2
@ChristopherFelton Acabo de reverter para sincronizar um sinal eletromiográfico RAW com a versão suavizada de si mesmo. Eu sei que poderia mudar o sinal, mas filtrar duas vezes acaba sendo menos problemático. Vale notar que a segunda passagem quase não altera a primeira passagem já filtrada ... Obrigado por observar!
heltonbiker
Ah sim. Para remover o atraso (atraso do grupo), bom ponto.
Christopher Felton
1

Não tenho direitos de comentário ...

@ endolith: Eu uso o mesmo que você, exceto usando o scipy.signal.filtfilt (B, A, x) onde x é o vetor de entrada a ser filtrado - por exemplo, numpy.random.normal (size = (N)) . filtfilt faz uma passagem para frente e para trás do sinal. Por uma questão de integridade (a maioria é igual a @endolith):

import numpy as np
import scipy.signal as sps

input = np.random.normal(size=(N)) # Random signal as example
bz, az = sps.butter(FiltOrder, Bandwidth/(SamplingFreq/2)) # Gives you lowpass Butterworth as default
output = sps.filtfilt(bz, az, input) # Makes forward/reverse filtering (linear phase filter)

filtfilt, como também sugerido por @heltonbiker, exige matrizes de coeficientes, acredito. Caso você precise executar a filtragem de banda em banda básica complexa, é necessária uma configuração mais envolvida, mas isso não parece ser um problema aqui.

Lars1
fonte