Como plotar manualmente a resposta de frequência de um filtro Butterworth passa-banda no MATLAB sem a função freqz?

Eu tenho o código abaixo que aplica um filtro passa-banda em um sinal. Eu sou um idiota no DSP e quero entender o que está acontecendo nos bastidores antes de prosseguir.

Para fazer isso, quero saber como plotar a resposta de frequência do filtro sem usar freqz.

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

Dadas as saídas, [b, a]como eu faria isso? Parece que seria uma tarefa simples, mas estou tendo dificuldades para encontrar o que preciso na documentação ou online.

Eu também gostaria de entender como fazer isso o mais rápido possível, por exemplo, usando um fftou outro algoritmo rápido.

Sabemos que, em geral, a função de transferência de um filtro é dada por:

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

Agora substitua $z=e^{j\omega}$ para avaliar a função de transferência no círculo unitário:

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

Assim, isso se torna apenas um problema de avaliação polinomial em um dado $\omega$ . Aqui estão os passos:

• Crie um vetor de frequências angulares $\omega = [0, \ldots,\pi]$ para a primeira metade do espectro (não é necessário ir até $2\pi$ ) e salve-o w.
• Pré-cálculo do expoente $e^{-j\omega}$ em todos eles e armazene-o em variável ze.
• Use a polyvalfunção para calcular os valores do numerador e denominador chamando polyval(b, ze):, divida-os e armazene-os H. Porque estamos interessados ​​em amplitude, então pegue o valor absoluto do resultado.
• Converta para a escala dB usando: $H_{dB}=20\log_{10} H$ - neste caso $1$ é o valor de referência.

Colocando tudo isso no 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)

Saída original de freqz:

E a saída do meu script:

E comparação rápida em escala linear - parece ótimo!

[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'})

Agora você pode reescrevê-lo em alguma função e adicionar poucas condições para torná-lo mais útil.

---

Outra maneira (proposta anteriormente é mais confiável) seria usar a propriedade fundamental, que a resposta de frequência de um filtro é uma Transformada de Fourier de sua resposta de impulso:

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

Portanto, você deve alimentar o sinal $\delta(t)$ , calcular a resposta do seu filtro e obter a FFT dele:

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));

Em comparação, isso produzirá o seguinte:

este é apenas um adendo à resposta de jojek, que é mais geral e perfeitamente boa quando se usa matemática de dupla precisão. quando há menos precisão, ocorre um "problema de cosseno" que surge quando a frequência na resposta de frequência é muito baixa (muito menor que Nyquist) e também quando as frequências ressonantes do filtro são muito baixas.

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

considere esta identidade trigonométrica:

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

$\sin^2 \left( \frac{\omega}{2} \right)$$\omega \to 0$ . tão pequeno que o termo e as informações de frequência desse termo se perdem quando o termo (ou é negativo) é adicionado a 1. esse é o caso mesmo do ponto flutuante, mas é menos um problema com flutuadores de precisão dupla. mas alguns de nós que colocamos essa função de resposta em frequência em marcha podem não ter o recurso de ponto flutuante de precisão dupla ou qualquer ponto flutuante.

$\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 tem resposta de frequência complexa

$$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 tem magnitude ao quadrado:

$$\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}$$

so, one can see that the magnitude frequency response $|H(e^{j\omega})|$ is an even symmetry function and depends only on the cosines $\cos(\omega)$ and $\cos(2\omega)$. for very low $\omega$, the values of those cosines are so close to $1$ that, with single-precision fixed or floating point, there are few bits remaining that differentiate those values from $1$. that is the "cosine problem".

using the trig identity above, you get for magnitude squared:

$$\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.