Pular para conteúdo

Lista 17 · Frequências e a FFT

Acompanha o capítulo 17. A transformada de Fourier com laços e com a FFT, a frequência mais forte de um sinal, a resolução, o limite de Nyquist e o aliasing.

Abrir no Colab

Baixar o caderno

Como esta lista funciona

Abra o caderno no Colab, escreva a sua solução no lugar do # TODO e rode a célula de teste logo abaixo. Ela mesma diz o que passou e o que falhou — você não precisa esperar correção de ninguém.

O código de cada exercício também está aqui embaixo, caso você prefira colar num editor local. Os testes comparam números com tolerância.

Exercícios marcados com ✏️ são à mão, no papel, com 4 casas decimais — como na parte em papel das avaliações.

Exercício 01

✏️ À mão. Um sinal tem 4 amostras, \(x = [1, 0, -1, 0]\), gravadas a 4 amostras por segundo (1 segundo de sinal).

a) Calcule \(C_k = \sum_n x_n \cos(2\pi k n / 4)\), \(S_k = \sum_n x_n \sin(2\pi k n / 4)\) e a amplitude \(A_k = \frac{2}{N}\sqrt{C_k^2 + S_k^2}\) para \(k = 0, 1, 2\). Que sinal é esse?

b) Um sensor grava 3 segundos a 1000 amostras por segundo. Qual a resolução em frequência? Qual a maior frequência que ele enxerga?

c) Um seno de 70 Hz é amostrado a 100 amostras por segundo. Em que frequência ele aparece na FFT?

Exercício 02

A sonda. Escreva amplitude_em(sinal, k), que calcula, com um laço, as somas \(C_k\) e \(S_k\) do capítulo e devolve a amplitude \(A_k = \frac{2}{N}\sqrt{C_k^2 + S_k^2}\).

amplitude_em([1, 0, -1, 0], 1)    # -> 1.0  (o ex01)

import numpy as np


def amplitude_em(sinal, k):
    """Amplitude da frequência k (em voltas por sinal inteiro), pela transformada com laço."""
    N = len(sinal)
    soma_cos = 0.0
    soma_sen = 0.0
    # TODO: para cada amostra n, acumule sinal[n] * cos(2 pi k n / N) e o mesmo com sin
    return 2 * np.sqrt(soma_cos ** 2 + soma_sen ** 2) / N


# --- testes (não altere) ---
import numpy as np


assert abs(amplitude_em([1, 0, -1, 0], 1) - 1) < 1e-12
assert abs(amplitude_em([1, 0, -1, 0], 2)) < 1e-12
t = np.linspace(0, 0.99, 100)
sinal = 3 * np.sin(2 * np.pi * 7 * t)
assert abs(amplitude_em(sinal, 7) - 3) < 1e-9, "um seno de amplitude 3 em 7 Hz"
assert amplitude_em(sinal, 8) < 1e-9
print("Exercício 02: todos os testes passaram!")

Exercício 03

Escreva amplitudes(sinal), que devolve o array das amplitudes de todas as frequências com a FFT: 2 * np.abs(np.fft.rfft(sinal)) / N. Confira que o resultado na posição 1 de [1, 0, -1, 0] é o mesmo do ex02.

import numpy as np


def amplitudes(sinal):
    """Amplitude de cada frequência, de 0 até fs/2, pela FFT."""
    # TODO: rfft, abs e o fator 2 / N
    pass


# --- testes (não altere) ---
import numpy as np


a = amplitudes([1, 0, -1, 0])
assert a is not None, "a função não devolveu nada"
assert len(a) == 3, "4 amostras: frequências 0, 1 e 2"
assert abs(a[1] - 1) < 1e-12 and abs(a[0]) < 1e-12 and abs(a[2]) < 1e-12
t = np.linspace(0, 0.99, 100)
a = amplitudes(np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 12 * t))
assert abs(a[5] - 1) < 1e-9 and abs(a[12] - 0.5) < 1e-9
print("Exercício 03: todos os testes passaram!")

Exercício 04

Um afinador. Escreva frequencia_dominante(sinal, fs), que devolve a frequência, em Hz, do pico mais forte do espectro. Antes da FFT, subtraia a média do sinal (senão a frequência 0 pode ganhar). Use np.fft.rfftfreq(N, 1 / fs) e np.argmax.

import numpy as np


def frequencia_dominante(sinal, fs):
    """A frequência (Hz) de maior amplitude, sem contar a média."""
    # TODO: tire a média, calcule o espectro e as frequências, devolva a do pico
    pass


# --- testes (não altere) ---
import numpy as np


fs = 1000
t = np.linspace(0, 0.999, 1000)
sinal = 5 + np.sin(2 * np.pi * 20 * t) + 0.3 * np.sin(2 * np.pi * 50 * t)
f = frequencia_dominante(sinal, fs)
assert f is not None, "a função não devolveu nada"
assert f != 0, "a média ganhou: subtraia np.mean(sinal) antes da FFT"
assert abs(f - 20) < 1e-9, f"veio {f}"
t = np.linspace(0, 1.998, 1000)
assert abs(frequencia_dominante(np.cos(2 * np.pi * 110 * t), 500) - 110) < 1e-9, "fs = 500: confira o rfftfreq"
print("Exercício 04: todos os testes passaram!")

Exercício 05

Escreva resolucao_e_nyquist(N, fs), que devolve a tupla (resolução, maior frequência visível) para \(N\) amostras gravadas a \(f_s\) amostras por segundo. A duração é \(T = N/f_s\).

resolucao_e_nyquist(4000, 8000)    # -> (2.0, 4000.0)  (o violão do capítulo)

Quantas amostras a 8000 por segundo seriam precisas para uma resolução de 0,5 Hz?

def resolucao_e_nyquist(N, fs):
    """(resolução em Hz, maior frequência visível em Hz) de N amostras a fs por segundo."""
    # TODO: a duração T, a resolução 1 / T e o limite fs / 2
    pass


# --- testes (não altere) ---
r = resolucao_e_nyquist(4000, 8000)
assert r is not None, "a função não devolveu nada"
assert abs(r[0] - 2) < 1e-12 and abs(r[1] - 4000) < 1e-12, f"veio {r}"
r = resolucao_e_nyquist(3000, 1000)
assert abs(r[0] - 1 / 3) < 1e-12 and abs(r[1] - 500) < 1e-12, f"veio {r}"
print("Exercício 05: todos os testes passaram!")

Exercício 06

Aliasing. Um seno de frequência \(f\) amostrado a \(f_s\) aparece na FFT numa frequência entre 0 e \(f_s/2\). A regra: tire de \(f\) quantos \(f_s\) couberem (o resto f % fs); se o resto passar de \(f_s/2\), a frequência aparente é \(f_s\) menos o resto. Escreva frequencia_aparente(f, fs).

frequencia_aparente(12, 15)     # -> 3   (o exemplo do capítulo)
frequencia_aparente(70, 100)    # -> 30  (o ex01)

Confira a regra com a FFT: amostre np.sin(2 * np.pi * 70 * t) a 100 por segundo e procure o pico.

def frequencia_aparente(f, fs):
    """Em que frequência um seno de f Hz aparece, amostrado a fs por segundo."""
    # TODO: o resto de f por fs e, se passar de fs / 2, fs menos o resto
    pass


# --- testes (não altere) ---
assert frequencia_aparente(12, 15) == 3
assert frequencia_aparente(70, 100) == 30
assert frequencia_aparente(20, 100) == 20, "abaixo de fs/2 não muda"
assert frequencia_aparente(260, 100) == 40
assert frequencia_aparente(100, 100) == 0, "amostrado na própria frequência: parece parado"
print("Exercício 06: todos os testes passaram!")

Exercício 07

Biologia — o relógio interno. A temperatura do corpo segue um ritmo circadiano: ela muda ao longo do dia mesmo sem luz. Num experimento, um voluntário ficou 10 dias numa sala sem janelas nem relógio, e a sua temperatura foi medida de hora em hora (240 valores). O ritmo continuou com 24 horas?

Escreva periodo_dominante(serie), que devolve o período, em número de amostras, do ciclo mais forte: tire a média, calcule a FFT com uma amostra por unidade de tempo (np.fft.rfftfreq(N, 1)) e devolva \(1/f\) do pico. O teste usa os dados do voluntário.

import numpy as np


def periodo_dominante(serie):
    """Período (em número de amostras) do ciclo mais forte de uma série."""
    # TODO: tire a média, a FFT, as frequências (uma amostra por unidade) e 1 / f do pico
    pass


# --- testes (não altere) ---
import numpy as np


temperatura = np.array([36.81, 36.88, 36.81, 36.70, 36.77, 36.49, 36.74, 36.73, 36.69, 36.86, 37.02, 36.90, 37.07, 37.38, 37.43, 37.54, 37.55, 37.43, 37.07, 37.36, 36.87, 37.07, 36.98, 36.99, 36.94, 36.74, 36.79, 36.35, 36.62, 36.29, 36.83, 36.71, 36.96, 37.19, 36.68, 36.99, 37.19, 37.35, 37.55, 37.55, 37.58, 37.47, 37.39, 37.22, 37.25, 37.06, 37.10, 37.21, 37.08, 36.74, 36.64, 36.65, 36.40, 36.27, 36.50, 36.87, 36.87, 36.91, 37.11, 36.88, 37.44, 37.22, 37.30, 37.76, 37.34, 37.38, 37.33, 37.07, 37.05, 37.12, 37.03, 36.96, 37.05, 37.07, 36.89, 36.23, 36.71, 36.66, 36.71, 36.64, 36.83, 36.94, 36.79, 37.22, 37.23, 37.21, 37.20, 37.66, 37.60, 37.62, 37.28, 37.57, 37.08, 37.33, 36.70, 37.01, 36.91, 36.41, 36.56, 36.36, 36.29, 36.43, 36.48, 36.84, 36.97, 37.01, 36.88, 36.96, 37.35, 37.16, 37.34, 37.37, 37.58, 37.59, 37.47, 37.28, 37.21, 37.00, 37.09, 36.82, 36.92, 36.92, 36.50, 36.66, 36.55, 36.36, 36.40, 36.54, 36.74, 37.02, 37.04, 37.22, 37.44, 37.14, 37.04, 37.23, 37.20, 37.35, 37.27, 37.40, 36.87, 37.00, 36.96, 36.95, 36.81, 37.10, 36.63, 36.55, 36.61, 36.58, 36.92, 36.61, 36.87, 36.89, 36.99, 37.17, 37.01, 36.92, 37.36, 37.30, 37.59, 37.49, 37.66, 37.29, 36.98, 37.09, 37.16, 36.98, 36.89, 36.75, 36.59, 36.40, 36.68, 36.60, 36.44, 36.95, 36.90, 37.00, 37.02, 37.21, 36.77, 37.21, 37.43, 37.44, 37.25, 37.27, 37.38, 37.43, 36.81, 36.95, 37.05, 37.00, 36.98, 36.97, 36.77, 36.33, 36.46, 36.38, 36.66, 37.02, 37.06, 37.32, 37.21, 37.20, 37.26, 37.12, 37.33, 37.55, 37.62, 37.56, 37.21, 37.06, 36.94, 37.32, 36.99, 37.10, 36.90, 36.77, 36.63, 36.49, 36.52, 36.67, 36.47, 36.80, 36.86, 36.91, 37.11, 37.28, 36.91, 37.24, 37.52, 37.59, 37.26, 37.54, 37.29, 37.15, 37.22, 36.89, 37.04, 37.10])
p = periodo_dominante(temperatura)
assert p is not None, "a função não devolveu nada"
assert abs(p - 24.0) < 1e-6, f"veio {p}"
t = np.linspace(0, 119, 120)
assert abs(periodo_dominante(np.sin(2 * np.pi * t / 12)) - 12) < 1e-6
print("Exercício 07: todos os testes passaram!")