Señal senoidal con filtro pasa banda en MATLAB

En este ejemplo generamos una señal senoidal en MATLAB y la pasamos por un filtro pasa banda FIR diseñado "a mano" a partir de su respuesta al impulso ideal. Es un ejercicio clásico de procesamiento digital de señales: permite ver cómo se construye un filtro con funciones sinc, cómo se analiza con la FFT y cómo se aplica con una convolución. En este post te explico el código paso a paso, corrijo algunos detalles del código original y te muestro cómo mejorarlo con una ventana.

Resultado

La figura muestra el espectro de la señal (en negro) y la respuesta en frecuencia del filtro (en azul). El filtro tiene una banda de paso entre 45 Hz y 150 Hz y 129 coeficientes. Como la señal es de 100 Hz, queda dentro de la banda de paso y sale del filtro prácticamente sin atenuación.

La idea: filtro pasa banda ideal

La respuesta al impulso de un filtro pasa bajos ideal con frecuencia de corte ω_c (en radianes por muestra) es:

h[n] = sen(ω_c·n) / (π·n), con h[0] = ω_c/π

Un pasa banda se obtiene restando dos pasa bajos: uno con corte alto ω_h menos otro con corte bajo ω_l. Como la respuesta ideal es infinita, se recorta a L+1 coeficientes centrados en n = 0.

Conversión de frecuencias

Las frecuencias digitales se expresan como fracción de π, donde π corresponde a la mitad de la frecuencia de muestreo (fs/2). Para convertir:

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

Con fs = 1000 Hz (fs/2 = 500 Hz):

• ω_l = 0,09π → f_l = 0,09 × 500 = 45 Hz
• ω_h = 0,3π → f_h = 0,3 × 500 = 150 Hz

En el comentario del código original aparecía "50*wh = fh", que no es correcto; la fórmula de arriba es la que da los 45–150 Hz que se ven en la figura.

Código corregido

clc; close all; clear all;

% --- Señal de prueba ---
fs = 1000;            % frecuencia de muestreo (Hz)
N  = 1000;            % número de muestras
t  = (0:N-1)/fs;      % vector de tiempo
data = sin(2*pi*100*t);   % senoidal de 100 Hz

% --- Diseño del filtro pasa banda ---
wl = 0.09*pi;         % corte bajo: 0.09*500 = 45 Hz
wh = 0.30*pi;         % corte alto: 0.30*500 = 150 Hz
L  = 128;             % orden (L+1 = 129 coeficientes)

bn = zeros(1, L+1);
for i = 1:L+1
    n = i - 1 - L/2;            % n va de -64 a 64 (simétrico)
    if n == 0
        bn(i) = wh/pi - wl/pi;  % valor central
    else
        bn(i) = sin(wh*n)/(pi*n) - sin(wl*n)/(pi*n);
    end
end
% bn = bn .* blackman(L+1)';   % opcional: aplicar ventana

% --- Espectros ---
H_data = abs(fft(data));
freq = (0:N/2-1)*fs/N;          % eje de frecuencia correcto
plot(freq, H_data(1:N/2), 'k'); hold on;
H = abs(fft(bn, N));
H = H * 1.2 * (max(H_data)/max(H));   % escalado solo para graficar
plot(freq, H(1:N/2), 'b');
xlabel('Frecuencia (Hz)'); ylabel('|H(f)|');
legend('Espectro de la señal', 'Filtro (escalado)');

% --- Filtrado por convolución ---
y = conv(data, bn);
figure;
subplot(2,1,1);
plot(t(1:N/2), data(1:N/2), 'k'); grid on;
xlabel('Tiempo (s)'); ylabel('Señal original');
subplot(2,1,2);
plot(t(1:N/2), y(1:N/2), 'k'); grid on;
xlabel('Tiempo (s)'); ylabel('Señal filtrada');

Qué se corrigió y por qué

• Simetría: en el original, n = i - L/2 hacía que n fuera de −63 a 65, un filtro ligeramente asimétrico. Con n = i - 1 - L/2 va de −64 a 64 y el filtro tiene fase lineal.
• Eje de frecuencias: debe empezar en 0 Hz: (0:N/2-1)*fs/N. El original usaba (1:N/2), que desplaza el eje un bin (1 Hz).
• Etiquetas: el original decía "EEG" porque el código venía de un ejercicio con señales de electroencefalografía; aquí la señal es una senoidal de prueba.
• Prealocación: bn = zeros(1, L+1) hace el código más rápido y claro.
• Comentarios: se corrigió la fórmula de conversión de frecuencias.

Paso a paso de lo que hace el código

1. Genera la señal: 1000 muestras de una senoidal de 100 Hz muestreada a 1000 Hz (1 segundo).
2. Diseña el filtro: calcula los 129 coeficientes restando dos funciones sinc.
3. Analiza en frecuencia: con la FFT muestra el espectro de la señal (un pico en 100 Hz) y la respuesta del filtro.
4. Filtra: la función conv aplica el filtro a la señal.
5. Compara: grafica la señal original y la filtrada.

Observación: la señal filtrada aparece con un pequeño "arranque" y retrasada L/2 = 64 muestras (64 ms). Es el retardo normal de un filtro FIR de fase lineal y el transitorio inicial de la convolución.

Mejora: aplicar una ventana

Cortar la respuesta ideal de golpe (ventana rectangular) produce rizado en la banda de paso y lóbulos laterales altos en la banda de rechazo (fenómeno de Gibbs). Si descomentas la línea de blackman, los lóbulos laterales bajan muchísimo, a cambio de una transición un poco más ancha. También puedes diseñar el mismo filtro con la función de MATLAB:

b = fir1(128, [0.09 0.30], 'bandpass', blackman(129));
y = filter(b, 1, data);

En fir1, las frecuencias se normalizan respecto a fs/2, así que 0.09 y 0.30 corresponden a 45 y 150 Hz.

Errores comunes

• Confundir la normalización: algunas funciones normalizan respecto a fs y otras respecto a fs/2.
• Olvidar que conv devuelve N + L muestras, no N.
• Graficar la FFT completa: la segunda mitad es el espejo de la primera para señales reales.

Preguntas frecuentes

¿Qué pasa si la señal es de 300 Hz?

Queda fuera de la banda de paso y sale muy atenuada. Pruébalo cambiando el 100 por 300.

¿Funciona en GNU Octave?

Sí; para fir1 y blackman instala y carga el paquete signal.

¿Hay otro ejemplo parecido?

Sí: el filtro pasa banda de 200 Hz en MATLAB. Y si estás empezando, mira cómo generar una onda senoidal en MATLAB.

Conclusión

Restando dos funciones sinc se obtiene un filtro pasa banda FIR; con 129 coeficientes y fs = 1000 Hz, los valores 0,09π y 0,3π dan una banda de 45 a 150 Hz. Corrigiendo la simetría, el eje de frecuencias y aplicando una ventana, el ejemplo queda listo para usar en clase o en tus proyectos.