Filtro Digital com Séries de Taylor e Cepstrum para Detecção de Assobio em RP2040

5 — FFT + log-magnitude (com Taylor) + IFFT ⇒ Cepstrum e detecção do assobio (já aciona LED)

Agora vamos fechar o “miolo” do sistema: pegar o frame filtrado, calcular cepstrum, encontrar o pico de quefrência e ligar/desligar o LED.

A partir daqui, o firmware já passa a funcionar como detector (mesmo que ainda possamos melhorar robustez, janelas e calibração nas próximas seções).

5.1 Estruturas complexas e utilitários mínimos

Vamos começar com um tipo complexo simples e algumas operações básicas.

typedef struct {
    float re;
    float im;
} cpx_t;

static inline cpx_t cpx_add(cpx_t a, cpx_t b) { return (cpx_t){a.re + b.re, a.im + b.im}; }
static inline cpx_t cpx_sub(cpx_t a, cpx_t b) { return (cpx_t){a.re - b.re, a.im - b.im}; }
static inline cpx_t cpx_mul(cpx_t a, cpx_t b) {
    return (cpx_t){ a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re };
}

typedef struct {
    float re;
    float im;
} cpx_t;

static inline cpx_t cpx_add(cpx_t a, cpx_t b) { return (cpx_t){a.re + b.re, a.im + b.im}; }
static inline cpx_t cpx_sub(cpx_t a, cpx_t b) { return (cpx_t){a.re - b.re, a.im - b.im}; }
static inline cpx_t cpx_mul(cpx_t a, cpx_t b) {
    return (cpx_t){ a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re };
}

5.2 FFT radix-2 iterativa (in-place) para N=256

5.2.1 Por que radix-2 e N=256

256 = \(2^8\), ideal para FFT radix-2.
Leve o suficiente para RP2040 em tempo real com (f_s=8kHz).
Frame curto ⇒ baixa latência.

5.2.2 Tabela de twiddles (pré-cálculo)

Vamos pré-calcular \(W_N^k = e^{-j 2\pi k/N}\) uma vez.

#include <math.h>
#define FRAME_N 256

static cpx_t g_tw[FRAME_N/2];

static void fft_init_twiddles(void)
{
    for (int k = 0; k < FRAME_N/2; k++) {
        float ang = -2.0f * (float)M_PI * (float)k / (float)FRAME_N;
        g_tw[k].re = cosf(ang);
        g_tw[k].im = sinf(ang);
    }
}

#include <math.h>
#define FRAME_N 256

static cpx_t g_tw[FRAME_N/2];

static void fft_init_twiddles(void)
{
    for (int k = 0; k < FRAME_N/2; k++) {
        float ang = -2.0f * (float)M_PI * (float)k / (float)FRAME_N;
        g_tw[k].re = cosf(ang);
        g_tw[k].im = sinf(ang);
    }
}

5.2.3 Bit-reversal (reordenação)

FFT iterativa precisa que a entrada esteja em ordem bit-reversed.

static inline uint32_t bitrev8(uint32_t x)
{
    // reverte 8 bits (pois N=256)
    x = ((x & 0xF0u) >> 4) | ((x & 0x0Fu) << 4);
    x = ((x & 0xCCu) >> 2) | ((x & 0x33u) << 2);
    x = ((x & 0xAAu) >> 1) | ((x & 0x55u) << 1);
    return x;
}

static void fft_bitreverse(cpx_t *a)
{
    for (uint32_t i = 0; i < FRAME_N; i++) {
        uint32_t j = bitrev8(i);
        if (j > i) {
            cpx_t tmp = a[i];
            a[i] = a[j];
            a[j] = tmp;
        }
    }
}

static inline uint32_t bitrev8(uint32_t x)
{
    // reverte 8 bits (pois N=256)
    x = ((x & 0xF0u) >> 4) | ((x & 0x0Fu) << 4);
    x = ((x & 0xCCu) >> 2) | ((x & 0x33u) << 2);
    x = ((x & 0xAAu) >> 1) | ((x & 0x55u) << 1);
    return x;
}

static void fft_bitreverse(cpx_t *a)
{
    for (uint32_t i = 0; i < FRAME_N; i++) {
        uint32_t j = bitrev8(i);
        if (j > i) {
            cpx_t tmp = a[i];
            a[i] = a[j];
            a[j] = tmp;
        }
    }
}

5.2.4 FFT direta (forward)

static void fft_forward(cpx_t *a)
{
    fft_bitreverse(a);

    for (uint32_t len = 2; len <= FRAME_N; len <<= 1) {
        uint32_t half = len >> 1;
        uint32_t step = (FRAME_N / len);

        for (uint32_t i = 0; i < FRAME_N; i += len) {
            uint32_t tw = 0;
            for (uint32_t j = 0; j < half; j++) {
                cpx_t u = a[i + j];
                cpx_t v = cpx_mul(a[i + j + half], g_tw[tw]);

                a[i + j]        = cpx_add(u, v);
                a[i + j + half] = cpx_sub(u, v);

                tw += step;
            }
        }
    }
}

static void fft_forward(cpx_t *a)
{
    fft_bitreverse(a);

    for (uint32_t len = 2; len <= FRAME_N; len <<= 1) {
        uint32_t half = len >> 1;
        uint32_t step = (FRAME_N / len);

        for (uint32_t i = 0; i < FRAME_N; i += len) {
            uint32_t tw = 0;
            for (uint32_t j = 0; j < half; j++) {
                cpx_t u = a[i + j];
                cpx_t v = cpx_mul(a[i + j + half], g_tw[tw]);

                a[i + j]        = cpx_add(u, v);
                a[i + j + half] = cpx_sub(u, v);

                tw += step;
            }
        }
    }
}

5.2.5 IFFT (inverse)

IFFT pode ser feita por:

Conjugar entrada
FFT forward
Conjugar saída e dividir por N

static void fft_inverse(cpx_t *a)
{
    // conjuga
    for (int i = 0; i < FRAME_N; i++) a[i].im = -a[i].im;

    fft_forward(a);

    // conjuga e normaliza
    for (int i = 0; i < FRAME_N; i++) {
        a[i].re =  a[i].re / (float)FRAME_N;
        a[i].im = -a[i].im / (float)FRAME_N;
    }
}

static void fft_inverse(cpx_t *a)
{
    // conjuga
    for (int i = 0; i < FRAME_N; i++) a[i].im = -a[i].im;

    fft_forward(a);

    // conjuga e normaliza
    for (int i = 0; i < FRAME_N; i++) {
        a[i].re =  a[i].re / (float)FRAME_N;
        a[i].im = -a[i].im / (float)FRAME_N;
    }
}

5.3 Montando o espectro: |X[k]| e log(|X[k]|) usando ln_approx (Taylor)

Nós já criamos ln_approx() na seção anterior, usando Séries de Taylor para aproximar o log.

Agora vamos usar isso no cepstrum:

Magnitude:
\[
|X[k]| = \sqrt{Re^2 + Im^2}
\]
Log-magnitude:
\[
L[k] = \ln(|X[k]| + \epsilon)
\]

No cepstrum real, a fase não importa. Fazemos IFFT do log-magnitude real.

Detalhe importante: em vez de usar sqrtf(), podemos usar (Re^2+Im^2) e fazer ln(sqrt(p)) = 0.5 ln(p). Assim economizamos custo.

static void build_log_magnitude(cpx_t *spec)
{
    const float eps = 1e-12f;

    for (int k = 0; k < FRAME_N; k++) {
        float p = spec[k].re*spec[k].re + spec[k].im*spec[k].im + eps;
        float ln_p = ln_approx(p);        // Taylor-based (seção 4)
        float ln_mag = 0.5f * ln_p;

        spec[k].re = ln_mag;
        spec[k].im = 0.0f;
    }
}

static void build_log_magnitude(cpx_t *spec)
{
    const float eps = 1e-12f;

    for (int k = 0; k < FRAME_N; k++) {
        float p = spec[k].re*spec[k].re + spec[k].im*spec[k].im + eps;
        float ln_p = ln_approx(p);        // Taylor-based (seção 4)
        float ln_mag = 0.5f * ln_p;

        spec[k].re = ln_mag;
        spec[k].im = 0.0f;
    }
}

5.4 Cepstrum: IFFT do log-magnitude e busca do pico de quefrência

Agora:

FFT no frame
log-magnitude
IFFT
procurar pico (q) numa faixa coerente

Como discutido:

Para (f_s=8000), assobio (800–3000 Hz) dá (q \approx 3) a (10).
Vamos procurar numa faixa um pouco maior para robustez: q_min=3, q_max=20.

static float compute_cepstrum_and_detect(const float *frame, int *best_q_out)
{
    static cpx_t buf[FRAME_N];

    // 1) carrega frame (real) em complexo
    for (int i = 0; i < FRAME_N; i++) {
        buf[i].re = frame[i];
        buf[i].im = 0.0f;
    }

    // 2) FFT
    fft_forward(buf);

    // 3) log-magnitude
    build_log_magnitude(buf);

    // 4) IFFT => cepstrum (real em buf[n].re)
    fft_inverse(buf);

    // 5) busca pico na faixa de quefrência
    const int q_min = 3;
    const int q_max = 20;

    float best = -1e9f;
    int best_q = q_min;

    for (int q = q_min; q <= q_max; q++) {
        float v = buf[q].re;  // cepstrum real
        if (v > best) {
            best = v;
            best_q = q;
        }
    }

    if (best_q_out) *best_q_out = best_q;
    return best; // “força” do pico
}

static float compute_cepstrum_and_detect(const float *frame, int *best_q_out)
{
    static cpx_t buf[FRAME_N];

    // 1) carrega frame (real) em complexo
    for (int i = 0; i < FRAME_N; i++) {
        buf[i].re = frame[i];
        buf[i].im = 0.0f;
    }

    // 2) FFT
    fft_forward(buf);

    // 3) log-magnitude
    build_log_magnitude(buf);

    // 4) IFFT => cepstrum (real em buf[n].re)
    fft_inverse(buf);

    // 5) busca pico na faixa de quefrência
    const int q_min = 3;
    const int q_max = 20;

    float best = -1e9f;
    int best_q = q_min;

    for (int q = q_min; q <= q_max; q++) {
        float v = buf[q].re;  // cepstrum real
        if (v > best) {
            best = v;
            best_q = q;
        }
    }

    if (best_q_out) *best_q_out = best_q;
    return best; // “força” do pico
}

5.5 Decisão e acionamento do LED com histerese

Para evitar “pisca-pisca”, usamos:

um limiar para ligar (TH_ON)
um limiar menor para desligar (TH_OFF)

Além disso, estimamos frequência detectada:

\[
f_0 \approx \frac{f_s}{q}
\]

#define FS_HZ 8000

static bool whistle_state = false;

static void whistle_update_led(float peak, int q, uint led_gpio)
{
    // Ajuste inicial: você vai calibrar na prática depois
    const float TH_ON  = 0.35f;
    const float TH_OFF = 0.25f;

    float f0 = (float)FS_HZ / (float)q;

    // sanity: só consideramos dentro do “mundo do assobio”
    bool f_ok = (f0 >= 700.0f && f0 <= 3500.0f);

    if (!whistle_state) {
        if (f_ok && peak > TH_ON) whistle_state = true;
    } else {
        if (!f_ok || peak < TH_OFF) whistle_state = false;
    }

    gpio_put(led_gpio, whistle_state ? 1 : 0);
}

#define FS_HZ 8000

static bool whistle_state = false;

static void whistle_update_led(float peak, int q, uint led_gpio)
{
    // Ajuste inicial: você vai calibrar na prática depois
    const float TH_ON  = 0.35f;
    const float TH_OFF = 0.25f;

    float f0 = (float)FS_HZ / (float)q;

    // sanity: só consideramos dentro do “mundo do assobio”
    bool f_ok = (f0 >= 700.0f && f0 <= 3500.0f);

    if (!whistle_state) {
        if (f_ok && peak > TH_ON) whistle_state = true;
    } else {
        if (!f_ok || peak < TH_OFF) whistle_state = false;
    }

    gpio_put(led_gpio, whistle_state ? 1 : 0);
}

5.6 Integração com o que já temos: processar frame quando encher

Você já tem:

filtro Taylor passa-faixa
frame buffer
ln_approx()

Agora, ao fechar um frame, chamamos cepstrum e atualizamos LED.

static void process_frame_if_ready(frame_buf_t *fb, uint led_gpio)
{
    if (!fb->full) return;
    fb->full = false;

    int q = 0;
    float peak = compute_cepstrum_and_detect(fb->data, &q);

    whistle_update_led(peak, q, led_gpio);
}

static void process_frame_if_ready(frame_buf_t *fb, uint led_gpio)
{
    if (!fb->full) return;
    fb->full = false;

    int q = 0;
    float peak = compute_cepstrum_and_detect(fb->data, &q);

    whistle_update_led(peak, q, led_gpio);
}

5.7 O que você já consegue testar agora

Com esse bloco pronto:

Assobie perto do microfone
O LED deve ligar e permanecer ligado enquanto o assobio é contínuo
Ao parar, ele deve apagar após a queda do pico

Se estiver sensível demais:

Aumente TH_ON e TH_OFF
Se estiver “duro” demais:
Reduza TH_ON e TH_OFF

Na próxima seção, vamos tornar isso bem mais robusto, adicionando:

janela (Hann) no frame
normalização de energia / gate de silêncio
média temporal do pico (suavização)
e então entregar o código completo e funcional da BitDogLab (main + ADC + GPIO + clock) tudo integrado.