Filtro pasa banda en MATLAB: ejemplo 2

Este es otro ejemplo de filtro pasa banda FIR en MATLAB diseñado a partir de su respuesta al impulso ideal (resta de dos funciones sinc), esta vez con una banda muy angosta. El código original cargaba una señal desde un archivo 60000.mat que no está disponible, y tenía un truco en la gráfica que hacía aparecer el filtro en una frecuencia equivocada. En este post te explico cómo calcular realmente la banda de paso, qué estaba mal en la gráfica y te dejo una versión corregida que puedes ejecutar sin archivos externos.

Los parámetros del filtro

El código original usa:

Frecuencia de muestreo: fs = 8000 Hz (en el código, la variable T = 8000 se asigna a fs).
ω_l = 0,745π y ω_h = 0,755π.
L = 128 (129 coeficientes).
N = 60 000 muestras (7,5 segundos de señal).

Para pasar de frecuencia digital a hercios se usa:

f = (ω / π) · (fs / 2)

Con fs/2 = 4000 Hz:

f_l = 0,745 × 4000 = 2980 Hz
f_h = 0,755 × 4000 = 3020 Hz

Es decir, un filtro centrado en 3000 Hz con solo 40 Hz de ancho de banda. (El comentario "50*wh = fh" del código original no corresponde a esta frecuencia de muestreo.)

¿Qué estaba mal en la gráfica?

En el original, la respuesta del filtro se calculaba con fft(bn, N/50) (1200 puntos) y se graficaba sobre un eje (1:N/50)*fs/N. Ese eje usa la resolución de una FFT de 60 000 puntos para una FFT de 1200 puntos, así que todas las frecuencias quedan divididas entre 50. Por eso el filtro aparecía alrededor de 60 Hz (3000 / 50) dentro del xlim([40 80]), cuando en realidad está en 3000 Hz. El eje correcto para una FFT de M puntos es (0:M-1)*fs/M.

Además, la escala H = H*90 era un ajuste manual para que las curvas se vieran parecidas, y las etiquetas "EEG" venían de otro ejercicio.

El problema de la banda angosta

Un filtro FIR de L+1 coeficientes tiene una resolución en frecuencia del orden de fs / L. Con 129 coeficientes y fs = 8000 Hz, eso es unos 62 Hz: más que el ancho de banda deseado de 40 Hz. El filtro resultante deja pasar una banda bastante más ancha y con mucho rizado. Para una banda de 40 Hz conviene usar del orden de 1000 coeficientes o más y aplicar una ventana.

Código corregido

clc; close all; clear all;

fs = 8000;                      % frecuencia de muestreo (Hz)
N  = 60000;                     % 7,5 s de señal
t  = (0:N-1)/fs;
% Señal de prueba: 3000 Hz (deseada) + 1000 Hz (a eliminar) + ruido
data = sin(2*pi*3000*t) + sin(2*pi*1000*t) + 0.3*randn(1,N);

% --- Filtro pasa banda 2980-3020 Hz ---
wl = 0.745*pi;                  % 0.745*4000 = 2980 Hz
wh = 0.755*pi;                  % 0.755*4000 = 3020 Hz
L  = 1024;                      % orden alto para banda angosta
n  = -L/2:L/2;                  % índices simétricos
bn = (sin(wh*n) - sin(wl*n)) ./ (pi*n);
bn(n==0) = (wh - wl)/pi;        % valor central
bn = bn .* blackman(L+1)';      % ventana para reducir lóbulos

% --- Respuesta en frecuencia con eje correcto ---
M = 2^16;
H = abs(fft(bn, M));
f = (0:M/2-1)*fs/M;
figure; plot(f, 20*log10(H(1:M/2) + eps)); grid on;
xlim([2800 3200]); ylim([-100 5]);
xlabel('Frecuencia (Hz)'); ylabel('Magnitud (dB)');
title('Filtro pasa banda 2980-3020 Hz');

% --- Espectro de la señal antes y después ---
y = filter(bn, 1, data);
Fx = abs(fft(data)); Fy = abs(fft(y));
fx = (0:N/2-1)*fs/N;
figure;
plot(fx, Fx(1:N/2), 'k'); hold on; plot(fx, Fy(1:N/2), 'b');
xlabel('Frecuencia (Hz)'); ylabel('|X(f)|');
legend('Original', 'Filtrada');

% --- Señales en el tiempo ---
figure;
subplot(2,1,1); plot(t(1:400), data(1:400), 'k'); grid on;
ylabel('Original');
subplot(2,1,2); plot(t(L+1:L+400), y(L+1:L+400), 'b'); grid on;
xlabel('Tiempo (s)'); ylabel('Filtrada');

Qué hace cada parte

Señal de prueba: en lugar del archivo .mat, se genera una mezcla de 3000 Hz, 1000 Hz y ruido, para comprobar que el filtro deja pasar solo los 3000 Hz.
Diseño vectorizado: con n = -L/2:L/2 el filtro queda simétrico (fase lineal) sin necesidad del ciclo for.
Ventana Blackman: reduce los lóbulos laterales, que con la ventana rectangular original eran altos.
Respuesta en dB: permite ver cuánto se atenúa fuera de la banda.
filter: aplica el filtro y devuelve una señal del mismo largo que la entrada. En la gráfica temporal se saltan las primeras L muestras, que corresponden al transitorio y al retardo de L/2 muestras.

Comparación con fir1

El mismo filtro se puede diseñar con la función de MATLAB:

b = fir1(1024, [0.745 0.755], 'bandpass', blackman(1025));
freqz(b, 1, 2^16, fs);

En fir1, las frecuencias se normalizan respecto a fs/2, igual que nuestras ω/π.

Errores comunes

Mezclar resoluciones de FFT distintas en un mismo eje de frecuencia.
Esperar un filtro muy selectivo con pocos coeficientes.
Escalar curvas a mano para que "coincidan", ocultando errores de eje.
Olvidar que conv devuelve N + L muestras y filter devuelve N.

Preguntas frecuentes

¿Cómo diseño un filtro alrededor de 60 Hz?

Calcula ω = f / (fs/2) × π. Con fs = 8000 Hz, 55–65 Hz da ω entre 0,01375π y 0,01625π; una banda tan baja y angosta necesita miles de coeficientes, o conviene bajar la frecuencia de muestreo primero.

¿Hay ejemplos más sencillos?

Sí: la señal senoidal con filtro pasa banda y el filtro pasa banda de 200 Hz.

¿Funciona en Octave?

Sí, cargando el paquete signal para blackman, fir1 y freqz.

Conclusión

Con fs = 8000 Hz, los valores 0,745π y 0,755π definen una banda de 2980 a 3020 Hz, no de 60 Hz como sugería la gráfica original. Usando un eje de frecuencia correcto, más coeficientes y una ventana, el filtro funciona como se espera.