Filtro de áudio DFT

Índice || Index

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.

Figura 1: Analisador de espectro do iTunes

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:

(link to audio)

Agora, a versão filtrada (filtro passa-banda 1000 a 4000Hz), que soa como aqueles radinhos de pilha:

(link to audio)

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. Neste caso, podemos utilizar DCT - Transformada Discreta de Cosseno, cujo resultado é uma matriz de números reais. A forma de onda original ainda pode ser obtida pela transformada inversa (IDCT) mas neste caso o sinal reconstituído terá fase diferente do original.