// Prueba de las funciones iir() y filter()


n = 3;                    // Orden del filtro
ftype = "bp";             // Pasabanda
fdesign = "butt";         // Butterworth
Fs = 44100;               // Tasa de muestreo 
fo = 1000;                // Frecuencia central
fi = fo/sqrt(2);          // Frecuencia inferior
fs = fo*sqrt(2);          // Frecuencia superior 
frq = [fi/Fs fs/Fs];      // Frecuencias normalizadas
delta = [0 0]             // Argumento irrelevante para 
                          // Butterworth pero necesario
                          // sintácticamente

// Diseñar el filtro digital 
[p, z, g] = iir(n, ftype, fdesign, frq, delta)

// Función de transferencia digital
h = iir(n, ftype, fdesign, frq, delta)

// Vector de frecuencias
f = 20*10.^[0:0.01:3];

// Vector de frecuencias angulares normalizado
omega = 2*%pi*f/Fs;

// Respuesta en frecuencia
RF = horner(h, exp(%i*omega));

// Prototipo analógico
// Función de transferencia del pasabajos
hs = analpf(n, fdesign, [0 0], 1)
// Frecuencia angular central
wo = 2*%pi*fo
// Ancho de banda relativo
B = (fs - fi)/fo
// Transformación PB --> PBn 
H = horner(hs, (wo/%s + %s/wo)/B)
// Respueste en frecuencia analógica
RFs = horner(H, %i*2*%pi*f);

// Gráfica de la respuesta en frecuencia
scf(1);
clf(1);
semilogx(f, 20*log10(abs(RFs)),"r", f, 20*log10(abs(RF)))
title("Respuesta en frecuencia de filtro pasabanda de Butterworth")
xlabel("f  (Hz)")
ylabel("|H|")
legend("Analógico","Digital")
gca().grid = [1 1];

// Carga de una señal de audio
// NOTA: sustituir por una señal del usuario)
x = wavread("r:\Mis documentos\Acústica\Cursos\Scilab\prueba_de_audio.wav");
x = x(1,:);

// Obtencion del numerador y denominador de
// la función de transferencia digital
NUM = h.num
DEN = h.den

// Filtrado digital
y = filter(NUM, DEN, x);

exec fftmean.sci;
// Obtención de los espectros FFT promedio
// de la señal original y filtrada
[X2, X, Xmin, Xmax, fo] = fftmean(x, Fs, N, "hann", N/2);
[Y2, Y, Ymin, Ymax, fo] = fftmean(y, Fs, N, "hann", N/2);


playsnd(x, Fs, 1)
playsnd(5*y, Fs)

scf(2);
clf(2);
//plot("ln", fo, Y2, fo, X2)
plot("ln", fo, 20*log10(X2), fo, 20*log10(Y2), "r")
title("Espectros de la señal original y filtrada")
legend("original", "filtrada")
xlabel("f  (Hz)")
ylabel("Y  (dB)")



