Funcionamiento filtro IIR

Estoy diseñando un filtro digital y tengo dudas sobre su funcionamiento. Se trata de una pregunta muy importante, y muy larga así que te la formulare en varias partes. Así que no tenga ninguna duda que cualquier respuesta suya sera valorada, por tantas partes como tenga la pregunta, creo que es lo más justo.
PARTE 1:
Encontré este documento, acerca de los filtros IIR:
http://www.musicdsp.org/files/EQ-Coefficients.pdf
El filtro digital que se pretende diseñar es tal como se explica en el apartado 3 del PDF. Es un filtro que amplifica o atenúa en la frecuencia central, pero que actúa como filtro pasa todo en el resto de frecuencias. Parte de la ecuación:
H = (s^2 + 2*K*apha*wo*s + wo^2)/(s^2 + 2*alpha*wo*s + wo*2)
Don de QUE es la amplificación, wo la frecuencia central y alpha es un parámetro que depende del ancho de banda y de QUE, así se indica en la página 6 del PDF.
Finalmente se llega a la conclusión que la función de transferencia es la que se muestra en la ecuación (11), página 6. Por mi parte yo hice esta function en Matlab para hallar los coeficientes de la función de transferencia en función de 's':
Los valores de entrada son G (amplficacion en dB), bw(bandwith) y wo(frecuencia central en rad/s).
function[num, den, K, alpha] = coef(G, bw, wo) 
K = 10^(G/20);     -- K es la amplificacion.
alpha = (1/sqrt(K))*sinh((bw*log(2))/2); -- este es el parametro alpha.
a0 = 1;                       -- a0 es el coeficiente que multiplica s^2 en el numerador.
a1 = 2*K*alpha*wo; -- a1 es el coeficiente que multiplica s en el numerador.
a2 = wo^2;               -- a2 es el ultimo coeficiente del numerador.
b0 = 1;                      -- b0 es el coeficiente que multiplica s^2 en el denominador.
b1 = 2*alpha*wo;    -- b1 es el coeficiente que multiplica s en el denominador.
b2 = wo^2;              -- b2 es el ultimo coeficiente del denominador.
num = [a0 a1 a2];    -- Agrupo los tres coeficientes del numerador en un vector.
den = [b0 b1 b2];    -- Agrupo los tres coeficientes del denominador en un vector.
g = tf(num, den);    -- Esto seria para visualizar la funcion de transferencia.
end
Mediante una segunda función creada en Matlab paso los coeficientes calculados en 's' en el dominio de 'z' mediante la transformada bilineal.
En esta function los valores de entrada son los vectores num, den calculados en la anterior function, y fs (frecuencia de muestreo). Mediante la instrucción "bilinear" hago la transformada.
function [numd, dend] = coefz(num, den, fs)
[numd, dend] = bilinear(num, den, fs);       -- paso los coeficientes calculados al dominio z
h = tf(numd, dend, 1/fs, 'variable', 'z^-1'); -- funcion de transferencia.
end
donde numd = [a0 a1 a2] y dend = [b0 b1 b2].
La funcion de transferencia resultante en z seria esta:
H (z) = (a0 + a1*z^-1 + a2*z^-2)/(b0 + b1*z^-1 + b2*z^-2)