Comment tracer manuellement la réponse en fréquence d'un filtre Butterworth passe-bande dans MATLAB sans fonction freqz?

J'ai un code comme ci-dessous qui applique un filtre passe-bande sur un signal. Je suis tout à fait un noob chez DSP et je veux comprendre ce qui se passe dans les coulisses avant de continuer.

Pour ce faire, je veux savoir comment tracer la réponse en fréquence du filtre sans utiliser freqz.

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

Compte tenu des sorties, [b, a]comment pourrais-je faire cela? Cela semble être une tâche simple, mais j'ai du mal à trouver ce dont j'ai besoin dans la documentation ou en ligne.

J'aimerais également comprendre comment procéder le plus rapidement possible, par exemple en utilisant un fftou un autre algorithme rapide.

On sait qu'en général la fonction de transfert d'un filtre est donnée par:

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

Remplacez maintenant $z=e^{j\omega}$ pour évaluer la fonction de transfert sur le cercle unitaire:

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

Cela devient donc seulement un problème d'évaluation polynôme à une donnée $\omega$ . Voici les étapes:

• Créez un vecteur de fréquences angulaires $\omega = [0, \ldots,\pi]$ pour la première moitié du spectre (pas besoin de monter jusqu'à $2\pi$ ) et enregistrez-le w.
• Pré-calculer l'exposant $e^{-j\omega}$ sur chacun d'eux et le stocker dans une variable ze.
• Utilisez la polyvalfonction pour calculer les valeurs du numérateur et du dénominateur en appelant polyval(b, ze):, divisez-les et enregistrez-les H. Parce que nous nous intéressons à l'amplitude, prenons alors la valeur absolue du résultat.
• Convertir en échelle dB en utilisant: $H_{dB}=20\log_{10} H$ - dans ce cas, $1$ est la valeur de référence.

Mettre tout cela dans le code:

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

Sortie originale de freqz:

Et la sortie de mon script:

Et une comparaison rapide à l'échelle linéaire - a fière allure!

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

Vous pouvez maintenant le réécrire dans une fonction et ajouter quelques conditions pour le rendre plus utile.

---

Une autre façon (précédemment proposée est plus fiable) serait d'utiliser la propriété fondamentale, que la réponse en fréquence d'un filtre est une transformée de Fourier de sa réponse impulsionnelle:

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

Par conséquent, vous devez alimenter le signal $\delta(t)$ votre système , calculer la réponse de votre filtre et en prendre la FFT:

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 comparaison, cela produira ce qui suit:

ce n'est qu'un addendum à la réponse de jojek qui est plus générale et parfaitement bonne lorsque les mathématiques en double précision sont utilisées. quand il y a moins de précision, il y a un "problème de cosinus" qui surgit lorsque la fréquence dans la réponse en fréquence est très basse (beaucoup plus basse que Nyquist) et aussi quand les fréquences de résonance du filtre sont très basses.

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

considérez cette identité trig:

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

qui a une réponse en fréquence complexe

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

qui a une magnitude au carré:

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

en utilisant l'identité trig ci-dessus, vous obtenez pour la grandeur au carré:

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