Wie zeichne ich den Frequenzgang eines Butterworth-Bandpassfilters in MATLAB manuell ohne Frequenzfunktion?

Ich habe Code wie unten, der ein Bandpassfilter auf ein Signal anwendet. Ich bin ein ziemlicher Neuling bei DSP und möchte verstehen, was sich hinter den Kulissen abspielt, bevor ich fortfahre.

Um dies zu tun, möchte ich wissen, wie man den Frequenzgang des Filters grafisch darstellt, ohne ihn zu verwenden freqz.

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

[b, a]Wie würde ich das tun, wenn ich die Ausgaben annehmen würde? Dies scheint eine einfache Aufgabe zu sein, aber es fällt mir schwer, das zu finden, was ich in der Dokumentation oder online benötige.

Ich möchte auch verstehen, wie das so schnell wie möglich geht, z. B. mit einem fftoder anderen schnellen Algorithmus.

Wir wissen, dass die Übertragungsfunktion eines Filters im Allgemeinen gegeben ist durch:

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

Setzen Sie nun $z=e^{j\omega}$ , um die Übertragungsfunktion auf dem Einheitskreis zu bewerten:

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

Somit wird dies nur ein Problem der Polynomauswertung bei einem gegebenen $\omega$ . Hier sind die Schritte:

• Erstellen Sie einen Vektor der Winkelfrequenzen $\omega = [0, \ldots,\pi]$ für die erste Hälfte des Spektrums (Sie müssen nicht auf $2\pi$ ) und speichern Sie ihn in w.
• Berechnen Sie den Exponenten $e^{-j\omega}$ und speichern Sie ihn in einer Variablen ze.
• Verwenden Sie die polyvalFunktion, um die Werte von Zähler und Nenner durch Aufrufen von: zu berechnen polyval(b, ze), zu teilen und in zu speichern H. Da wir an der Amplitude interessiert sind, nehmen wir den absoluten Wert des Ergebnisses.
• In dB-Skala umwandeln mit: $H_{dB}=20\log_{10} H$ - in diesem Fall ist $1$ der Referenzwert.

All das in Code einfügen:

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

Ursprüngliche Ausgabe von freqz:

Und die Ausgabe meines Skripts:

Und schneller Vergleich im linearen Maßstab - sieht gut aus!

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

Jetzt können Sie es in eine Funktion umschreiben und einige Bedingungen hinzufügen, um es nützlicher zu machen.

---

Ein anderer Weg (der zuvor vorgeschlagen wurde, ist zuverlässiger) wäre die Verwendung der fundamentalen Eigenschaft, dass der Frequenzgang eines Filters eine Fourier-Transformation seines Impulsgangs ist:

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

$\delta(t)$

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

Im Vergleich ergibt sich folgendes:

Dies ist nur ein Nachtrag zu Jojeks Antwort, die allgemeiner gehalten ist und bei Verwendung von Mathematik mit doppelter Genauigkeit perfekt funktioniert. Wenn die Genauigkeit geringer ist, tritt ein "Kosinusproblem" auf, das auftritt, wenn entweder die Frequenz im Frequenzgang sehr niedrig ist (viel niedriger als in Nyquist), oder wenn die Resonanzfrequenzen des Filters sehr niedrig sind.

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

Betrachten Sie diese Trigger-Identität:

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

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

Also, was ich getan habe, ist, die oben angegebene Triggeridentität zu verwenden und alle Cosinus-Terme zu entfernen und sie im Wesentlichen durch Terme zu ersetzen, die wie aussehen$\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}}$$

das hat einen komplexen Frequenzgang

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

welches die quadratische Größe hat:

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

Mit der obigen Triggeridentität erhalten Sie das Quadrat der Größe:

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