Pytanie matematyczne wynikające z zastosowania transformacji dwuliniowej

Jest to związane z książką kucharską i próbowałem rozwiązać ją może dwie dekady temu, poddałem się i przypomniałem sobie o nierozwiązanym problemie. Ale to jest cholernie proste, ale wciąż jestem zalany w błocie.

Jest to prosty filtr pasmowy (BPF) o częstotliwości rezonansowej $\Omega_0$ i rezonans $Q$:

$$H(s) = \frac{\frac{1}{Q}\frac{s}{\Omega_0}}{\left(\frac{s}{\Omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\Omega_0} + 1}$$

Przy częstotliwości rezonansowej

$$|H(j\Omega)| \le H(j\Omega_0) = 1$$

a górne i dolne pasy są zdefiniowane w ten sposób

$$|H(j\Omega_U)|^2 = \left|H\left(j\Omega_0 2^{BW/2} \right)\right|^2 = \tfrac12$$

$$|H(j\Omega_L)|^2 = \left|H\left(j\Omega_0 2^{-BW/2} \right)\right|^2 = \tfrac12$$

Nazywamy to „bandażami połowy mocy” . Ponieważ jesteśmy audio, określamy szerokość pasma w oktawach, aw świecie analogowym to pasmo w oktawach,$BW$, odnosi się do $Q$ tak jak:

$$\frac{1}{Q} = \frac{2^{BW} - 1}{\sqrt{2^{BW}}} = 2 \sinh\left( \tfrac{\ln(2)}{2} BW \right)$$

Używamy transformacji dwuliniowej (z wstępnie wypaczoną częstotliwością rezonansową), która odwzorowuje:

$$\begin{align} \frac{s}{\Omega_0} & \leftarrow \frac{1}{\tan(\omega_0/2)} \, \frac{1-z^{-1}}{1+z^{-1}} \\ \\ \frac{j \Omega}{\Omega_0} & \leftarrow \frac{j\tan(\omega/2)}{\tan(\omega_0/2)} \\ \end{align}$$

najmu $z=e^{j \omega}$ i $s=j\Omega$.

Rezonansowa częstotliwość kątowa filtra analogowego wynosi$\Omega_0$, oraz z kompensacją dopasowania częstotliwości do częstotliwości rezonansowej w zrealizowanym filtrze cyfrowym, kiedy $\omega=\omega_0$ (częstotliwość rezonansowa zdefiniowana przez użytkownika) $\Omega=\Omega_0$.

Więc jeśli analogowa częstotliwość kątowa jest

$$\frac{\Omega}{\Omega_0} = \frac{\tan(\omega/2)}{\tan(\omega_0/2)}$$

następnie jest mapowany na cyfrową częstotliwość kątową jako

$$\omega = 2 \arctan\left( \frac{\Omega}{\Omega_0} \, \tan(\omega_0/2) \right)$$

Teraz górne i dolne pasma w świecie analogowym są

$$\Omega_U = \Omega_0 2^{BW/2}$$ $$\Omega_L = \Omega_0 2^{-BW/2}$$

i w dziedzinie częstotliwości cyfrowych są

$$\begin{align} \omega_U &= 2 \arctan\left( \frac{\Omega_U}{\Omega_0} \, \tan(\omega_0/2) \right) \\ &= 2 \arctan\left(2^{BW/2} \, \tan(\omega_0/2) \right) \\ \end{align}$$

$$\begin{align} \omega_L &= 2 \arctan\left(\frac{\Omega_L}{\Omega_0} \, \tan(\omega_0/2) \right) \\ &= 2 \arctan\left(2^{-BW/2} \, \tan(\omega_0/2) \right) \\ \end{align}$$

Zatem rzeczywista różnica w częstotliwości rejestrowania pasm (która jest faktyczną przepustowością w filtrze cyfrowym) wynosi:

$$\begin{align} bw & = \log_2(\omega_U) - \log_2(\omega_L) \\ & = \log_2\Bigg(2\arctan\left( 2^{ BW/2} \, \tan(\omega_0/2) \right) \Bigg) - \log_2\Bigg(2\arctan\left( 2^{-BW/2} \, \tan(\omega_0/2) \right) \Bigg) \ \end{align}$$

lub

$$\ln(2) \, bw = \ln\left(\arctan\left( e^{\ln(2)BW/2} \, \tan(\omega_0/2) \right) \right) - \ln\left(\arctan\left( e^{-\ln(2)BW/2} \, \tan(\omega_0/2) \right) \right)$$

Ma to funkcjonalną formę

$$f(x) = \ln\left(\arctan\left(\alpha \, e^{x} \right) \right) - \ln\left(\arctan\left( \alpha \, e^{-x}\right) \right)$$

gdzie $f(x) \triangleq \ln(2) \, bw$, $x \triangleq \frac{\ln(2)}{2}BW$ i $\alpha \triangleq \tan(\omega_0/2)$

Chcę tylko odwrócić $f(x)$(ale wiem, że nie mogę tego zrobić z ładną, zamkniętą postacią). Dokonałem już przybliżenia pierwszego rzędu i chcę podnieść go do przybliżenia trzeciego rzędu. I stało się to rodzajem kopulującego samica psa, mimo że powinno być proste.

Teraz ma to coś wspólnego z formułą inwersji Lagrange'a i chcę przejść na jeszcze jeden termin.

Wiemy to z góry $f(x)$ jest funkcją nieparzystej symetrii:

$$f(-x) = -f(x)$$

To znaczy $f(0)=0$ a wszystkie warunki parzystej kolejności serii Maclaurin będą wynosić zero:

$$y = f(x) = a_1 x + a_3 x^3 + ...$$

Funkcja odwrotna jest również nieparzystą symetrią, przechodzi przez zero i może być wyrażona jako szereg Maclaurin

$$x = g(y) = b_1 y + b_3 y^3 + ...$$

a jeśli wiemy co $a_1$ i $a_3$ są z $f(x)$, to mamy dobry pomysł co $b_1$ i $b_3$ musi być:

$$b_1=\frac{1}{a_1} \quad \quad \quad \quad b_3 = -\frac{a_3}{a_1^4}$$

Teraz jestem w stanie obliczyć pochodną $f(x)$ i oceniam to na zero i dostaję

$$\\ a_1 = \frac{2 \alpha}{(1+\alpha^2) \arctan(\alpha)} = \frac{\sin(\omega_0)}{\omega_0/2}$$ $$\\ b_1 = \frac{(1+\alpha^2) \arctan(\alpha)}{2 \alpha} = \frac{\omega_0/2}{\sin(\omega_0)} \\ $$

Ale mam dziwkę czasu $a_3$ i dlatego $b_3$. Czy ktoś może to zrobić? Zadowoliłbym się nawet solidnym wyrazem trzeciej pochodnej$f(x)$ oceniono na $x=0$.

Uzupełnienie mojej części do tego pytania: Oto nieco krótka odpowiedź oparta na ręcznym rozszerzeniu funkcji nieparzystej $f(x)$

$$\begin{align*} f(x)&=\ln\left(\arctan\left(\alpha e^x\right)\right)-\ln\left(\arctan\left(\alpha e^{-x}\right)\right)\\ &=f_1x+f_3x^3+O\left(x^5\right)\tag{1} \end{align*}$$ w serię do trzeciego rzędu. Więcej szczegółów można znaleźć na stronie mathSE .

Na początku skupiamy się na lewostronnym okresie

$$\ln\left(\arctan\left(\alpha e^x\right)\right)$$ z $f(x)$ i zacznij od

Rozszerzenie serii $\arctan$:

Otrzymujemy

$$\begin{align*} \arctan(\alpha e^x)&=\sum_{n=0}^\infty \frac{(-1)^n}{2n+1}\alpha^{2n+1} e^{(2n+1)x}\\ &=\cdots\\ &=\sum_{j=0}^\infty\frac{1}{j!}\sum_{n=0}^\infty(-1)^n(2n+1)^{j-1}\alpha^{2n+1}x^j\tag{2} \end{align*}$$

Teraz wywodzimy z (2) współczynników do $x^3$. Korzystanie ze współczynnika operatora$[x^k]$ oznaczać współczynnik $x^k$ w serii, którą otrzymujemy

$$\begin{align*} [x^0]\arctan\left(\alpha e^x\right)&=\sum_{n=0}^\infty \frac{(-1)^n}{2n+1}\alpha^{2n+1}=\arctan \alpha\\ [x^1]\arctan\left(\alpha e^x\right)&=\sum_{n=0}^\infty (-1)^n\alpha^{2n+1}=\frac{\alpha}{1+\alpha ^2}\\ [x^2]\arctan\left(\alpha e^x\right)&=\frac{1}{2}\sum_{n=0}^\infty (-1)^n(2n+1)\alpha^{2n+1}\\ &=\cdots =\frac{\alpha}{2}\cdot\frac{d}{d\alpha}\left(\frac{\alpha}{1+\alpha^2}\right) =\frac{\alpha(1-\alpha^2)}{2\left(1+\alpha^2\right)^2}\\ [x^3]\arctan\left(\alpha e^x\right)&=\frac{1}{6}\sum_{n=0}^\infty (-1)^n(2n+1)^2\alpha^{2n+1}\\ &=\frac{\alpha^2}{6}\sum_{n=0}^\infty (-1)^n(2n+1)(2n)\alpha^{2n-1} +\frac{\alpha}{6}\sum_{n=0}^\infty (-1)^n(2n+1)\alpha^{2n}\\ &=\cdots =\left(\frac{\alpha^2}{6}\cdot\frac{d^2}{d\alpha^2}+\frac{\alpha}{6}\cdot\frac{d}{d\alpha}\right)\left(\frac{\alpha}{1+\alpha^2}\right)\\ &=\cdots =\frac{\alpha^5-6\alpha^3+\alpha}{6(1+\alpha^2)^3} \end{align*}$$

Wnioskujemy

$$\begin{align*} \arctan\left(\alpha e^x\right)&=\arctan(\alpha)+\frac{\alpha}{1+\alpha^2}x+\frac{\alpha(1-\alpha^2)}{2\left(1+\alpha^2\right)^2}x^2\\ &\qquad+\frac{\alpha^5-6\alpha^3+\alpha}{6(1+\alpha^2)^3}x^3+O(x^4)\tag{3} \end{align*}$$

Moce w szeregach logarytmicznych:

W celu uzyskania współczynników szeregu logarytmicznego

$$\begin{align*} \ln\left(\arctan(\alpha e^x)\right)&=\sum_{n=1}^\infty\frac{(-1)^{n+1}}{n}(\arctan(\alpha e^x)-1)^n \end{align*}$$ piszemy wyrażenie (3) jako $$\begin{align*} \arctan\left(\alpha e^x\right)&=a_0+a_1x+a_2x^2+a_3x^3+O(x^4) \end{align*}$$ i rozważamy $$\begin{align*} \ln\left(\arctan(\alpha e^x)\right)&=\sum_{n=1}^\infty\frac{(-1)^{n+1}}{n}((a_0-1)+a_1x+a_2x^2+a_3x^3)^n+O(x^4)\tag{4} \end{align*}$$

Teraz ustawiamy $A(x)=(a_0-1)+a_1x+a_2x^2+a_3x^3$ i wyodrębnij współczynniki $x^0$ do $x^3$ od

$$\begin{align*} \left(A(x)\right)^n&=((a_0-1)+a_1x+a_2x^2+a_3x^3)^n\\ &=\sum_{j=0}^n\binom{n}{j}(a_0-1)^j\left(a_1x+a_2x^2+a_3x^3\right)^{n-j}\\ &=\sum_{j=0}^n\binom{n}{j}(a_0-1)^j\sum_{k=0}^{n-j}\binom{n-j}{k}a_1^kx^k\left(a_2x^2+a_3x^3\right)^{n-j-k}\tag{5}\\ \end{align*}$$

Uzyskujemy od (5)

$$\begin{align*} [x^0]&\left(A(x)\right)^n=\cdots=(a_0-1)^n\\ [x^1]&\left(A(x)\right)^n=\cdots=a_1n(a_0-1)^{n-1}\\ [x^2]&\left(A(x)\right)^n=\dots=a_2n(a_0-1)^{n-1}+\frac{1}{2}n(n-1)a_1^2(a_0-1)^{n-2}\\ [x^3]&\left(A(x)\right)^n=\cdots=na_3(a_0-1)^{n-1}+a_1a_2n(n-1)(a_0-1)^{n-2}\\ &\qquad\qquad\qquad\qquad+\frac{1}{6}n(n-1)(n-2)a_1^3(a_0-1)^{n-3}\tag{6} \end{align*}$$

Rozszerzanie szeregowe logarytmu:

Obliczamy przy użyciu (6) współczynników $\ln(\arctan(\alpha e^x))$ pod względem $a_j, 0\leq j\leq 3$

$$\begin{align*} [x^0]&\ln(\arctan(\alpha e^x))=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^0]A(x)\\ &=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^0](a_0-1)^n\\ &=\ln(a_0-1)\\ [x^1]&\ln(\arctan(\alpha e^x))=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^1]A(x)\\ &=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^0]a_1n(a_0-1)^{n-1}\\ &=a_1\sum_{n=0}^\infty(-1)^n(a_0-1)^n\\ &=\frac{a_1}{a_0}\\ [x^2]&\ln(\arctan(\alpha e^x))=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^2]A(x)\\ &=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}\left(a_2n(a_0-1)^{n-1}+\frac{1}{2}n(n-1)a_1^2(a_0-1)^{n-2}\right)\\ &=\cdots\\ &=\left(a_2+\frac{a_1^2}{2}\frac{d}{da_0}\right)\left(\frac{1}{a_0}\right)\\ &=\frac{a_2}{a_0}-\frac{a_1^2}{2a_0^2}\\ [x^3]&\ln(\arctan(\alpha e^x))=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^3]A(x)\\ &=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}\left(na_3(a_0-1)^{n-1}+a_1a_2n(n-1)(a_0-1)^{n-2}\right.\\ &\qquad\qquad\qquad\qquad\left.+\frac{1}{6}n(n-1)(n-2)a_1^3(a_0-1)^{n-3}\right)\\ &=\cdots\\ &=\left(a_3+a_1a_2\frac{d}{da_0}+\frac{a_1^3}{6}\frac{d^2}{da_0^2}\right)\left(\frac{1}{a_0}\right)\\ &=\frac{a_3}{a_0}-\frac{a_1a_2}{a_0^2}+\frac{a_1^3}{3a_0^3}\tag{7} \end{align*}$$

Rozszerzenie serii $f(x)$:

Teraz nadszedł czas na żniwa. W końcu uzyskujemy z (3) i (7) szanując to$f(x)$ to jest dziwne

$$\begin{align*} \color{blue}{f(x)}&\color{blue}{=\ln\left(\arctan\left(\alpha e^x\right)\right)-\ln\left(\arctan\left(\alpha e^{-x}\right)\right)}\\ &=\cdots\\ &=\frac{2a_1}{a_0}x+2\left(\frac{a_3}{a_0}-\frac{a_1a_2}{a_0^2}+\frac{a_1^3}{3a_0^3}\right)x^3+O(x^5)\\ &\color{blue}{=\frac{2\alpha}{(1+\alpha^2)\arctan(\alpha)}x +\frac{\alpha}{3(1+\alpha^2)^3\arctan(\alpha)}\Bigg(\alpha^4-6\alpha^2+1}\\ &\qquad\color{blue}{\left.-\frac{3\alpha(1-\alpha^2)}{\arctan(\alpha)}+\frac{2\alpha^2}{\left(\arctan(\alpha)\right)^2}\right)x^3+O(x^5)} \end{align*}$$

(Konwertowanie komentarza na odpowiedź).

Używając Wolfram Alpha, $f'''(x)$ w $x=0$ ocenia:

$$\begin{align} \\ f'''(0) = & -\frac{6 \alpha^2}{(\alpha^2 + 1)^2 (\arctan(\alpha))^2} \ + \ \frac{2 \alpha}{(\alpha^2 + 1) \arctan(\alpha)} \\ & \quad \quad + \frac{16 \alpha^5}{(\alpha^2 + 1)^3 \arctan(\alpha)} \ + \ \frac{12 \alpha^4}{(\alpha^2 + 1)^3 (\arctan(\alpha))^2} \\ & \quad \quad \quad \quad - \frac{16 \alpha^3}{(\alpha^2 + 1)^2 \arctan(\alpha)} \ + \ \frac{4 \alpha^3}{(\alpha^2 + 1)^3 (\arctan(\alpha))^3} \\ \\ &= \frac{2 (\alpha^4 - 6 \alpha^2 + 1) \alpha}{(\alpha^2 + 1)^3 \arctan(\alpha)} + \frac{6 (\alpha^2 - 1) \alpha^2}{(\alpha^2 + 1)^3 (\arctan(\alpha))^2} + \frac{4 \alpha^3}{(\alpha^2 + 1)^3 (\arctan(\alpha))^3} \\ \end{align}$$

http://www.wolframalpha.com/input/?i=evaluate+d3%2Fdx3++(+ln+(arctan+(a+exp(+x)))+-+ln+(arctan(a+exp(-+x) )) +) + w + x% 3D0

Możemy również dwukrotnie sprawdzić, czy odpowiada to odpowiedzi Markusa tutaj .

Jego współczynnik $x^3$ wydaje się być

$$\frac{\alpha}{3(1+\alpha^2)^3\arctan(\alpha)}\left(\alpha^4-6\alpha^2+1-\frac{3\alpha(1-\alpha^2)}{\arctan(\alpha)}+\frac{2\alpha^2}{\left(\arctan(\alpha)\right)^2}\right).$$

Jeśli pomnożymy to przez 6 i przestawimy niektóre czynniki, otrzymamy:

$$\frac{2\alpha(\alpha^4-6\alpha^2+1)}{(1+\alpha^2)^3\arctan(\alpha)} - \frac{6\alpha^2(1-\alpha^2)}{(1+\alpha^2)^3(\arctan(\alpha))^2} + \frac{4 \alpha^3}{(1+\alpha^2)^3 (\arctan(\alpha))^3}$$

który pasuje!

Problem postawiony w pytaniu wydaje się nie mieć rozwiązania w formie zamkniętej. Jak wspomniano w pytaniu i pokazano w innych odpowiedziach, wynik można przekształcić w serię, którą można osiągnąć za pomocą dowolnego symbolicznego narzędzia matematycznego, takiego jak Mathematica. Jednak warunki stają się dość skomplikowane i brzydkie i nie jest jasne, jak dobre jest to przybliżenie, jeśli uwzględniamy warunki do trzeciego rzędu. Ponieważ nie możemy uzyskać dokładnej formuły, może być lepiej obliczyć rozwiązanie numerycznie, co w przeciwieństwie do przybliżenia da (prawie) dokładny wynik.

Jednak nie o to chodzi w mojej odpowiedzi. Proponuję inną drogę, która daje dokładne rozwiązanie, zmieniając sformułowanie problemu. Po zastanowieniu się przez chwilę okazuje się, że jest to specyfikacja częstotliwości środkowej$\omega_0$oraz określenie szerokości pasma jako stosunku (lub równoważnie w oktawach), co powoduje matematyczną trudność. Istnieją dwa sposoby wyjścia z dylematu:

1. określić szerokość pasma filtra czasu dyskretnego jako różnicę częstotliwości$\Delta\omega=\omega_2-\omega_1$, gdzie $\omega_1$ i $\omega_2$ są odpowiednio dolną i górną krawędzią pasmową filtra czasu dyskretnego.
2. przepisz stosunek $\omega_2/\omega_1$i zamiast $\omega_0$ określ jedną z dwóch częstotliwości krawędziowych $\omega_1$ lub $\omega_2$.

W obu przypadkach możliwe jest proste rozwiązanie analityczne. Ponieważ pożądane jest określenie szerokości pasma filtra czasu dyskretnego jako stosunku (lub, równoważnie, w oktawach), opiszę drugie podejście.

Zdefiniujmy częstotliwości brzegowe $\Omega_1$ i $\Omega_2$ filtra czasu ciągłego według

$$|H(j\Omega_1)|^2=|H(j\Omega_2)|^2=\frac12\tag{1}$$

z $\Omega_2>\Omega_1$, gdzie $H(s)$ jest funkcją przesyłania filtra pasmowego drugiego rzędu:

$$H(s)=\frac{\Delta\Omega s}{s^2+\Delta\Omega s+\Omega_0^2}\tag{2}$$

z $\Delta\Omega=\Omega_2-\Omega_1$, i $\Omega_0^2=\Omega_1\Omega_2$. Zauważ, że$H(j\Omega_0)=1$, i $|H(j\Omega)|<1$ dla $\Omega\neq\Omega_0$.

Używamy transformacji dwuliniowej do mapowania częstotliwości krawędzi $\omega_1$ i $\omega_2$ filtra czasu dyskretnego na częstotliwości brzegowe $\Omega_1$ i $\Omega_2$filtra czasu ciągłego. Bez utraty ogólności możemy wybrać$\Omega_1=1$. Dla naszych celów transformacja dwuliniowa przyjmuje wówczas postać

$$s = \frac{1}{\tan\left(\frac{\omega_1}{2}\right)}\frac{z-1}{z+1}\tag{3}$$

odpowiadające następującej zależności między częstotliwościami w czasie ciągłym a częstotliwościami w czasie dyskretnym:

$$\Omega=\frac{\tan\left(\frac{\omega}{2}\right)}{\tan\left(\frac{\omega_1}{2}\right)}\tag{4}$$

Od $(4)$ otrzymujemy $\Omega_2$ przez ustawienie $\omega=\omega_2$. Z$\Omega_1=1$ i $\Omega_2$ obliczone z $(4)$, otrzymujemy funkcję przesyłania analogowego filtra prototypowego od $(2)$. Zastosowanie transformacji dwuliniowej$(3)$, otrzymujemy funkcję transferu filtru pasmowego z czasem dyskretnym:

$$H_d(z)=g\cdot\frac{z^2-1}{z^2+az+b}\tag{5}$$

z

$$\begin{align}g&=\frac{\Delta\Omega c}{1+\Delta\Omega c+\Omega_0^2c^2}\\a&=\frac{2(\Omega_0^2c^2-1)}{1+\Delta\Omega c+\Omega_0^2c^2}\\b&=\frac{1-\Delta\Omega c+\Omega_0^2c^2}{1+\Delta\Omega c+\Omega_0^2c^2}\\c&=\tan\left(\frac{\omega_1}{2}\right)\end{align}\tag{6}$$

Podsumowanie:

Szerokość pasma filtru z czasem dyskretnym można określić w oktawach (lub ogólnie jako stosunek), a parametry analogowego filtra prototypowego można dokładnie obliczyć, tak aby osiągnąć określoną szerokość pasma. Zamiast częstotliwości środkowej$\omega_0$, określamy krawędzie pasma $\omega_1$ i $\omega_2$. Częstotliwość środkowa zdefiniowana przez$|H_d(e^{j\omega_0})|=1$ jest wynikiem projektu.

Niezbędne kroki są następujące:

1. Określ pożądany stosunek krawędzi pasma $\omega_2/\omega_1$i jedną z krawędzi pasma (co oczywiście jest równoznaczne z prostym określeniem) $\omega_1$ i $\omega_2$).
2. Wybierać $\Omega_1=1$ i ustal $\Omega_2$ od $(4)$. Obliczać$\Delta\Omega=\Omega_2-\Omega_1$ i $\Omega_0^2=\Omega_1\Omega_2$ analogowego filtra prototypowego $(2)$.
3. Oceń stałe $(6)$ aby uzyskać funkcję przenoszenia w czasie dyskretnym $(5)$.

Zauważ, że przy bardziej powszechnym podejściu gdzie $\omega_0$ i $\Delta\omega=\omega_2-\omega_1$ są określone, rzeczywiste krawędzie pasma $\omega_1$ i $\omega_2$są wynikiem procesu projektowania. W proponowanym rozwiązaniu można określić krawędzie pasma i$\omega_0$jest wynikiem procesu projektowania. Zaletą tego drugiego podejścia jest to, że szerokość pasma można określić w oktawach, a rozwiązanie jest dokładne, tzn. Wynikowy filtr ma dokładnie określoną szerokość pasma w oktawach.

Przykład:

Określmy szerokość pasma jednej oktawy i wybieramy dolną krawędź pasma jako $\omega_1=0.2\pi$. Daje to górną krawędź pasma$\omega_2=2\omega_1=0.4\pi$. Krawędzie pasma analogowego filtra prototypowego to$\Omega_1=1$ i od $(4)$ (z $\omega=\omega_2$) $\Omega_2=2.2361$. To daje$\Delta\Omega=\Omega_2-\Omega_1=1.2361$ i $\Omega_0^2=\Omega_1\Omega_2=2.2361$. Z$(6)$ dostajemy za funkcję transferu w czasie dyskretnym $(5)$

$$H_d(z)=0.24524 \cdot \frac{z^2-1}{z^2-0.93294z+0.50953}$$

która osiąga dokładnie szerokość pasma 1 oktawę i określone krawędzie pasma, jak pokazano na poniższym rysunku:

Numeryczne rozwiązanie pierwotnego problemu:

Z komentarzy rozumiem, że ważne jest, aby móc dokładnie określić częstotliwość środkową $\omega_0$ dla którego $|H_d(e^{j\omega_0})|=1$jest spełniony. Jak wspomniano wcześniej, nie jest możliwe uzyskanie dokładnego rozwiązania w formie zamkniętej, a rozwój serii daje dość niewygodne wyrażenia.

Dla jasności chciałbym podsumować możliwe opcje z ich zaletami i wadami:

1. podaj pożądaną szerokość pasma jako różnicę częstotliwości $\Delta\omega=\omega_2-\omega_1$i określ $\omega_0$; w tym przypadku możliwe jest proste rozwiązanie w formie zamkniętej.
2. określ krawędzie pasma $\omega_1$ i $\omega_2$(lub równoważnie szerokość pasma w oktawach i jedna z krawędzi pasma); prowadzi to również do prostego rozwiązania w formie zamkniętej, jak wyjaśniono powyżej, ale częstotliwości środkowej$\omega_0$ jest wynikiem projektu i nie można go określić.
3. podaj żądaną szerokość pasma w oktawach i częstotliwość środkową $\omega_0$ (as asked in the question); no closed form solution is possible, nor is there (for the time being) any simple approximation. For this reason I think it's desirable to have a simple and efficient method for obtaining a numerical solution. This is what is explained below.

When $\omega_0$ is specified we use a form of the bilinear transform with a normalization constant that is different from the one used in $(3)$ and $(4)$:

$$\Omega=\frac{\tan\left(\frac{\omega}{2}\right)}{\tan\left(\frac{\omega_0}{2}\right)}\tag{7}$$

We define $\Omega_0=1$. Denote the specified ratio of band edges of the discrete-time filter as

$$r=\frac{\omega_2}{\omega_1}\tag{8}$$

With $c=\tan(\omega_0/2)$ we get from $(7)$ and $(8)$

$$r=\frac{\arctan(c\Omega_2)}{\arctan(c\Omega_1)}\tag{9}$$

With $\Omega_1\Omega_2=\Omega_0^2=1$, $(9)$ can be rewritten in the following form:

$$f(\Omega_1)=r\arctan(c\Omega_1)-\arctan\left(\frac{c}{\Omega_1}\right)=0\tag{10}$$

For a given value of $r$ this equation can be solved for $\Omega_1$ with a few Newton iterations. For this we need the derivative of $f(\Omega_1)$:

$$f'(\Omega_1)=c\left(\frac{r}{1+c^2\Omega_1^2}+\frac{1}{c^2+\Omega_1^2}\right)\tag{11}$$

With $\Omega_0=1$, we know that $\Omega_1$ must be in the interval $(0,1)$. Even though it's possible to come up with smarter initial solutions, it turns out that the initial guess $\Omega_1^{(0)}=0.1$ works well for most specs, and will result in very accurate solutions after only $4$ iterations of Newton's method:

$$\Omega_1^{(n+1)}=\Omega_1^{(n)}-\frac{f(\Omega_1^{(n)})}{f'(\Omega_1^{(n)})}\tag{12}$$

With $\Omega_1$ obtained with a few iterations of $(12)$ we can determine $\Omega_2=1/\Omega_1$ and $\Delta\Omega=\Omega_2-\Omega_1$, and and we use $(5)$ and $(6)$ to compute the coefficients of the discrete-time filter. Note that the constant $c$ is now given by $c=\tan(\omega_0/2)$.

Example 1:

Let's specify $\omega_0=0.6\pi$ and a bandwidth of $0.5$ octaves. This corresponds to a ratio $r=\omega_2/\omega_1=2^{0.5}=\sqrt{2}=1.4142$. With an initial guess of $\Omega_1=0.1$, $4$ iterations of Newton's method resulted in a solution $\Omega_1=0.71$, from which the coefficients of the discrete-time can be computed as explained above. The figure below shows the result:

The filter was calculated with this Matlab/Octave script:

% specifications
bw = 0.5;    % desired bandwidth in octaves
w0 = .6*pi;  % resonant frequency

r = 2^(bw);  % ratio of band edges
W1 = .1;     % initial guess (works for most specs)
Nit = 4;     % # Newton iterations
c = tan(w0/2);

% Newton
for i = 1:Nit,
    f = r*atan(c*W1) - atan(c/W1);
    fp = c * ( r/(1+c^2*W1^2) + 1/(c^2+W1^2) );
    W1 = W1 - f/fp
end

W1 = abs(W1);
if (W1 >= 1), error('Failed to converge. Reduce value of initial guess.'); end

W2 = 1/W1;
dW = W2 - W1;

% discrete-time filter
scale = 1 + dW*c + W1*W2*c^2;
b = ( dW*c/scale) * [1,0,-1];
a = [1, 2*(W1*W2*c^2-1)/scale, (1-dW*c+W1*W2*c^2)/scale ];

Example 2:

Dodaję kolejny przykład, aby pokazać, że ta metoda może również poradzić sobie ze specyfikacjami, w przypadku których większość przybliżeń da niesensowne wyniki. Często dzieje się tak, gdy pożądana szerokość pasma i częstotliwość rezonansowa są duże. Zaprojektujmy filtr z$\omega_0=0.95\pi$ i $bw=4$oktawy. Cztery iteracje metody Newtona z wstępnym odgadnięciem$\Omega_1^{(0)}=0.1$ dają końcową wartość $\Omega_1=0.00775$, tj. w paśmie analogowego prototypu $\log_2(\Omega_2/\Omega_1)=\log_2(1/\Omega_1^2)\approx 14$oktawy. Odpowiedni filtr czasu dyskretnego ma następujące współczynniki, a jego odpowiedź częstotliwościowa pokazano na wykresie poniżej:

b = 0,90986 * [1,0, -1];
a = [1,00000 0,17806 -0,81972];

Wynikowe krawędzie połowy mocy są $\omega_1=0.062476\pi$ i $\omega_2 = 0.999612\pi$, które są rzeczywiście dokładnie $4$ oktawy (tj. współczynnik wynoszący $16$) osobno.

okej, obiecałem wystawić nagrodę i dotrzymam obietnicy. ale muszę wyznać, że mógłbym się trochę oderwać od zadowolenia z trzeciej pochodnej$f(x)$. tak naprawdę chcę dwa współczynniki$g(y)$.

więc nie zdawałem sobie sprawy z tego, że istnieje ten język Wolfram jako alternatywa dla matematyki lub Derive, i nie zdawałem sobie sprawy, że może tak łatwo obliczyć trzecią pochodną i uprościć wyrażenie.

a ten facet Markus z matematyki SE opublikował tę odpowiedź (która, jak sądzę, będzie musiała być ilością grunge, która moim zdaniem będzie potrzebna).

$$\begin{align*} y = f(x) &= \ln\left(\arctan\left(\alpha e^x\right)\right) - \ln\left(\arctan\left(\alpha e^{-x}\right)\right)\\ \\ &\approx a_1 x \ + \ a_3 x^3 \\ \\ &=\frac{2\alpha}{(1+\alpha^2)\arctan(\alpha)} x + \frac{\alpha}{3(1+\alpha^2)^3\arctan(\alpha)}\Bigg(\alpha^4-6\alpha^2+1\\ &\qquad\left.-\frac{3\alpha(1-\alpha^2)}{\arctan(\alpha)}+\frac{2\alpha^2}{\left(\arctan(\alpha)\right)^2}\right) x^3 \\ \\ \end{align*}$$

so i put together the third-order approximation to the inverse:

$$\begin{align*} x = g(y) &\approx b_1 y \ + \ b_3 y^3 \\ \\ &= \frac{1}{a_1} y \ - \ \frac{a_3}{a_1^4} y^3 \\ \\ &= \frac{(1+\alpha^2)\arctan(\alpha)}{2\alpha} y - \frac{(1+\alpha^2)(\arctan(\alpha))^3}{48 \alpha^3}\Bigg(\alpha^4-6\alpha^2+1 \\ &\qquad\left.-\frac{3\alpha(1-\alpha^2)}{\arctan(\alpha)}+\frac{2\alpha^2}{\left(\arctan(\alpha)\right)^2}\right) y^3 \\ \\ &= \frac{(1+\alpha^2)\arctan(\alpha)}{2\alpha} y - \frac{(1+\alpha^2)(\arctan(\alpha))^3}{48 \alpha}\Bigg(\alpha^2-6+\alpha^{-2} \\ &\qquad\left.-\frac{3(1-\alpha^2)}{\alpha\arctan(\alpha)}+\frac{2}{\left(\arctan(\alpha)\right)^2}\right) y^3 \\ \\ &= y \left(\arctan(\alpha)\frac{\alpha + \alpha^{-1}}{2}\right) \Bigg(1 \ + \\ & \left((\arctan(\alpha))^2\left(1-\frac{\alpha^2+\alpha^{-2}}{6} \right) - \arctan(\alpha)\frac{\alpha-\alpha^{-1}}{2}-\frac{1}{3}\right) \frac{y^2}{4} \Bigg) \\ \\ \end{align*}$$

i was kinda hoping someone else would do this. recall $y=f(x) \triangleq \ln(2) \, bw$, $g(y) = x \triangleq \frac{\ln(2)}{2}BW$ and $\alpha \triangleq \tan(\omega_0/2)$

$$\begin{align*} x = g(y) &\approx y \left(\arctan(\alpha)\frac{\alpha + \alpha^{-1}}{2}\right) \Bigg(1 \ + \\ & \left((\arctan(\alpha))^2\left(1-\frac{\alpha^2+\alpha^{-2}}{6} \right) - \arctan(\alpha)\frac{\alpha-\alpha^{-1}}{2}-\frac{1}{3}\right) \frac{y^2}{4} \Bigg) \\ \\ \frac{\ln(2)}{2}BW &\approx (\ln(2) bw) \left(\arctan(\alpha)\frac{\alpha + \alpha^{-1}}{2}\right) \Bigg(1 \ + \\ & \left((\arctan(\alpha))^2\left(1-\frac{\alpha^2+\alpha^{-2}}{6} \right) - \arctan(\alpha)\frac{\alpha-\alpha^{-1}}{2}-\frac{1}{3}\right) \frac{(\ln(2) bw)^2}{4} \Bigg) \\ \\ \end{align*}$$

i have three convenient trig identities:

$$\frac12 \left(\alpha + \alpha^{-1} \right) = \frac12 \left(\tan(\omega_0/2) + \frac{1}{\tan(\omega_0/2)} \right) = \frac{1}{\sin(\omega_0)} \\ $$

$$\frac12 \left(\alpha - \alpha^{-1} \right) = \frac12 \left(\tan(\omega_0/2) - \frac{1}{\tan(\omega_0/2)} \right) = -\frac{1}{\tan(\omega_0)} \\ $$

$$\frac12 \left(\alpha^2 + \alpha^{-2} \right) = \frac12 \left(\tan^2(\omega_0/2) + \frac{1}{\tan^2(\omega_0/2)} \right) = \frac{1}{\sin^2(\omega_0)} + \frac{1}{\tan^2(\omega_0)} \\ = \frac{2}{\sin^2(\omega_0)} - 1$$

"finally" we got:

$$BW \approx bw \, \frac{\omega_0}{\sin(\omega_0)} \left(1 \ + \ \frac{(\ln(2))^2}{24} \left( 2(\omega_0^2 - 1) - \left(\frac{\omega_0}{\sin(\omega_0)}\right)^2 + 3\frac{\omega_0}{\tan(\omega_0)} \right) (bw)^2 \right)$$

this ain't so bad. fits on a single line. if someone sees an error or a good way to simplify further, please lemme know.

with the power series approximation from the comment above,

$$BW \approx bw \, \frac{\omega_0}{\sin(\omega_0)} \left(1 \ + \ (\ln(2))^2 \left( \tfrac{1}{36}\omega_0^2 - \tfrac{1}{180}\omega_0^4 - \tfrac{2}{2835}\omega_0^6 \right) (bw)^2 \right)$$

oto kilka wyników ilościowych. zaplanowałem określoną przepustowość$bw$dla filtra cyfrowego na osi x i wynikowej cyfrowej przepustowości na osi y. jest pięć wykresów od zielonego do czerwonego reprezentujących częstotliwość rezonansową$\omega_0$ znormalizowany przez Nyquist:

$\frac{\omega_0}{\pi} =$ [0,0002 0,2444 0,4880 0,7320 0,9759]

więc częstotliwość rezonansowa przechodzi od prawie stałego prądu stałego do prawie Nyquista.

tutaj nie ma żadnej rekompensaty (lub wstępnego wypaczenia) za przepustowość:

oto prosta rekompensata za pierwsze zamówienie, którą książka kucharska zrobiła przez cały czas:

oto odszkodowanie trzeciego rzędu, które właśnie rozwiązaliśmy tutaj:

chcemy, aby wszystkie linie leżały bezpośrednio na głównej przekątnej.

I nie pomylił się w przypadku trzeciego rzędu i skorygować go w tej wersji. to ma wyglądać zbliżenia do trzeciego rzędu$g(y)$ jest nieco lepszy niż przybliżenie pierwszego rzędu dla małych $bw$.

więc podważyłem współczynnik terminu trzeciego rzędu (chcę pozostawić ten sam termin pierwszego rzędu), zmniejszając jego efekt. wynika to z pomnożenia tylko terminu trzeciego rzędu przez 50%:

zmniejsza to do 33%:

a to zmniejsza termin trzeciego rzędu do 25%:

ponieważ celem funkcji odwrotnej jest cofnięcie określonej funkcji, celem tego wszystkiego jest doprowadzenie krzywych funkcji złożonej do położenia jak najbliżej głównej przekątnej. nie jest tak źle nawet do 75% Nyquista dla częstotliwości rezonansowej$\omega_0$ i przepustowość 3 oktaw $bw$. ale nie o wiele lepiej, aby naprawdę warto było to zrobić w kodzie „Współczynnik gotowania”, który jest wykonywany za każdym razem, gdy użytkownik obraca pokrętło lub przesuwa suwak.