¿Cómo graficar manualmente la respuesta de frecuencia de un filtro Butterworth de paso de banda en MATLAB sin la función freqz?

Tengo el siguiente código que aplica un filtro de paso de banda a una señal. Soy un novato en DSP y quiero entender lo que está sucediendo detrás de escena antes de continuar.

Para hacer esto, quiero saber cómo trazar la respuesta de frecuencia del filtro sin usar freqz.

[b, a] = butter(order, [flo fhi]);
filtered_signal = filter(b, a, unfiltered_signal)

Dados los resultados, [b, a]¿cómo haría esto? Parece que sería una tarea simple, pero me cuesta encontrar lo que necesito en la documentación o en línea.

También me gustaría entender cómo hacer esto lo más rápido posible, por ejemplo, usando un fftu otro algoritmo rápido.

Sabemos que, en general, la función de transferencia de un filtro viene dada por:

$$H(z)=\dfrac{\sum_{k=0}^{M}b_kz^{-k}}{\sum_{k=0}^{N}a_kz^{-k}}$$

Ahora sustituya para evaluar la función de transferencia en el círculo unitario:$z=e^{j\omega}$

$$H(e^{j\omega})=\dfrac{\sum_{k=0}^{M}b_ke^{-j\omega k}}{\sum_{k=0}^{N}a_ke^{-j\omega k}}$$

Por lo tanto, esto se convierte solo en un problema de evaluación polinómica en un determinado . Aquí están los pasos:$\omega$

• Cree un vector de frecuencias angulares para la primera mitad del espectro (no es necesario subir a ) y guárdelo .$\omega = [0, \ldots,\pi]$$2\pi$w
• Pre-calcular el exponente en todos ellos y almacenarlo en variable .$e^{-j\omega}$ze
• Use la polyvalfunción para calcular los valores de numerador y denominador llamando: polyval(b, ze)dividirlos y almacenarlos H. Como estamos interesados ​​en la amplitud, tomemos el valor absoluto del resultado.
• Convierta a escala dB utilizando: - en este caso es el valor de referencia.$H_{dB}=20\log_{10} H$$1$

Poniendo todo eso en el código:

%% Filter definition
a = [1 -0.5 -0.25]; % Some filter with lot's of static gain
b = [1 3 2];

%% My freqz calculation
N = 1024; % Number of points to evaluate at
upp = pi; % Evaluate only up to fs/2
% Create the vector of angular frequencies at one more point.
% After that remove the last element (Nyquist frequency)
w = linspace(0, pi, N+1); 
w(end) = [];
ze = exp(-1j*w); % Pre-compute exponent
H = polyval(b, ze)./polyval(a, ze); % Evaluate transfer function and take the amplitude
Ha = abs(H);
Hdb  = 20*log10(Ha); % Convert to dB scale
wn   = w/pi;
% Plot and set axis limits
xlim = ([0 1]);
plot(wn, Hdb)
grid on

%% MATLAB freqz
figure
freqz(b,a)

Salida original de freqz:

Y la salida de mi script:

Y una comparación rápida en escala lineal: ¡se ve genial!

[h_f, w_f] = freqz(b,a);
figure
xlim = ([0 1]);
plot(w, Ha) % mine
grid on
hold on
plot(w_f, abs(h_f), '--r') % MATLAB
legend({'my freqz','MATLAB freqz'})

Ahora puede reescribirlo en alguna función y agregar algunas condiciones para que sea más útil.

---

Otra forma (propuesta anteriormente es más confiable) sería usar la propiedad fundamental, que la respuesta de frecuencia de un filtro es una Transformada de Fourier de su respuesta de impulso:

$$H(\omega)=\mathcal{F}\{h(t)\}$$

Por lo tanto, debe alimentar la señal $\delta(t)$ su sistema , calcular la respuesta de su filtro y tomar la FFT de la misma:

d = [zeros(1,length(w_f)) 1 zeros(1,length(w_f)-1)];
h = filter(b, a, d);
HH = abs(fft(h));
HH = HH(1:length(w_f));

En comparación, esto producirá lo siguiente:

esto es solo un apéndice a la respuesta de jojek que es más general y perfectamente bueno cuando se usa matemática de doble precisión. cuando hay menos precisión, hay un "problema de coseno" que surge cuando la frecuencia en la respuesta de frecuencia es muy baja (mucho más baja que Nyquist) y también cuando las frecuencias resonantes del filtro son muy bajas.

$|H(e^{j\omega})|^2$$|H(e^{-j\omega})| = |H(e^{j\omega})|$

considere esta identidad trigonométrica:

$$\cos(\omega) \ = \ 1 - 2 \sin^2 \left( \frac{\omega}{2} \right)$$

$\sin^2 \left( \frac{\omega}{2} \right)$$\omega \to 0$

$\sin^2 \left( \frac{\omega}{2} \right)$

$$H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{a_0 + a_1 z^{-1} + a_2 z^{-2}}$$

que tiene una respuesta de frecuencia compleja

$$H(e^{j\omega}) = \frac{b_0 + b_1 e^{-j\omega} + b_2 e^{-j2\omega}}{a_0 + a_1 e^{-j\omega} + a_2 e^{-j2\omega}}$$

que tiene magnitud al cuadrado:

$$\begin{align} |H(e^{j\omega})|^2 &= \frac{|b_0 + b_1 e^{-j\omega} + b_2 e^{-j2\omega}|^2}{|a_0 + a_1 e^{-j\omega} + a_2 e^{-j2\omega}|^2} \\ \\ &= \frac{\big(b_0 + b_1\cos(\omega) + b_2\cos(2\omega)\big)^2 + \big(b_1\sin(\omega) + b_2\sin(2\omega)\big)^2}{\big(a_0 + a_1\cos(\omega) + a_2\cos(2\omega)\big)^2 + \big(a_1\sin(\omega) + a_2\sin(2\omega)\big)^2} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)\cos(\omega) + 2b_0b_2\cos(2\omega)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)\cos(\omega) + 2a_0a_2\cos(2\omega)} \\ \end{align}$$

$|H(e^{j\omega})|$$\cos(\omega)$$\cos(2\omega)$$\omega$$1$$1$

usando la identidad trigonométrica anterior, obtienes una magnitud al cuadrado:

$$\begin{align} |H(e^{j\omega})|^2 &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)\cos(\omega) + 2b_0b_2\cos(2\omega)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)\cos(\omega) + 2a_0a_2\cos(2\omega)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2b_0b_2\left(1 - 2 \sin^2(\omega)\right)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2a_0a_2\left(1 - 2 \sin^2(\omega)\right)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2b_0b_2\left(2\cos^2(\omega) - 1\right)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2a_0a_2\left(2\cos^2(\omega) - 1\right)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2b_0b_2\left(2\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right)^2 - 1\right)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right) + 2a_0a_2\left(2\left(1 - 2 \sin^2 \left( \tfrac{\omega}{2}\right)\right)^2 - 1\right)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)(1 - 2\phi) + 2b_0b_2\left(2(1 - 2\phi )^2 - 1\right)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)(1 - 2\phi) + 2a_0a_2\left(2(1 - 2\phi)^2 - 1\right)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1(b_0+b_2)(1 - 2\phi) + 2b_0b_2(1 - 8\phi + 8\phi^2)}{a_0^2+a_1^2+a_2^2 + 2a_1(a_0+a_2)(1 - 2\phi) + 2a_0a_2(1 - 8\phi + 8\phi^2)} \\ \\ &= \frac{b_0^2+b_1^2+b_2^2 + 2b_1b_0+2b_1b_2 - 4(b_1b_0+b_1b_2)\phi + 2b_0b_2 - 16b_0b_2\phi + 16b_0b_2\phi^2}{a_0^2+a_1^2+a_2^2 + 2a_1a_0+2a_1a_2 - 4(a_1a_0+a_1a_2)\phi + 2a_0a_2 - 16a_0a_2\phi + 16a_0a_2\phi^2} \\ \\ &= \frac{\big(b_0^2+b_1^2+b_2^2 + 2b_1b_0+2b_1b_2+2b_0b_2\big) - 4(b_1b_0+b_1b_2-4b_0b_2)\phi + 16b_0b_2\phi^2}{\big(a_0^2+a_1^2+a_2^2 + 2a_1a_0+2a_1a_2+2a_0a_2\big) - 4(a_1a_0+a_1a_2-4a_0a_2)\phi + 16a_0a_2\phi^2} \\ \\ &= \frac{\tfrac14\big(b_0^2+b_1^2+b_2^2 + 2b_1b_0+2b_1b_2+2b_0b_2\big) - (b_1b_0+b_1b_2-4b_0b_2)\phi + 4b_0b_2\phi^2}{\tfrac14\big(a_0^2+a_1^2+a_2^2 + 2a_1a_0+2a_1a_2+2a_0a_2\big) - (a_1a_0+a_1a_2-4a_0a_2)\phi + 4a_0a_2\phi^2} \\ \\ &= \frac{\left(\frac{b_0+b_1+b_2}{2}\right)^2 - \phi \big(4b_0b_2(1-\phi) + b_1(b_0+b_2)\big)}{\left(\frac{a_0+a_1+a_2}{2}\right)^2 - \phi \big(4a_0a_2(1-\phi) + a_1(a_0+a_2)\big)} \\ \end{align}$$

where $\phi \triangleq \sin^2 \left( \frac{\omega}{2} \right)$

if your gear is intending to plot this as dB, it comes out as

$$20 \log_{10}|H(e^{j\omega})| \ = \ 10 \log_{10}\left( \left(\tfrac{b_0+b_1+b_2}{2}\right)^2 - \phi \big(4b_0b_2(1-\phi) + b_1(b_0+b_2)\big) \right) \\ - 10 \log_{10}\left(\left(\tfrac{a_0+a_1+a_2}{2}\right)^2 - \phi \big(4a_0a_2(1-\phi) + a_1(a_0+a_2)\big) \right)$$

so your division turns into subtraction, but you have to be able to compute logarithms to some base or another. numerically, you will have much less trouble with this for low frequencies than doing it the apparent way.