Filtro de áudio DFT
Filtro de áudio DFT
Ok, fiquei cansado de usar os filtros do Audacity, então é hora de construir nosso próprio filtro passa-baixas. Até onde sei, há três métodos básicos de desenhar um filtro digital, como o filtro passa-baixas que precisamos para a modulação AM. O mais comum é o filtro FIR. Outro método é o filtro IIR.
E o terceiro método é a Transformada Discreta de Fourier (DFT). Neste artigo, vou utilizar o método mais fácil, que é justamente DFT. Outros artigos abordam FIR e IIR.
DFT é um conceito matemático; é a versão "digital" da transformada de Fourier. A implementação mais comumente utilizada chama-se FFT (transformada rápida de Fourier). A a quase totalidade das bibliotecas matemáticas oferece FFT. A biblioteca Python-Numeric não é exceção.
DFT converte um sinal no domínio do tempo, por exemplo um arquivo WAV, para uma representação no domínio da freqüência. O analisador de espectro (figura mais acima) é quase a mesma coisa, embora o DFT costume ser calculado com muito mais resolução que as meras 18 barras (por canal) da figura.
Ou seja, o DFT permite dissecar um áudio em graves, médios, agudos, até por nota musical — bem como inúmeras outras aplicações.
DFT não é um processamento unidirecional. Ele pode ser invertido. Dada uma representação no domínio da freqüência, podemos executar a transformada inversa (IDFT) e obter a onda no domínio do tempo correspondente. Então, nós podemos pegar um áudio qualquer, transformar com DFT, manipular o espectro (aumentar graves, diminuir agudos, etc.) e converter de volta para áudio WAV.
Em nosso caso, queremos construir filtros passa-baixas e passa-altas, então só precisamos eliminar do espectro as freqüências que não desejamos. A rigor, um corte brusco não é a melhor forma de configurar um filtro, pois causa fortes distorções de fase nas freqüências próximas aos cortes, mas é bom o suficiente para demonstração.
Converter de/para o domínio das freqüências parece trabalho demais. De fato, os métodos FIR e IIR permitem fazer essa filtragem diretamente no domínio do tempo, e normalmente seriam as escolhas certas em aplicações práticas. Porém eles não são tão fáceis de entender, nem tão versáteis.
Neste post, utilizarei a voz de uma cantora de que gosto muito: Núbia Lafayette. Outra razão de ter escolhido este áudio é porque eu possuo o CD, então ninguém pode reclamar que estou fazendo pirataria. Vamos ouvir o original:
Agora, a versão filtrada (filtro passa-banda 1000 a 4000Hz), que soa como aqueles radinhos de pilha:
E o código que fez a mágica:
#!/usr/bin/env python
LOWPASS = 4000 # Hz
HIGHPASS = 1000 # Hz
SAMPLE_RATE = 44100 # Hz
import wave, struct, math
from numpy import fft
FFT_LENGTH = 2048
# The higher the better, but this seems enough
OVERLAP = 512
FFT_SAMPLE = FFT_LENGTH - OVERLAP
# Nyquist rate = bandwidth available for a given sample rate
NYQUIST_RATE = SAMPLE_RATE / 2.0
# Convert frequencies from Hz to our digital sampling units
LOWPASS /= (NYQUIST_RATE / (FFT_LENGTH / 2.0))
HIGHPASS /= (NYQUIST_RATE / (FFT_LENGTH / 2.0))
zeros = [ 0 for x in range(0, OVERLAP) ]
# Builds filter mask. Note that this filter is BAD, a
# good filter must have a smooth curve, respecting a
# dB/octave attenuation ramp!
mask = []
for f in range(0, FFT_LENGTH / 2 + 1):
rampdown = 1.0
if f > LOWPASS:
rampdown = 0.0
elif f < HIGHPASS:
rampdown = 0.0
mask.append(rampdown)
def bound(sample):
# takes care of "burnt" samples, due some mistake in filter design
if sample > 1.0:
print "!",
sample = 1.0
elif sample < -1.0:
print "!",
sample = -1.0
return sample
original = wave.open("NubiaCantaDalva.wav", "r")
filtered = wave.open("NubiaFiltered.wav", "w")
filtered.setnchannels(1)
filtered.setsampwidth(2)
filtered.setframerate(SAMPLE_RATE)
n = original.getnframes()
original = struct.unpack('%dh' % n, original.readframes(n))
# scale from 16-bit signed WAV to float
original = [s / 32768.0 for s in original]
saved_td = zeros
for pos in range(0, len(original), FFT_SAMPLE):
time_sample = original[pos : pos + FFT_LENGTH]
# convert frame to frequency domain representation
frequency_domain = fft.fft(time_sample, FFT_LENGTH)
l = len(frequency_domain)
# mask positive frequencies (f[0] is DC component)
for f in range(0, l/2+1):
frequency_domain[f] *= mask[f]
# mask negative frequencies
# http://stackoverflow.com/questions/933088/clipping-fft-matrix
for f in range(l-1, l/2, -1):
cf = l - f
frequency_domain[f] *= mask[cf]
# convert frame back to time domain
time_domain = fft.ifft(frequency_domain)
# modified overlap-add logic: previously saved samples
# prevail in the beginning, and they are ramped down
# in favor of this frame's samples, to avoid 'clicks'
# in either end of frame.
for i in range(0, OVERLAP):
time_domain[i] *= (i + 0.0) / OVERLAP
time_domain[i] += saved_td[i] * (1.0 - (i + 0.00) / OVERLAP)
# reserve last samples for the next frame
saved_td = time_domain[FFT_SAMPLE:]
# do not write reserved samples right now
time_domain = time_domain[:FFT_SAMPLE]
# scale back into WAV 16-bit and write
time_domain = [ bound(sample) * 32767 for sample in time_domain ]
filtered.writeframes(struct.pack('%dh' % len(time_domain),
*time_domain))
Existe uma pegadinha neste código que merece ser mencionada. Toda filtragem, mesmo que executada no domínio da freqüência, causa distorção de fase. A filtragem DFT trabalha o áudio em blocos fixos, mas você não pode simplesmente reproduzir um bloco após o outro. A descontinuidade é audível como estalos, bem desagradáveis. (Se você não acredita, mude a variável OVERLAP para zero e veja o que acontece.)
Ambas as soluções (overlap-save e overlap-add) que encontrei em livros especializados não pareceram resolver o problema. Modifiquei o método overlap-add, o que aparentemente resolveu o caso. Minha solução faz fade-in do bloco "n+1" ao mesmo tempo que faz fade-out do bloco "n", juntando os dois de forma suave. Talvez haja solução melhor, mas esta gerou áudio de boa qualidade e não encontrei outra forma na literatura.
Outra coisa que deve ser lembrada por qualquer um que vá brincar com FFT, é que a representação do domínio da freqüência é linear. Em geral, estamos acostumados a uma representação logarítmica da mesma, porque corresponde às notas musicais. Isto significa que a resolução de freqüência de uma matriz gerada pelo FFT é (musicalmente) muito mais precisa nas oitavas baixas.
O FFT gera dois "lobos" de freqüência: positivas e negativas. Isto parece muito estranho, já que não existem freqüências negativas no mundo real. Porém o FFT tem esse cacoete matemático. Em nosso filtro, nós tínhamos duas escolhas: simplesmente zerar as freqüências negativas (o que causaria uma redução de volume do áudio final em 50%), ou então filtrá-las da mesma forma que as freqüências positivas (que é o que fazemos).
A matriz gerada pelo FFT consiste de números complexos. Isto soa assustador, mas tem um propósito válido: um número complexo pode representar amplitude e fase de forma compacta. Por exemplo, se encontramos uma amostra 1+0j em 3000Hz, isto significa que o áudio original tem um componente de 3000Hz, amplitude=1 e em fase com cos(x). Se a amostra fosse 0.707+0.707j, isto significaria que o componente ainda tem amplitude=1 (módulo do número complexo) porém a fase é de 45 graus em relação a cos(x).
Algumas aplicações, como o bonito analisador de espectro da figura inicial, não precisam da fase do sinal, só estamos interessados na energia contida em cada freqüência. Neste caso, podemos calcular o FFT, considerar apenas o lobo positivo, e obter os módulos ou valores absolutos das amostras complexas (lembrando que o valor absoluto de um número complexo é sempre real e não-negativo).
Um sinal real, como um áudio, sinal de rádio ou qualquer sinal do mundo real, possui transformada de Fourier com amplitude par. Isto quer dizer que um sinal com freqüência "f", processado pelo FFT, aparece tanto no lobo positivo quanto no lobo negativo, nas freqüências +f e −f.
(A propósito, função "par" possui a propriedade f(x)=f(−x).)
Por outro lado, um sinal real possui transformada com fase ímpar, ou seja, a fase do lobo negativo é oposta à do lado positivo. (Função "impar" é aquela em que f(x)=−f(−x)).
A razão matemática disto é que um sinal real pode ser construído por dois sinais complexos conjugados. Direto de um livro qualquer de cálculo, temos a seguinte relação entre trigonometria e exponenciais:
cos(2πft+θ) = (ei(2πft+θ) + e−i(2πft+θ)) / 2
onde "f" é a freqüência θ é a fase.
Fourier provou que todo sinal pode ser construído por uma soma de senóides e cossenóides, então a relação acima nos permite dizer que qualquer sinal pode ser construído pela soma de pares de exponenciais complexos (cada par com seus próprios valores de "w" e θ).
A transformada de Fourier de
ei(2πft+θ)
é zero para todo valor de F(x) exceto quando x=f. Neste ponto singular F(f), o valor é um número complexo de fase θ e módulo extremamente alto (um impulso unitário) por ser um tom puro.
A transformada da parcela conjugada e−i(2πft+θ) é zero exceto em F(-f). Neste ponto singular, o módulo também é um impulso unitário mas a fase θ é negativa.
Ou seja, a amplitude dos dois lobos da transformada de Fourier é igual (função par) porém a fase é oposta (função ímpar). Esta é a "cara" da transformada de um sinal real, conforme queríamos demonstrar.
Um sinal complexo também pode ser construído por uma soma de exponenciais na forma ei(2πft+θ) porém neste caso não é garantido que virão em pares conjugados. Portanto, a respectiva transformada não tem nenhuma simetria obrigatória entre lobo positivo e negativo.
Se desejamos que uma transformada inversa de Fourier produza um sinal real, precisamos cuidar que o espectro tenha essas características (amplitude par e fase ímpar), do contrário o sinal produzido será complexo, sem significado no mundo real. Sempre podemos converter números complexos na marra, usando apenas a parte real ou obtendo os valores absolutos, mas isso causará algum tipo de distorção, que pode ou não destruir o sinal.
Sinais complexos não podem existir no mundo material, mas podem existir em nossa imaginação ou na memória de um computador. E eles podem ser bastante úteis. Por exemplo, a demodulação I/Q, utilizada em modems e rádios, produz um sinal complexo que representa a freqüência e a fase do sinal original, que pode ser filtrado e transformado como se fosse um sinal real, pois no fim do dia o sinal complexo é apenas um par de sinais reais que correm em paralelo.