filtro passa-baixo e FFT para iniciantes com Python

23

Eu sou novo no processamento de sinais e principalmente na FFT; portanto, não tenho certeza se estou fazendo a coisa correta aqui e estou um pouco confuso com o resultado.

Eu tenho uma função real discreta (dados de medição) e quero configurar um filtro passa-baixo nisso. A ferramenta escolhida é o Python com o pacote numpy. Eu sigo este procedimento:

  • calcular o fft da minha função
  • cortar altas frequências
  • executar o fft inverso

Aqui está o código que estou usando:

import numpy as np
sampling_length = 15.0*60.0 # measured every 15 minutes
Fs = 1.0/sampling_length
ls = range(len(data)) # data contains the function
freq = np.fft.fftfreq(len(data), d = sampling_length)
fft = np.fft.fft(data)
x = freq[:len(data)/2] 
for i in range(len(x)):
if x[i] > 0.005: # cut off all frequencies higher than 0.005
    fft[i] = 0.0
    fft[len(data)/2 + i] = 0.0
inverse = np.fft.ifft(fft)

Esse é o procedimento correto? O resultado inversecontém valores complexos, que me confundem.

Até B
fonte
1
Quando eu estava aprendendo FFT, achei este post muito útil. glowingpython.blogspot.com/2011/08/...
David Poole

Respostas:

23

O fato de o resultado ser complexo é de se esperar. Quero destacar algumas coisas:

Você está aplicando um filtro no domínio da frequência da parede de tijolos aos dados, tentando zerar todas as saídas FFT que correspondem a uma frequência maior que 0,005 Hz e depois transformando inversamente para obter um sinal no domínio do tempo novamente. Para que o resultado seja real, a entrada na FFT inversa deve ser conjugada simétrica . Isso significa que, por um comprimento - FFT,N

X[k]=X[Nk],k=1,2,,N21(Neven)

X[k]=X[Nk],k=1,2,,N2(Nodd)
  • Observe que, para pares, e não são iguais em geral, mas são reais. Para ímpar , deve ser real.NX[0]NX[0]X[N2]NX[0]

Vejo que você tentou fazer algo assim no seu código acima, mas não está totalmente correto. Se você aplicar a condição acima no sinal transmitido para a FFT inversa, deverá receber um sinal real.

Meu segundo ponto é mais filosófico: o que você está fazendo funcionará, pois suprimirá o conteúdo do domínio da frequência que você não deseja. No entanto, isso não é tipicamente o modo como um filtro passa-baixo seria implementado na prática. Como mencionei antes, o que você está fazendo é essencialmente aplicar um filtro que tenha uma resposta de magnitude na parede de tijolos (ou seja, perfeitamente retangular). A resposta ao impulso de um filtro desse tipo tem uma forma . Como a multiplicação no domínio da frequência é equivalente a (no caso de usar a DFT, circular) convolução no domínio do tempo, essa operação é equivalente a convolver o sinal do domínio do tempo com uma função .s i n csinc(x)sinc

Por que isso é um problema? Lembre-se de como é a função no domínio do tempo (abaixo da imagem, vergonhosamente emprestada da Wikipedia):sinc

trama da função sinc

A função possui amplo suporte no domínio do tempo; decai muito lentamente à medida que você se afasta no tempo do lóbulo principal. Para muitas aplicações, essa não é uma propriedade desejável; quando você envolve um sinal com um , os efeitos dos lóbulos laterais decadentes lentamente serão aparentes na forma no domínio do tempo do sinal de saída filtrado. Esse tipo de efeito costuma ser chamado de toque . Se você sabe o que está fazendo, há alguns casos em que esse tipo de filtragem pode ser apropriado, mas no caso geral, não é o que você deseja.s i n csincsinc

Existem meios mais práticos de aplicar filtros passa-baixo, nos domínios de tempo e frequência. A resposta de impulso finito e os filtros de resposta de impulso infinito podem ser aplicados diretamente usando sua representação na equação de diferença . Ou, se o seu filtro tiver uma resposta de impulso suficientemente longa, você poderá obter benefícios de desempenho usando técnicas de convolução rápidas baseadas na FFT (aplicando o filtro multiplicando no domínio da frequência em vez de convolução no domínio do tempo), como a sobreposição salvar e sobrepor-adicionar métodos.

Jason R
fonte
A função sinc é a filtragem ideal, não? É para isso que todos os outros filtros estão buscando, mas não alcançando. É ruim para o processamento de imagens, porque as imagens não são suavizadas primeiro, então produz um toque que parece horrível, mas para áudio ou outros sinais que foram filtrados antes da amostragem, não é o melhor filtro que você pode obter?
endolith 29/11
1
Sim, meu resultado não foi conjugado simétrico. Corrigi o código, agora tudo funciona bem. Obrigado!
2
3
@endolith - um Sinc é um interpolador ideal para certains tipos de interpolação, mas está longe de ser ideal como um filtro para a maioria dos tipos de requisitos de filtro comuns, tais como nivelamento de resposta banda de passagem, a rejeição banda stop, e etc.
hotpaw2
+1 para a boa explicação sobre "por que as pessoas não implementar o filtro como o PO faz"
Sibbs Gambling
Você tem que usar um sinc em janelas. Se você não tem restrições de tempo, esse é o filtro ideal, muito melhor que o Chebichev.