FILTRADO DE AUDIO EN MATLAB
Un filtro digital es un sistema que, dependiendo de
las variaciones de las señales de entrada en el tiempo
y amplitud, se realiza un procesamiento matemático
sobre dicha señal; generalmente mediante el uso de
la Transformada rápida de Fourier; obteniéndose en
la salida el resultado del procesamiento matemático
o la señal de salida.
Los filtros digitales tienen como entrada una señal
analógica o digital y en su salida tienen otra señal
analógica o digital, pudiendo haber cambiado en
amplitud, frecuencia o fase dependiendo de las
características del filtro digital.
El filtrado digital es parte del procesado de señal
digital. Se le da la denominación de digital más por
su funcionamiento interno que por su dependencia
del tipo de señal a filtrar, así podríamos llamar filtro
digital tanto a un filtro que realiza el procesado de
señales digitales como a otro que lo haga de señales
analógicas.
Comunmente se usa para atenuar o amplificar
algunas frecuencias, por ejemplo se puede
implementar un sistema para controlar los tonos
graves y agudos del audio del estéreo del auto.
La gran ventaja de los filtros digitales sobre los
analógicos es que presentan una gran estabilidad de
funcionamiento en el tiempo.
El filtrado digital consiste en la realización interna de
un procesado de datos de entrada.
En general el proceso de filtrado consiste en el
muestreo digital de la señal de entrada, el
procesamiento considerando el valor actual de
entrada y considerando las entradas anteriores. El
último paso es la reconstrucción de la señal de salida.
En general la mecánica del procesamiento es:
1. Tomar las muestras actuales y algunas muestras
anteriores (que previamente habían sido
almacenadas) para multiplicadas por unos
coeficientes definidos.
2. También se podría tomar valores de la salida en
instantes pasados y multiplicarlos por otros
coeficientes.
3. Finalmente todos los resultados de todas estas
multiplicaciones son sumados, dando una salida para
el instante actual.
El procesamiento interno y la entrada del filtro serán
digitales, por lo que puede ser necesario una
conversión analógica-digital o digital-analógica para
uso de filtros digitales con señales analógicas.
Un tema muy importante es considerar las
limitaciones del filtro de entrada debido a Teorema
de muestreo de Nyquist-Shannon que en pocas
palabras; si quiero procesar hasta una frecuencia de
10KHz, debo muestrear a por lo menos 20 KHz.
Los filtros digitales se usan frecuentemente para
tratamiento digital de la imagen o para tratamiento
del sonido digital. [1]
Pulsa el botón descargar para bajar los programas:
Advertencia: los scripts reproducen el audio.
Mantener bajo el volumen del PC.
%% Procesamiento de una señal de audio usando MATLAB
%% Crear señal de audio
% Frecuencia fundamental
f0=1e3; % 1KHz
% Amplitud
a=4; % V=4
% Frecuencia de muestreo
fs=44.1e3; % Frecuencia de una señal de audio
% Tiempo de duración en segundos
T=1.5;
L = round(T*fs); % Número de muestras
% Frecuencia normalizada
fn=f0/fs;
y = a*sin(2*pi*fn*(0:L-1))+0.5*a*sin(2*pi*2*fn*(0:L-1));
% Graficar la señal original
subplot(411)
plot((0:L-1)/fs,y)
title('SEÑAL ORIGINAL')% Título
xlabel('Tiempo (s)') % Etiqueta del eje X
ylabel('Amplitud (V)') % Etiqueta del eje Y
xlim([0 10/1000]) % Límite de la señal
%% Grabar y reproducir la señal de audio
%wavwrite(y,fs,'audio')
% wavplay(y,fs)
%% FFT de la señal
subplot(412)
% Llamado a la función que calcula la FFT
fft_signal(y,fs);title('ESPECTRO DE LA SEÑAL ORIGINAL')
xlim([0 2500])
%% Filtrado de la señal
% Frecuencia normalizada
fNorm = 1500 / (fs/2);
% Cálculo de los coeficientes del filtro (filtro pasa bajas)
[b,a] = butter(10, fNorm, 'low');
% Filtrado de la señal
y_Low = filtfilt(b, a, y);
% Graficación de la señal en el tiempo
subplot(413)
plot((0:L-1)/fs,y_Low)
title('SEÑAL FILTRADA')
xlabel('Tiempo (s)')
ylabel('Amplitud (V)')
xlim([0 10/1000])
% Graficación de la señal en frecuencia
subplot(414)
% Llamado a la función que calcula la FFT
fft_signal(y_Low,fs);title('ESPECTRO DE LA SEÑAL FILTRADA')
xlim([0 2500])
%% Gráficas del filtro
% Respuesta en frecuencia del filtro
[H,w]=freqz(b,a,512,1);
figure(2)
%Trazado de la respuesta en Magnitud
subplot(221)
plot(w,20*log10(abs(H)));
grid on;
title ('Filtro pasa-altos, Respuesta en magnitud');
xlabel('frecuencia');
ylabel('H(f) db')
xlim([0 0.4])
% Respuesta en fase
subplot(222)
plot(w,angle(H));
grid on;
title ('Filtro pasa-altos, Respuesta en fase');
xlabel('frecuencia')
ylabel('ángulo de H rad')
xlim([0 0.4])
%Respuesta al impulso
subplot(223)
[y_eje,t]= impz(b,a,60);
stem(t,y_eje);
title ('Filtro pasa-altos, Respuesta al impulso');
%Ploteo de los polos y ceros
z= roots(b); % Ceros
p = roots(a); % Polos
subplot(224)
zplane(z,p)
title('Polos y ceros')
legend('Ceros','Polos')
%% Reproducción de audio de entrada y salida
pause(2)
disp('Audio de entrada')
wavplay(0.5*y,fs)
disp('Audio de salida (señal filtrada)')
wavplay(0.5*y_Low,fs)
%% Procesamiento de una señal de audio usando MATLAB
%% Selección del tipo de filtrado
% 1 -> Pasa bajo
% 2 -> Pasa alto
% 3 -> Pasa banda
tipo=3
%% Crear señal de audio
% Frecuencia fundamental
f0=1e3; % 1KHz
% Amplitud
a=2; % V=4
% Frecuencia de muestreo
fs=44.1e3; % Frecuencia de una señal de audio CD
% Tiempo de duración en segundos
T=2;
% Vector de tiempo
t=linspace(0,T,T*fs);
% Creación de la señal
% Primer señal (tono 1)
s1=a*sin(2*pi*f0*t);
% Segunda señal (tono 2)
s2=0.75*a*sin(2*pi*(1.5*f0)*t);
% Tercera señal (tono 3)
s3=0.5*a*sin(2*pi*(2*f0)*t);
% Señal compuesta (suma de dos tonos)
y = s1 + s2 + s3;
% Graficar la señal original
subplot(411)
plot(t,y)
title('SEÑAL ORIGINAL')% Título
xlabel('Tiempo (s)') % Etiqueta del eje X
ylabel('Amplitud (V)') % Etiqueta del eje Y
xlim([0 20/1000]) % Límite de la señal
%% Grabar y reproducir la señal de audio
% wavwrite(0.1*y,fs,'audio')
% wavplay(0.1*y,fs)
%% FFT de la señal
subplot(412)
% Llamado a la función que calcula la FFT
fft_signal(y,fs);title('ESPECTRO DE LA SEÑAL ORIGINAL')
xlim([0 2500])
%% Filtrado de la señal
switch tipo
case 1
% Cálculo de los coeficientes del filtro (filtro pasa bajas)
% Este filtrado deja solo la señal de 1000 Hz
% Frecuencia normalizada
titutlo='FILTRO PASA BAJAS';
fNorm = 1200 / (fs/2);
[b,a] = butter(10, fNorm, 'low');
case 2
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Cálculo de los coeficientes del filtro (filtro pasa bajas)
% Este filtrado deja solo la señal de 2000 Hz
% Frecuencia normalizada
titutlo='FILTRO PASA ALTAS';
fNorm = 1900 / (fs/2);
[b,a] = butter(10, fNorm, 'high');
otherwise
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Cálculo de los coeficientes del filtro (filtro pasa banda)
% Este filtrado deja solo la señal de 1500 Hz
% Frecuencias normalizadas
titutlo='FILTRO PASA BANDA';
fNorm_1 = 1000 / (fs/2); % 1250
fNorm_2 = 2000 / (fs/2); % 1700
[b_alta,a_alta] = butter(10, fNorm_1, 'high');
[b,a] = butter(10, fNorm_2, 'low');
y_alta=filtfilt(b_alta, a_alta, y);
y=y_alta;
% Wp = [1490 1510]/(fs/2); Ws = [1410 1560]/(fs/2);
% Rp = 3; Rs = 40;
% [n,Wn] = buttord(Wp,Ws,Rp,Rs)
% [b,a] = butter(n,Wn);
% y_Low = filtfilt(b,a,y);
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end
% Filtrado de la señal
y_Low = filtfilt(b, a, y);
% Graficación de la señal en el tiempo
subplot(413)
plot(t,y_Low)
title('SEÑAL FILTRADA')
xlabel('Tiempo (s)')
ylabel('Amplitud (V)')
xlim([0 20/1000])
% Graficación de la señal en frecuencia
subplot(414)
% Llamado a la función que calcula la FFT
fft_signal(y_Low,fs);title('ESPECTRO DE LA SEÑAL FILTRADA')
xlim([0 2500])
%% Gráficas del filtro
% Respuesta en frecuencia del filtro
[H,w]=freqz(b,a,512,fs);
figure(2)
%Trazado de la respuesta en Magnitud
subplot(221)
plot(w,20*log10(abs(H)));
grid on;
title ([titutlo, 'Respuesta en magnitud']);
xlabel('Frecuencia (Hz)');
ylabel('H(f) db')
% xlim([0 0.4])
% Respuesta en fase
subplot(222)
plot(w,angle(H));
grid on;
title ([titutlo,'Respuesta en fase']);
xlabel('Frecuencia (Hz)')
ylabel('ángulo de H rad')
% xlim([0 0.4])
%Respuesta al impulso
subplot(223)
[y_eje,tt]= impz(b,a,60);
stem(tt,y_eje);
title ([titutlo,'Respuesta al impulso']);
%Ploteo de los polos y ceros
z= roots(b); % Ceros
p = roots(a); % Polos
subplot(224)
zplane(z,p)
title('Polos y ceros')
legend('Ceros','Polos')
%% Reproducción de audio de entrada y salida
pause(2)
disp('Audio de entrada')
% Se multiplica por 0.1 para atenuar la salida del tono por la bocina
wavplay(0.1*y,fs)
disp('Audio de salida (señal filtrada)')
% Se multiplica por 0.1 para atenuar la salida del tono por la bocina
wavplay(0.1*y_Low,fs)
%% Procesamiento de una señal de audio usando MATLAB
%% Selección del tipo de filtrado
% 1 -> Pasa bajo
% 2 -> Pasa alto
% 3 -> Pasa banda
tipo=1;
%% Crear señal de audio
% Frecuencia fundamental
f0=8e3; % 8KHz
% Amplitud
a=3; % V=4
% Frecuencia de muestreo
fs=44.1e3; % Frecuencia de una señal de audio CD
% Tiempo de duración en segundos
T=1.5;
% Vector de tiempo
t=linspace(0,T,T*fs);
% Creación de la señal
% Primer señal (tono 1)
s1=a*sin(2*pi*f0*t);
% Segunda señal (tono 2)
s2=0.75*a*sin(2*pi*(1.5*f0)*t);
% Tercera señal (tono 3)
s3=0.5*a*sin(2*pi*(2*f0)*t);
% Señal compuesta (suma de dos tonos)
y = s1 + s2 + s3;
% Graficar la señal original
subplot(411)
plot(t,y)
title('SEÑAL ORIGINAL')% Título
xlabel('Tiempo (s)') % Etiqueta del eje X
ylabel('Amplitud (V)') % Etiqueta del eje Y
xlim([0 20/f0]) % Límite de la señal
%% Grabar y reproducir la señal de audio
% wavwrite(0.1*y,fs,'audio')
% wavplay(0.1*y,fs)
%% FFT de la señal
subplot(412)
% Llamado a la función que calcula la FFT
fft_signal(y,fs);title('ESPECTRO DE LA SEÑAL ORIGINAL')
xlim([0 f0*3])
%% Filtrado de la señal
switch tipo
case 1
% Cálculo de los coeficientes del filtro (filtro pasa bajas)
% Este filtrado deja solo la señal de 1000 Hz
% Frecuencia normalizada
titutlo='FILTRO PASA BAJAS';
fNorm = 8.5e3 / (fs/2);
[b,a] = butter(10, fNorm, 'low');
case 2
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Cálculo de los coeficientes del filtro (filtro pasa bajas)
% Este filtrado deja solo la señal de 2000 Hz
% Frecuencia normalizada
titutlo='FILTRO PASA ALTAS';
fNorm = 15e3 / (fs/2);
[b,a] = butter(10, fNorm, 'high');
otherwise
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Cálculo de los coeficientes del filtro (filtro pasa banda)
% Este filtrado deja solo la señal de 1500 Hz
% Frecuencias normalizadas
titutlo='FILTRO PASA BANDA';
Wp = [11.5e3 12.5e3]/(fs/2); Ws = [11e3 13e3]/(fs/2);
Rp = 3; Rs = 40;
[n,Wn] = buttord(Wp,Ws,Rp,Rs)
[b,a] = butter(n,Wn);
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end
% Filtrado de la señal
y_Low = filtfilt(b, a, y);
% Graficación de la señal en el tiempo
subplot(413)
plot(t,y_Low)
title('SEÑAL FILTRADA')
xlabel('Tiempo (s)')
ylabel('Amplitud (V)')
xlim([0 20/f0])
% Graficación de la señal en frecuencia
subplot(414)
% Llamado a la función que calcula la FFT
fft_signal(y_Low,fs);title('ESPECTRO DE LA SEÑAL FILTRADA')
xlim([0 3*f0])
%% Gráficas del filtro
% Respuesta en frecuencia del filtro
[H,w]=freqz(b,a,512,fs);
figure(2)
%Trazado de la respuesta en Magnitud
subplot(221)
plot(w,20*log10(abs(H)));
grid on;
title ([titutlo, 'Respuesta en magnitud']);
xlabel('Frecuencia (Hz)');
ylabel('H(f) db')
% xlim([0 0.4])
% Respuesta en fase
subplot(222)
plot(w,angle(H));
grid on;
title ([titutlo,'Respuesta en fase']);
xlabel('Frecuencia (Hz)')
ylabel('ángulo de H rad')
% xlim([0 0.4])
%Respuesta al impulso
subplot(223)
[y_eje,tt]= impz(b,a,60);
stem(tt,y_eje);
title ([titutlo,'Respuesta al impulso']);
%Ploteo de los polos y ceros
z= roots(b); % Ceros
p = roots(a); % Polos
subplot(224)
zplane(z,p)
title('Polos y ceros')
legend('Ceros','Polos')
%% Reproducción de audio de entrada y salida
pause(2)
disp('Audio de entrada')
% Se multiplica por 0.1 para atenuar la salida del tono por la bocina
wavplay(0.1*y,fs)
disp('Audio de salida (señal filtrada)')
% Se multiplica por 0.1 para atenuar la salida del tono por la bocina
wavplay(0.1*y_Low,fs)
%% Procesamiento de una señal de audio usando MATLAB
%% Selección del tipo de filtrado
% 1 -> Pasa bajo
% 2 -> Pasa alto
% 3 -> Pasa banda
tipo=3;
%% Leer la señal de audio
[y,fs]=wavread('Yabu_mono');
% y -> muestras de la señal
% fs-> frecuencia de muestreo
% Graficar la señal original
subplot(411)
T=length(y)/fs;
t=linspace(0,T,T*fs);
plot(t,y)
title('SEÑAL ORIGINAL')% Título
xlabel('Tiempo (s)') % Etiqueta del eje X
ylabel('Amplitud (V)') % Etiqueta del eje Y
% xlim([0 20/]) % Límite de la señal
%% Reproducir la señal de audio
% wavplay(0.1*y,fs)
%% FFT de la señal
subplot(412)
% Llamado a la función que calcula la FFT
fft_signal(y,fs);title('ESPECTRO DE LA SEÑAL ORIGINAL')
xlim([0 4e3])
%% Filtrado de la señal
switch tipo
case 1
% Cálculo de los coeficientes del filtro (filtro pasa bajas)
% Este filtrado deja solo la señal por debajo de 500 Hz
titulo='FILTRO PASA BAJAS';
% Frecuencia normalizada
fNorm = 1000 / (fs/2);
[b,a] = butter(10, fNorm, 'low');
case 2
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Cálculo de los coeficientes del filtro (filtro pasa altas)
% Este filtrado deja solo la señal por encima de 500 Hz
titulo='FILTRO PASA ALTAS';
% Frecuencia normalizada
fNorm = 1000 / (fs/2);
[b,a] = butter(10, fNorm, 'high');
otherwise
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Cálculo de los coeficientes del filtro (filtro pasa banda)
% Este filtrado deja solo la señal de 2KHz a 3KHz
% Frecuencias normalizadas
titulo='FILTRO PASA BANDA';
Wp = [2e3 3e3]/(fs/2); Ws = [1.5e3 3.5e3]/(fs/2);
Rp = 3; Rs = 40; % Rizado de la banda de paso y de parada (s)
[n,Wn] = buttord(Wp,Ws,Rp,Rs);% Orden del filtro y frecuencia de corte óptima
[b,a] = butter(n,Wn); % Coeficientes del filtro
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end
% Filtrado de la señal
y_Low = filtfilt(b, a, y);
% Graficación de la señal en el tiempo
subplot(413)
plot(t,y_Low)
title('SEÑAL FILTRADA')
xlabel('Tiempo (s)')
ylabel('Amplitud (V)')
% xlim([0 20/f0])
% Graficación de la señal en frecuencia
subplot(414)
% Llamado a la función que calcula la FFT
fft_signal(y_Low,fs);title('ESPECTRO DE LA SEÑAL FILTRADA')
xlim([0 4e3])
%% Gráficas del filtro
% Respuesta en frecuencia del filtro
[H,w]=freqz(b,a,512,fs);
figure(2)
%Trazado de la respuesta en Magnitud
subplot(221)
plot(w,20*log10(abs(H)));
grid on;
title ([titulo, ' Respuesta en magnitud']);
xlabel('Frecuencia (Hz)');
ylabel('H(f) db')
xlim([0 5e3])
% Respuesta en fase
subplot(222)
plot(w,angle(H));
grid on;
title ([titulo,' Respuesta en fase']);
xlabel('Frecuencia (Hz)')
ylabel('ángulo de H rad')
xlim([0 5e3])
%Respuesta al impulso
subplot(223)
[y_eje,tt]= impz(b,a,60);
stem(tt,y_eje);
title ([titulo,' Respuesta al impulso']);
xlabel('n')
ylabel('h[n]')
%Ploteo de los polos y ceros
z= roots(b); % Ceros
p = roots(a); % Polos
subplot(224)
zplane(z,p)
title('Polos y ceros')
legend('Ceros','Polos')
%% Reproducción de audio de entrada y salida
pause(2)
disp('Audio de entrada')
% Se multiplica por 0.2 para atenuar la salida del tono por la bocina
% wavplay(0.2*y,fs)
disp('Audio de salida (señal filtrada)')
% Se multiplica por 0.2 para atenuar la salida del tono por la bocina
% wavplay(0.2*y_Low,fs)
Referencias
