通过在参数曲线上分布极点和零点来进行滤波器设计

#1 楼

在整个答案中，我将使用数学符号表示，即在频域中表示滤波器幅度响应的数学等价形式。为此，将使用$ x $代替$ j \ omega $，以更好地反映@Olli关于寻找数学参数曲线以近似滤波器的问题。由于这不是滤波器设计，因此拐角频率被归一化为单位，因此是$ x $而不是$ \ omega / \ omega_p $。


我不确定这是否是您正在寻找答案，但是任何过滤器都可以通过通用传递函数来表示：

$$ H ^ 2（x）= \ frac {1} {1+ \ epsilon_p ^ 2 R ^ 2（x）} $$

其中$ \ epsilon_p = \ sqrt {10 ^ {A_p / 10} -1} $和$ R（x）$是特征衰减函数。 $ A_p $是通带衰减/纹波，单位为dB，但对于Cauer / Elliptic，反Pascal或反Chebyshev（也称为“ Chebyshev Type II”），也可以位于阻带中。后者表示为：

$$ H ^ 2（x）= \ frac {1} {1+ \ displaystyle \ frac {1} {\ epsilon_s ^ 2 T_N ^ 2（x）} } $$

对于巴特沃思来说，如您所见：

$$ R（x）= x ^ N $$

对于Chebyshev它是$ R（x）= T_N（x）$或Chebyshev多项式（$ x \ leq1 $的$ \ cos $ / $ \ rm acos $和$ x> 1的$ \ cosh $ / $ \ rm acosh $ $），对于Elliptic，它是：

$$ R（x）= \ mathrm {cd} \ left（N \ frac {K_1} {K} \ mathrm {cd} ^ {-1}（ x，k），k_1 \ right）$$

在Lutovac的书中，通过椭圆过滤器的精确等效函数，有一些极其简单的表示形式。例如，可以通过以下方式精确表示二阶传递函数：

$$ R（x）= \ frac {\ left（\ sqrt {1-k ^ 2} +1 \ right）x ^ 2-1} {\ left（\ sqrt {1-k ^ 2} -1 \ right）x ^ 2 + 1} $$

其中唯一的依赖关系是模数$ k $ 。

这些是已知类型，对于不太为人所知的类型，例如，Legendre，$ R（x）= P_N（x）$，其中$ P_N（x）$是Legendre多项式，对于Pascal滤波器是Pascal多项式的平移和归一化版本，它是：

$$ \ binom {\ frac {N + 1} {2} x + \ frac {N-1} {2}} {N} $$

清单继续。有些近似地有所不同，例如，高斯是$ \ lvert H（x）\ rvert ^ 2 = \ exp（-x ^ 2）$，用MacLaurin系列进行了扩展，对于Bessel来说也是一样，拉普拉斯表达式$ \ exp（-s）$的分母形式为：

$$ a_i = \ frac {（2N-1）！} {2 ^ {Ni} i！（Ni）！ } $$

还有更多奇特的方法来推导传递函数，例如Papoulis（Optimum L）和Halpern，它们都使用Legendre多项式来积分响应，从而传递函数高滤光片选择性单调减小。对于Papoulis，它是：

$$ R（x ^ 2）= \ int_ {i = -1} ^ {2x ^ 2-1} {\ left（\ sum_ {i = 0} ^ k {a_i \ cdot P_i（x）} \ right）^ 2} $$

其中$ k $是$ \ lfloor（N-1）/ 2 \ rfloor $，而$ a_i $是一些巧妙选择的术语，取决于$ N $或$ k $两者都是奇/偶。

如上所述，所有这些都不使用频域来表示，如$ x $是数学的$ x $，是实数，而不是虚构的$ j \ omega $。求解根可以通过在用$ j \ omega $替换$ x $时简单地找到传递函数的极点（和零点）来完成，从而找出$ H（s）H（-s）$并选择Hurwitz多项式，或简单地在$ x $中找到数学表达式的根（请参见下面第二条注释中的链接）。这将产生根旋转90度的根，这意味着要做的就是在它们之间切换实部和虚部，然后选择右手侧。


我想在这一点上，重要的是要说不存在过滤器，因为人们在地图上投掷飞镖以标出极点，在仔细考虑他们的目标之后。

例如，随着质量的不断提高，巴特沃斯滤波器的出现是因为需要一种设计简单，衰减单调增加的滤波器。 Linkwitz-Riley只是伪装成巴特沃斯而已，巧妙地掩盖了具有相同转折频率的低通和高通，从而得到平坦的响应，对音频应用很有用。被设计为具有更好的衰减，但以通带或阻带的纹波为代价。勒让德（Legendre），超球形，帕斯卡（Pascal）（以及其他）可以最大程度地减小纹波，从而改善群时延，但会以稍微降低衰减为代价。 ，同时改善了拐角频率附近的衰减，但以通带下垂为代价。

Cauer /椭圆滤波器在通带和阻带中都使用了纹波，以便将所需的阶数最小化相同或更好的衰减。

所有这些都在频域中，大多数滤波器都在频域中。换句话说，贝塞尔滤波器是由于需要近似模拟延迟，因此随着阶次增加，它们趋向于向\\ exp（-j \ omega）$收敛，而高斯滤波器是为零过冲而创建的，因此它们约$ \ exp（-x ^ 2）$随顺序增加。

当然，正如有人建议的那样，您也可以撒上杆子，看看​​结果如何，可以将它们配置为星形或某种蜂窝状图案，选择自己喜欢的lemniscate，但这不是您想要过滤器的方法出来。当然，您可能会得到一个异乎寻常的响应，甚至可能适用于谁知道谁，在百万分之一的情况下，但这实际上只是一个特殊的情况。可行的方法是首先提出一个设计目标，然后看看如何通过物理上可实现的滤波器实现该目标。即使这意味着想出一个“谁知道哪里”适用的过滤器。 :-)


鉴于@Olli的最新答案，请考虑设计用于0.9@fp=1, 0.1@fs=5的Butterworth滤波器的简单情况。计算是这样的：

\ begin {align}
Ap＆=-20 \ log_ {10}（0.9）= 0.91515 \ textrm {dB} \\
As＆ = -20 \ log_ {10}（0.1）= 20 \ textrm {dB} \\
\ epsilon_p＆= \ sqrt {10 ^ {A_p / 10} -1} = 0.48432 \\
\ epsilon_s＆ = \ sqrt {10 ^ {A_s / 10} -1} = 9.94987 \\
F＆= \ frac {f_s} {f_p} = \ frac 51 = 5 \\
N＆= \ frac {\ log {\ frac {\ epsilon_s} {\ epsilon_p}}} {\ log {F}} = 1.878
\ end {align}

$ N $被四舍五入计算，因此$ N = 2 $。这意味着，如果将滤波器的响应与通带匹配，则阻带@fs的衰减会更高。使用第一个公式，衰减@fs为：

$$ H（f_s）= \ frac {1} {\ sqrt {1 + 0.48432 ^ 2 * 5 ^ {2N}}} == 0.08231 <0.1 $$

如果必须匹配阻带以具有0.1@fs，则必须应用频率校正：

\ begin {align }
\ omega _ {\ text {scale}}＆= \ left（\ frac {\ epsilon_s} {\ epsilon_p} \ right）^ {1 / N} \ frac {f_p} {f_s} = 9.94987 ^ { 0.5} = 0.9065 \\
H（5 * \ omega _ {\ text {scale}}）＆= 0.1
\ end {align}

因此，$ \ omega _ {\ text {scale}} $的范围从$ 1 $到$ 0.9065 $不等，您将获得介于两个极端之间的所有无限可能。你能做到吗？是。这值得么？即使您可能找到一两个参数，一般的答案仍然是“否”。这怎么可能？因为已经获得了巴特沃思滤波器的初始响应，所以您事先知道了一个滤波器的解析表达式，该滤波器具有单调递减的频率衰减，从而导致从传递函数的分母中找出极点。


鉴于@Olli的最新答案，有一些事情需要说明。首先，所有这些都是关于滤波器设计的，无论您如何看待：从数学还是从物理可实现性的角度来看。

如果它是数学的，那么关于它的理论有一些有趣的部分，即从同一滤波器获得不同的阶数而无需重新设计原始滤波器。

但是从物理可实现性的角度来看，整个过程意味着一些额外的，不需要的工作，（应该）导致相同的结果，而这恰恰是关于增加/减少过滤器顺序以获得新过滤器的部分。一。我的论据如下。

任何滤波器的核心都是用来过滤不需要的频率，无论它们是电，机械还是其他物理量。他们的目的是修改频谱（或群时延或时间响应）。如果需要这样的设备，则不能仅通过放入任何类型的过滤器来设计该设备，“只要将其放在这里，它就会过滤掉东西”；它的设计通常涉及很多。但是所有这些过程必须从需求开始。也就是说，首先必须有一个特定的目标，“让我们过滤掉高于$ 100 \ textrm {Hz} $的所有内容”，或“让只有红外光通过”，或者类似的东西，首先要确定参数为作为一个简单的例子，如果需要滤除低于$ 300 \ textrm {Hz} $和高于$ 3000 \ textrm {Hz} $的频率，则不会如果只使用具有这些转折频率的带通滤波器，则还必须指定衰减，无论是需要还是接受通带或阻带中的波纹，或者相位是否线性，群延迟如何影响所有因此，首先，有一些需要设计过滤器的特定参数。

指定了参数之后，将如何设计过滤器？让我们假设需要一个12阶椭圆低通滤波器，并且有可能将一个低阶滤波器增加到一个高阶滤波器（请参阅@Olli的答案）。假设将4阶转换为12阶的过程是完美无缺的，有一种方法可以指定4阶滤波器的设计参数，从而在转换后最终得到12阶满足那些条件。如果您愿意的话，请使用“有预谋的思维”。随之而来的问题是：如何设计四阶滤波器？答案只能通过已知的设计方法来解决。并且，如果有其他方法要出现或尚未发明，则必须首先应用这些方法以设计该四阶滤波器。此后才能计算12阶。从一开始就假设，即使转换过程完美无瑕，也仅意味着最终滤波器（整个设计试图收敛到的12阶）需要两个设计步骤：一个是四阶，第二个，对于12阶，整个过程成为不必要的麻烦，因为首先可以简单地设计4阶方法使用12阶滤波器。

让我们再进一步一点，并假设更多。生成的12阶极点将位于椭圆上，零点位于虚轴上。它们之间的距离将由控制椭圆滤波器的基本椭圆函数精确定义。就像@Olli所希望的那样，假设有一种方法来定义这些曲线，这样就可以通过简单地使用这些（无论是否有参数）曲线从一开始就很容易地设计滤波器。电极放置完成。到现在为止还挺好。但是必须首先计算这些曲线，并且解开曲线所用的参数是用于滤波器设计的确切参数，这些参数可以通过已知或未知的其他方法生成滤波器。而且，仍然需要进行计算，并且最有可能的是，这些参数曲线的基本定义必须是椭圆形，一种或另一种形式，否则就不会出现椭圆过滤器[note＃1]。这意味着整个过程将仅仅是椭圆滤波器的另一种设计方法，因为椭圆滤波器的极点已经具有闭合的形式表达式。

不要误会我。如果可以以一种方式设计一个滤波器，则可以以另一种方式设计它。这只是“尚未被人们知道”的方式之一。发明家真是太好了。但是，如果这种设计方法隐含了额外的步骤才能收敛到与其他方法相同的结果，那么这似乎不是可行的方法。而且请注意：在谈论过滤器设计时，我不会使用名称或描述性标签，而只是使用通用名称，因为只要结果正确且方法不正确，那么使用哪种方法并不重要。 t妨碍了设计过程。

[note＃1]：仅遵循一条通用曲线来放置极点是不够的，我将举两个与Butterworth滤波器有关的示例，其中将极点放置在等距角度的圆上。 Chebyshev I型滤波器的极点以Butterworth的角度放置在椭圆上，但投影在虚轴上，直到它们与椭圆相交。修改极点之间的距离将导致非等距行为，从而使滤波器成为非切比雪夫类型。同样，最小Q椭圆滤波器的极点放置在下面的圆上，但这并不意味着它是巴特沃思（即使波纹对于椭圆滤波器而言是最小的），因为角度之间的距离不相等。对于最后一个，这里是两个8阶Butterworth和最小Q椭圆的比较：



总体上，尽管这个问题带来了真正的兴趣，但我担心它有只不过是一个理论值，充其量只是一个教育意义，因为它不能完全满足过滤器设计的要求。当然，如果它应该被证明具有实际价值，我很高兴被证明是错误的，因为这将意味着有一种新的滤波器设计方法，可能比已经存在的方法更好。

#2 楼

虽然我凭直觉感到自己了解要求的内容，但我还是很难表达。我不确定这是否是由于我自己的局限性，或者问题确实是困难的还是不适当地的。我感觉它不适。因此，这是我的尝试：


目标是构建过滤器。也就是说，计算一组有理形式的系数：

$$ H（s）= \ frac {B（s）} {A（s）} = \ frac {\ sum_ {m = 0} ^ {M} b_m \ cdot s ^ m} {s ^ N + \ sum_ {n = 0} ^ {N-1} a_n \ cdot s ^ n} $$

（请注意，它不必在s平面上，也可以在z平面上，而且可以考虑使用更简单的形式（例如$ H（s）$仅具有极点）。现在让我们在s平面上运行，让我们也保留提名者。）


数字滤波器的特征在于它们的频率和相位响应，这两者可以完全由$ a_n，b_m $系数的值（或s平面上的位置）。到目前为止，讨论似乎集中在频率响应上，因此让我们暂时考虑一下。
给定一组一些$ a_n，b_m $和一些点\\ sigma + j \ omega $平面上，得出该点频率响应的几何方法是形成“零向量”（从零的位置朝向特定点）和“极向量”（对于极点类似），求和它们的大小和形式该比率如上等式中所述。
要问“由某些最优性准则定义的哪些滤波器类型具有由参数曲线定义的无穷子集”，就是要问“ ...的一对一些参数曲线$ A（s，\ Theta），B（s，\ Theta）$，其位置也会产生幅度响应曲线，在$ \ Theta $上具有特定的所需特性（例如，斜率，波纹，其他）。 Theta $是...参数的参数。
此时，请注意：一方面，我们正在寻找满足两个约束的$ A（s），B（s）$。首先，它们必须满足参数的约束（简单），其次，必须满足幅度响应特性（困难）所指定的约束。 -摆姿势是因为除了直接评估之外，没有将频率响应约束与参数$ A（s，\ Theta），B（s，\ Theta）$连接的分析方法。换句话说，目前无法在频率响应曲线上指定一些约束，并且无法通过该约束向后工作并找到满足这些约束的参数。我们可以走另一条路，但不能倒退。
因此，目前，我（实际上认为）可以接受$ A（s，\ Theta），B（s，\ Theta）$的某种特定形式，然后检查它们的性能如何作为滤波器，或者迭代地在其参数允许的范围内反复移动其系数，以在它们的$ \ Theta的特定范围内挤压出它们可以提供的最佳性能。 $。但是，我们可能会发现，给定椭圆的特性（例如），给定的参数化迭代方案可能会选择“弯曲”系数，使其尽可能接近某些“椭圆”区域特性。这就是为什么我在前面提到我们可能会发现将一个复杂的参数分解为“椭圆和”或“具有已知特征的曲线和”的原因。也许在这里需要第三个约束，即阅读“远离已知的$ A（s），B（s）$配置”，换句话说，对开始看起来像椭圆的解决方案进行惩罚（但仍处于迭代方案中）。 br />
最后，如果到目前为止这条路径还不太错误，那么我们就接近于遗传算法，用于滤波器设计，或者其他一些明智的“暗中射击”技术，通过这些技术可以得出满足特定条件的滤波器系数， 。以上只是一个例子，有更多类似的出版物。

希望有帮助。

#3 楼

我认为Butterworth滤波器定义为在$ \ omega = 0 $时最大平坦的全极点滤波器（对于LPF原型，意味着$ | H（j \ omega）| $的最可能导数）并不是特别引人注目在$ \ omega = 0 $处为零），其s平面极点在半径$ \ omega_0 $的左半圆上等距分布。

从“最大平坦”和“否零”，则可以导出

$$ | H（j \ omega）| ^ 2 = \ frac {1} {1 + \ left（\ frac {\ omega} {\ omega_0} \ right ）^ {2N}} $$

用于$ N $阶Butterworth。

so

$$ | H（s） | ^ 2 = \ frac {1} {1 + \ left（\ frac {s} {j \ omega_0} \ right）^ {2N}} $$

$ s = p_n $是一个分母为零时的极点。

$$ 1 + \ left（\ frac {p_n} {j \ omega_0} \ right）^ {2N} = 0 $$

or

$$ \ left（\ frac {p_n} {j \ omega_0} \ right）^ {2N} = -1 $$

$$ p_n ^ { 2N} =-（j \ omega_0）^ {2N} $$

$$ | p_n | = \ omega_0 $$

$$ 2N \ cdot \ arg \ {p_n \} =-\ pi + 2N \ cdot \ frac {\ pi} {2} + 2 \ pi n $$

$$ \ arg \ {p_n \} = \ frac {\ pi} {2} + \ frac {\ pi} {N} \ left（n-\ tfrac12 \ right）$$


对于$ N $阶Tchebyshev（类型1，是全极点），它是这样的：

$$ | H（j \ omega） | ^ 2 = \ frac {1} {1 + \ epsilon ^ 2 T_N ^ 2 \ left（\ frac {\ omega} {\ omega_c} \ right）} $$

其中$$ T_N （x）\ triangleq \ begin {cases}
\ cos \ big（N \，\ arccos（x）\ big）和\ text {if} | x | \ le 1 \\
\ cosh \ big（N \，\ operatorname {arccosh}（x）\ big）和＆text {if} x \ ge 1 \\
（-1）^ N \，\ cosh \ big（N \，\ operatorname {arccosh}（-x）\ big）和\ text {if} x \ le -1
\ end {cases} $$

是$ N $阶Tchebyshev多项式并满足递归：

$$ \ begin {align}
T_0（x）＆= 1 \\
T_1（x）＆= x \\
T_ {n + 1}（x）＆= 2xT_n（x）-T_ {n-1}（x）\ quad \ quad \ forall n \ in \ mathbb {Z} \ ge 1
\ end {align} $$

，$ \ omega_c $是“通带截止”频率，不要与-3 dB频率$ \混淆omega_0 $。 （但两者是相关的。）

通带纹波参数为$ \ epsilon = \ sqrt {10 ^ {\ tfrac {dB_ \ text {ripple}} {10}}-1} $

再次进行分析扩展：

$$ | H（s）| ^ 2 = \ frac {1} {1 + \ epsilon ^ 2 T_N ^ 2 \ left（\ frac {s} {j \ omega_c} \ right）} $$

，当分母为零时，$ s = p_n $也是极点。

$$ 1 + \ epsilon ^ 2 T_N ^ 2 \ left（\ frac {p_n} { j \ omega_c} \ right）= 0 $$

或

$$ T_N \ left（\ frac {p_n} {j \ omega_c} \ right）= \ pm \ frac {j} {\ epsilon} $$

（因为$ \ cos（\ theta）= \ cosh（j \ theta）$，我们可以使用$ \ cos（）$或$ \ $ T_N（）$的cosh（）$表达式$

$$ \ cosh \ big（N \，\ operatorname {arccosh} \ left（\ frac {p_n} {j \ omega_c} \ right） \ big）= \ pm \ frac {j} {\ epsilon} $$

$$ N \，\ operatorname {arccosh} \ left（\ frac {p_n} {j \ omega_c} \ right ）= \ operatorname {arccosh} \ left（\ pm \ frac {j} {\ epsilon} \ right）$$

因为$$ y = \ cosh（x）= \ tfrac12（e ^ x + e ^ {-x}）$$
和$$ x = \ operatorname {arccosh}（y）= \ log \ left（y \ pm \ sqrt {y ^ 2-1} \ right）$ $

t hen

$$ N \ log \ left（\ frac {p_n} {j \ omega_c} \ pm \ sqrt {\ left（\ frac {p_n} {j \ omega_c} \ right）^ 2 -1} \ right）= \ log \ left（\ pm \ frac {j} {\ epsilon} \ pm \ sqrt {\ left（\ pm \ frac {j} {\ epsilon} \ right）^ 2-1} \ right）$$

$$ N \ log \ left（\ frac {\ Re（p_n）+ j \ Im（p_n）} {j \ omega_c} \ pm \ sqrt {\ left（ \ frac {\ Re（p_n）+ j \ Im（p_n）} {j \ omega_c} \ right）^ 2-1-1} \ right）= \ log \ left（\ pm j \ left（\ frac {1} { \ epsilon} \ pm \ sqrt {\ frac {1} {\ epsilon ^ 2} +1} \ right）\ right）$$

$$ N \ log \ left（\ frac {- j \ Re（p_n）+ \ Im（p_n）} {\ omega_c} \ pm \ sqrt {\ left（\ frac {-j \ Re（p_n）+ \ Im（p_n）} {\ omega_c} \ right）^ 2-1} \ right）\\ = \ log \ left（\ pm j \ left（\ frac {1} {\ epsilon} \ pm \ sqrt {\ frac {1} {\ epsilon ^ 2} +1} \对）\ right）$$

哦，亲爱的我可能在12小时之内就不会爆这个东西了。

我认为我太懒了，无法解决这个问题。如果有人想拿起它，请随时。在复数值的矩形和极坐标表示法之间进行lota转换。记得什么时候

$$ w = \ pm \ sqrt {\ z \} $$
然后
$$ | w | = + \ sqrt {| z |} $$
和
$$ \ begin {align} \ arg \ {w \}＆= \ frac12 \ arg \ {z \} + \ arg \ { \ pm 1 \} \\
＆= \ frac12 \ arg \ {z \} + \ frac {\ pi} {2}（1 \ pm 1）\ end {align} $$

并记住
$$ \ log（z）= \ log | z | + j \ arg \ {z \} + j 2 \ pi n \ quad n \ in \ mathbb {Z} $$

您可以加上$ 2 \ pi $的任何整数倍（例如“ $ 2 \ pi n $“）到任何$ \ arg \ {\ cdot \} $（选择右侧的$ \ log（）$，这就是为$ p_n $获取不同极点的方式）。

如果您喜欢带有复杂变量的数学自慰，请把自己踢出去。

#4 楼

看来，在此讨论中的大多数参与者都不知道哪种过滤器可能是他们真正的解决方案！即Paynter滤波器是由麻省理工学院教授，​​Philbrick Reseach的合伙人Henry M.Paynter开发的。
它们是“运行”平均滤波和处理非确定性输入信号的最佳方法，效果更好。比贝塞尔-汤姆森我将它们用于生理，医学和声纳应用。 “设计有源低通滤波器”，作者Peter D. Hansen
给出了二阶，四阶和六阶滤波器的极点表。我对8阶计算了相同的结果。

#5 楼

12阶椭圆到4阶椭圆

（我没有资格获得赏金。）我试图在Octave中提出问题3的反例，但令人惊讶的是我没有。如果对问题3的回答为是，则根据问题的图5，应在阶数为4的椭圆滤波器和阶数为12的椭圆滤波器之间共享特定的极点和零点。 图1.顺序为$ N = 12 $和$ N = 4 $的椭圆滤波器之间的极点和零点，蓝色，并按参数曲线$ f（\ alpha）$的递增参数$ \ alpha $编号。

让我们设计一个具有一些任意参数的12阶椭圆滤波器：1 dB通带纹波，-90 dB阻带纹波，截止频率0.1234，s平面而不是z平面：
图2.使用ellip设计的12阶椭圆滤波器的幅频响应。

pkg load signal;
[b12, a12] = ellip (12, 0.1, 90, 0.1234, "s");
ra12 = roots(a12);
rb12 = roots(b12);
freqs(b12, a12, [0:10000]/10000);

图3.使用ellip设计的12阶椭圆滤波器的极点（红色）和零点（蓝色）。水平轴：实部，垂直轴：虚部。

让我们通过重用12阶滤波器的选择极点和零点来构造4阶滤波器，如图1所示。虚部的极点和零点就足够了：

scatter(vertcat(real(ra12), real(rb12)), vertcat(imag(ra12), imag(rb12)));

图4.所有极点和零点与某些零点相同的四阶滤波器的幅频响应按照图1所示，是12阶滤波器的放大倍数。放大后可以看到滤波器的特征：3.14 dB通带等效波纹，-27.69 dB阻带等效波纹，截止频率0.1234。据我了解，具有等于极点和零点数量所允许的波纹的等波纹通带和等波纹阻带是说滤波器是椭圆形的充分条件。但是，让我们尝试通过用图3所示的特性设计ellip设计4阶椭圆滤波器并比较两个4阶滤波器之间的零点和零点来确认这一点：
反对：

[~, ira12] = sort(imag(ra12));
[~, irb12] = sort(imag(rb12));
ra4 = [ra12(ira12)(2), ra12(ira12)(5), ra12(ira12)(8), ra12(ira12)(11)];
rb4 = [rb12(irb12)(2), rb12(irb12)(5), rb12(irb12)(8), rb12(irb12)(11)];
freqs(poly(rb4), poly(ra4), [0:10000]/10000);

图5. ellip设计的4阶滤波器（交叉点）之间极点（红色）和零（蓝色）位置的比较）和一个与12阶滤波器共享某些极点和零位置的4阶滤波器（圆形）。水平轴：实部，垂直轴：虚部。

两个滤波器之间的极点和零点重合到小数点后三位，这是表征12阶滤波器的精度。结论是，至少在这种特定情况下，通过均匀地将它们分布在相同的参数曲线上，至少可以获得四阶椭圆滤波器和十二阶椭圆滤波器的极点和零点。 。滤波器不是Butterworth或Chebyshev I或II型滤波器，因为通带和阻带都具有波纹。

4阶椭圆到12阶椭圆

相反，可以从4阶ellip滤波器的极点和零点拟合的一对连续函数可以近似得出12阶滤波器的极点和零点？

如果我们复制四个极点（图5）并翻转复制品实部的符号，我们会得到各种椭圆形。当我们绕椭圆旋转时，经过的极点位置会给出一个周期性的离散序列。通过对其离散傅立叶变换（DFT）进行零填充，它是周期性带限内插的理想选择。在产生的$ 24 $极中，那些正数为正的极将被丢弃，将极数减半至$ 12 $。插值的倒数不是零，而是插值，否则插值的完成方式与极点相同。我们从与前面相同的ellip设计的4阶滤波器开始（大约与图4相同）：

[b4el, a4el] = ellip (4, 3.14, 27.69, 0.1234, "s");
rb4el = roots(b4el);
ra4el = roots(a4el);
scatter(vertcat(real(ra4), real(rb4)), vertcat(imag(ra4), imag(rb4)));

图6. 12阶的幅值频率响应滤波器的极点和零点从与4阶滤波器的曲线相匹配的曲线中采样。

图2响应的模型不够准确，无法发挥作用。阻带的效果很好，但通带却倾斜了。频带边缘频率近似正确。尽管如此，这仍然表明考虑参数曲线的潜力，每个曲线仅用4个自由度来描述。

让我们看看极点和零点如何与$ N = 12 $ ellip-生成的滤波器的极点和零点匹配：

scatter(vertcat(real(ra4el), real(rb4el)), vertcat(imag(ra4el), imag(rb4el)), "blue", "x");

图7. ellip设计的12阶滤波器（十字）和从4阶滤波器派生的12阶滤波器（圆圈）之间的极点（红色）和零（蓝色）位置的比较。水平轴：实部，垂直轴：虚部。

内插极点相差很大，但零点的匹配相对较好。应该研究较大的$ N $作为起点。

6阶椭圆到18阶椭圆

进行与上述相同的操作，但从6阶开始并内插到18阶，显示出看似良好的幅度频率响应，但仔细检查时通带仍存在问题：

pkg load signal;
[b4el, a4el] = ellip (4, 3.14, 27.69, 0.1234, "s");
rb4el = roots(b4el);
ra4el = roots(a4el);
rb4eli = 1./rb4el;
[~, ira4el] = sort(imag(ra4el));
[~, irb4eli] = sort(imag(rb4eli));
ra4eld = vertcat(ra4el(ira4el), -ra4el(ira4el));
rb4elid = vertcat(rb4eli(irb4eli), -rb4eli(irb4eli));
ra12syn = -interpft(ra4eld, 24)(12:23);
rb12syn = -1./interpft(rb4elid, 24)(12:23);
freqs(poly(rb12syn), poly(ra12syn), [0:10000]/10000);

图8.顶部）6阶ellip生成的滤波器，底部）从6阶滤波器派生的18阶滤波器。放大后，通带只有两个最大值，纹波约为1 dB。阻带几乎等于2.5 dB的变化。

我对通带的麻烦的猜测是，带限制的插值在（的实际部分）不能很好地工作

椭圆滤波器的精确曲线

结果证明，对于$ NN_z = NN_p = N $的椭圆滤波器，它们为问题1和2提供了正例。C. Sidney Burrus，数字信号处理和数字滤波器设计（草稿）。 OpenStax CNX。 2012年11月18日给出了一个足够通用的$ NN_z = NN_p = N $椭圆滤波器的传递函数的零点和极点，以Jacobi椭圆正弦$ \ operatorname {sn}（t，k）。$表示运算符{sn}（t，k）=-\运算符{sn}（-t，k），$ Burrus等式。 3.136可以用零$ s_ {zi}，$ $ i = 1 \ dots N $重写为：


$$ s_ {zi} = \ frac {j} {k \运算符名称{sn} \ big（K + K（2i + 1）/ N，k \ big）}，\ tag {1} $$

其中$ K $是$ \的四分之一周期实际$ t $的operatorname {sn}（t，k）$和$ 0 \ le k \ le 1 $可以看作是滤波器参数化的自由度。它控制相对于通带宽度的过渡带宽度。识别$（2i + 1）/ N = 2 \ alpha $（请参阅问题的等式2）其中$ \ alpha $是参数曲线的参数：

$$ f_z（\ alpha ）= \ frac {j} {k \ operatorname {sn}（K + 2K \ alpha，k）} \ tag {2}，$$

Burrus等式。 3.146给出了左上角的四分之一平面极点，包括一个奇数$ N $的实极点。可以为所有极点$ s_ {pi}，$ $ i = 1 \ dots N $以及任何$ N $重写为：

$$ s_ {pi} =
\ frac {\ operatorname {cn} \ big（K + K（2i + 1）/ N，k \ big）\ operatorname {dn} \ big（K + K（2i + 1）/ N，k \ big）\运算符名称{sn} \ big（\ nu_0，\ sqrt {1-k ^ 2} \ big）\\
\ times \ operatorname {cn} \ big（\ nu_0，\ sqrt {1-k ^ 2} \ big）
+ j \ operatorname {sn} \ big（K + K（2i + 1）/ N，k \ big）\ operatorname {dn} \ big（\ nu_0，\ sqrt {1-k ^ 2} \ big）}
{1- \ operatorname {dn} ^ 2 \ big（K + K（2i + 1）/ N，k \ big）\ operatorname {sn} ^ 2 \ big（\ nu_0 ，\ sqrt {1-k ^ 2} \ big）}，\ tag {3} $$

其中$ \ operatorname {dn}（t，k）= \ sqrt {1-k ^ 2 \ operatorname {sn} ^ 2（t，k）} $是Jacobi椭圆函数之一。一些源将$ k ^ 2 $作为所有这些函数的第二个参数，并将其称为模数。我们有$ k $并将其称为模数。变量$ 0 <\ nu_0
通过与零相同的替代，极的参数曲线可以写为：

$$ f_p（\ alpha） =
\ frac {\ operatorname {cn} \ big（K + 2K \ alpha，k \ big）\ operatorname {dn} \ big（K + 2K \ alpha，k \ big）\ operatorname {sn} \ big（\ nu_0，\ sqrt {1-k ^ 2} \ big）\\
\ times \ operatorname {cn} \ big（\ nu_0，\ sqrt {1-k ^ 2} \ big）
+ j \ operatorname {sn} \ big（K + 2K \ alpha，k \ big）\ operatorname {dn} \ big（\ nu_0，\ sqrt {1-k ^ 2} \ big）}
{1 -\ operatorname {dn} ^ 2 \ big（K + 2K \ alpha，k \ big）\ operatorname {sn} ^ 2 \ big（\ nu_0，\ sqrt {1-k ^ 2} \ big）}。 \ tag {4} $$

让我们绘制从Burrus示例3.4复制的$ k $和$ \ nu_0 $（代码中的v0）值的Octave函数和曲线。 />

[b12, a12] = ellip (12, 0.1, 90, 0.1234, "s");
ra12 = roots(a12);
rb12 = roots(b12);
scatter(vertcat(real(ra12), real(rb12)), vertcat(imag(ra12), imag(rb12)), "blue", "x");
scatter(vertcat(real(ra12syn), real(rb12syn)), vertcat(imag(ra12syn), imag(rb12syn)));

图9. Burrus示例3.4的$ f_z（\ alpha）$和$ f_p（\ alpha）$，经分析扩展到期间$ \ alpha = 0 \ dots2 $。本示例中的三个极点（红叉）和三个零点（蓝色圆圈，一个无穷大且未显示）是相对于$ \ alpha $在$ \ alpha = 1/6，$ $ \ alpha = 3 /根据等式，来自这些函数的6，$和$ \ alpha = 5/6，$。问题2。有了扩展名，$ \ operatorname {Im} \ big（f_z（\ alpha）\ big）$（未显示）的倒数会非常缓慢地振荡，从而可以像前面部分中那样通过截断的傅立叶级数来近似。其他周期性扩展函数也很平滑，但是用这种方式近似起来不太容易。

[b6el, a6el] = ellip (6, 0.03, 30, 0.1234, "s");
rb6el = roots(b6el);
ra6el = roots(a6el);
rb6eli = 1./rb6el;
[~, ira6el] = sort(imag(ra6el));
[~, irb6eli] = sort(imag(rb6eli));
ra6eld = vertcat(ra6el(ira6el), -ra6el(ira6el));
rb6elid = vertcat(rb6eli(irb6eli), -rb6eli(irb6eli));
ra18syn = -interpft(ra6eld, 36)(18:35);
rb18syn = -1./interpft(rb6elid, 36)(18:35);
freqs(poly(rb18syn), poly(ra18syn), [0:10000]/10000);

图10. Burrus示例3.4的参数曲线。水平轴：实部，垂直轴：虚部。该视图未显示参数曲线的速度，因此三个极点（红叉）和三个零点（蓝色圆圈，一个无穷大且未显示）似乎没有均匀分布在曲线上，即使它们是相对于参数曲线的参数$ \ alpha $。

由Burrus给出的精确极点和零公式进行的椭圆滤波器设计完全等效于从精确的$ f_p（\ alpha）$和$ f_z（\ alpha）$中采样，因此方法是等效的并且可用。问题1仍是不限成员名额。可能其他类型的过滤器具有由$ f_p（\ alpha）$和$ f_z（\ alpha）$和$ N_z / N_p $定义的无限子集。在近似椭圆参数曲线的方法中，那些不依赖于确切函数形式的方法可能会转移到其他滤波器类型，我认为最有可能推广到一般椭圆滤波器的方法，例如一般等波纹滤波器的一些子集。对于它们，极点和零点的确切公式可能是未知的或难以理解的。

回到等式。 2，对于奇数$ N $，我们有零个值$ \ alpha = 0.5 $，它通过$ \ operatorname {sn} \ big（2K，k \ big）= 0 $将其发送到无穷大。两极（方程式4）不会发生这种情况。我已经更新了问题，在计数$ NN_z $（或$ NN_p $）中包含这样的零（和极点，以防万一）。在$ k = 0 $时，所有零都根据$ f_z（\ alpha）$变为无穷大，这看起来像是类型I的切比雪夫滤波器。

我认为问题3刚刚解决，答案是“是”。如此看来，我们可以涵盖所有椭圆滤波器的情况，而不会与$ NN_z = NN_p $冲突，使用这些新定义。

评论

---

$ \ begingroup $
Olli，无论如何您都不能给自己赏金。您的500点数将永远消失。只是不要浪费它们，就像我在EE.SE页面上一次不小心那样。
$ \ endgroup $
–罗伯特·布里斯托-约翰逊
16年11月1日在20:56

---

$ \ begingroup $
评论不用于扩展讨论；此对话已移至聊天。
$ \ endgroup $
–jojek♦
16年2月2日，9：48

---

$ \ begingroup $
是的，它们仍然存在，这是奇数阶的特殊情况，这时传递函数还有一个额外的单实数极点$ re /（s + re）$。对于仅从零开始的有理函数，您具有：$$ R（x）= \ frac {\ prod_i ^ {n / 2} {x ^ 2-zero_i ^ 2}} {\ prod_j ^ {n / 2 } {x ^ 2-k / zero_j ^ 2}} $$，其中$ k = fs / fp $。对于奇数阶，$ R（x）= R（x）* x $。这将使滤波器具有非归一化的增益，因此应按$ R（0）$进行缩放。极点来自展开$ 1 + \ epsilon ^ 2 R ^ 2（x）$并找到分母的根，然后选择左侧，并形成传递函数。
$ \ endgroup $
–有关的公民
16年3月3日18:00

---

$ \ begingroup $
我不确定我是否这样说。要建立传递函数，并不一定要遵循这本书，先制作$ H（s）H（-s）$，然后是左侧极点，然后根据@A_A的公式进行有理传递函数。从数学上和实际的结果来看，是从$ 1 + \ epsilon ^ 2R ^ 2（x）$找到根后（注意：$ x $，而不是$ j \ omega $或$ s $），只需选择根具有正实部和正或负imagpart（不是两者）。即对于$ N = 4 $，将有4对/ 8极；选择后，您有2个不同的极点。然后简单地说：$$ N（s）= \ prod_i ^ {N / 2} {| p_i | ^ 2} $$ ...
$ \ endgroup $
–有关的公民
16年11月4日在7:26

---

$ \ begingroup $
（用于全极点滤波器），其中$ p = \ sigma + j \ omega $，而$$ N（s）= \ prod_i ^ {N / 2} {s ^ 2 + | z_i | ^ 2 } $$，其中$ z = j \ mu $（用于零极点滤波器），而分母：$$ D（s）= \ prod_j ^ {N / 2} {s ^ 2 + 2 * Re（p_j） * s + | p_j | ^ 2} $$和$$ H（s）= \ frac {N（s）} {D（s）} $$。这将是低通原型。
$ \ endgroup $
–有关的公民
16年11月4日在7:31

---

#6 楼

我将在此处添加一些注释，如果有人要计算传递函数的第N个$ N $量级的极限$ N \ to \ infty $，且传递函数具有多个$ N $极点和零点，则可能会有用。任意参数曲线。通过使用较大的$ N $并将极点和零点均匀地分布在参数曲线的参数上，可以近似得出该值。不幸的是，在可实现传递函数的极点和零点位置处，近似值始终具有dB级的无限误差。从这个意义上讲，更好的构件是沿其长度具有均匀的极点或零分布的线段。仅考虑$ N $个零，它们分布在起点$ x_0 + y_0i $和终点$ x_1 + y_1i的线段上：$

$$ \ lim_ {N \ to \ infty} | H （0）| ^ {1 / N} = \ prod_0 ^ 1 \ Bigg | \ left（x_0 + y_0i \ right）（1- \ alpha）+ \ left（x_1 + y_1i \ right）\ alpha \ Bigg | ^ { d \ alpha} \\
= \ prod_0 ^ 1 \ left（\ sqrt {\ left（x_0（1-\ alpha）+ x_1x \ right）^ 2 + \ left（y_0（1-\ alpha）+ y_1 \ alpha \ right）^ 2} \ right）^ {d \ alpha} \\
= e ^ {\ displaystyle \ int_0 ^ 1 \ log \ left（\ sqrt {\ left（x_0（1-\ alpha）+ x_1 \ alpha \ right）^ 2 + \ left（y_0（1-\ alpha）+ y_1 \ alpha \ right）^ 2} \ right）d \ alpha} \\
= e ^ {\ left（\ displaystyle \ frac {\ left（x_0y_1-x_1y_0 \ right）\ operatorname {atan2} \ left（x_0y_1-x_1y_0，x_0x_1 + y_0y_1 \ right）} {\ left（x_0-x_1 \ right）^ 2 + \ left （y_0-y_1 \ right）^ 2}-1 \ right）} \\
\ times \ left（x_0 ^ 2 + y_0 ^ 2 \ right）^ {\ left（\ displaystyle \ frac {x_1 \ left （x_0-x_1 \ right）+ y_0 \ left（y_0-y_1 \ right）+ \ left（x_0-x_1 \ right）^ 2} {2 \ left（\ left（x_0-x_1 \ right）^ 2 + \ left （y_0-y_1 \ right）^ 2 \ right）} \ right）} \\
\ times \ left（x_1 ^ 2 + y_1 ^ 2 \ right）^ {\ left（\ displaystyle \ frac {x_0 \ left（x_1-x_0 \ right）+ y_1 \ left（y_1-y_0 \ right）+ \ left（x_1-x_0 \ right）^ 2} {2 \ left（\ left （x_1-x_0 \ right）^ 2 + \ left（y_1-y_0 \ right）^ 2 \ right）} \ right）} $$

某些特殊情况需要单独处理。如果$ x_0 = 0 $和$ y_0 = 0 $，则必须使用限制：
$$ = e ^ {-1} \ sqrt {x_1 ^ 2 + y_1 ^ 2} $$

或者相反，如果$ x_1 = 0 $和$ y_1 = 0 $：
$$ = e ^ {-1} \ sqrt {x_0 ^ 2 + y_0 ^ 2} $$

或如果线段的长度为零，$ x_0 = x_1 $并且$ y_0 = y_1 $，则我们只有一个常规零：
$$ = \ sqrt {x_0 ^ 2 + y_0 ^ 2} $$

要对$ H（z）$或$ H（s）$的不同自变量值进行求值，只需从行的起点和终点减去该值即可。

这是什么在复平面上看起来像：
图1.传递函数的大小，具有单个零。绿松石色表示1 dB步长，黄色表示10 dB步长。

图2.带$的传递函数的第N个$ N $量级的极限$ N \ to \ infty $ N $零均匀分布在线段上。线段处有折痕，但该值永远不会像规则的，可实现的零那样变为零。在足够的距离处，这看起来像是常规零。颜色代码与图1中的相同。

图3.使用离散零的近似图2：具有5个零的多项式的均方根的第5个根均匀地分布在直线上分割。在每个零的位置，该值为零，因为$ 0 ^ {1/5} =0。$

1和2是使用此处理草图生成的，其源代码为：

float[] dragPoints;
int dragPoint;
float dragPointBackup0, dragPointBackup1;
boolean dragging, activated;
PFont fnt;
PImage bg;
float pi = 2*acos(0.0);
int appW, appH;
float originX, originY, scale;

int numDragPoints = 2;

void setup() {
  appW = 600;
  appH = 400;
  originX = appW/2;
  originY = appH/2;
  scale = appH*7/16;
  size(600, 400);
  bg = createImage(appW, appH, RGB);
  dragging = false;
  dragPoint = -666;
  dragPoints = new float[numDragPoints*2]; 
  dragPoints[0] = originX-appW*0.125;
  dragPoints[1] = originY+appH*0.125;
  dragPoints[2] = originX+appW*0.125;
  dragPoints[3] = originY-appH*0.125;
  fnt = createFont("Arial",16,true);
  ellipseMode(RADIUS);
  activated = false;
}

void findDragPoint() {
  int cutoff = 49;
  int oldDragPoint = dragPoint;
  float dragPointD = 666666666;
  dragPoint = -666;
  for (int t = 0; t < numDragPoints; t++) {
    float d2 = (mouseX-dragPoints[t*2])*(mouseX-dragPoints[t*2]) + (mouseY-dragPoints[t*2+1])*(mouseY-dragPoints[t*2+1]);
    if (d2 <= dragPointD) {
       dragPointD = d2;
       if (dragPointD < cutoff) {
         dragPoint = t;
       }
    }
  }
  if (dragPoint != oldDragPoint) {
    loop();
  }
}

void mouseMoved() {
  if (activated) {
    if (!dragging) {
      findDragPoint();
      loop();
    }
  }
}

void mouseClicked() {
  if (dragPoint < 0) {
    activated = !activated;
    if (activated) {
      findDragPoint();      
    }
  }
  loop();
}

void mousePressed() {  
  if (dragPoint >= 0) {
    dragging = true;
    dragPointBackup0 = dragPoints[dragPoint*2];
    dragPointBackup1 = dragPoints[dragPoint*2+1];
  } else {
    dragging = false; // Not needed?
  }
  loop();
}

void mouseDragged() {
  if (!activated) {
    dragPoint = -666;
    activated = true;
    findDragPoint();
  }
  if (dragging) {
    int x = mouseX;
    int y = mouseY;
    if (x < 5) {
      x = 5;
    } else if (x >= appW - 5) {
      x = appW - 6;
    }
    if (y < 5) {
      y = 5;
    } else if (y >= appH - 5) {
      y = appH - 6;
    }
    dragPoints[dragPoint*2] = x;
    dragPoints[dragPoint*2+1] = y;
    loop();
  }  
}

void mouseReleased() {
  if (activated && dragging) {
    dragging = false;
    loop();
  }
}

float sign(float value) {
  if (value > 0) {
    return 1.0;
  } else if (value < 0) {
    return -1.0;
  } else {
    return 0;
  }
}

void draw() {
  for(int y = 0; y < appH; y++) {
    for(int x = 0; x < appW; x++) {
      float x0 = (dragPoints[0]-x)/scale;
      float y0 = (dragPoints[1]-y)/scale;
      float x1 = (dragPoints[2]-x)/scale;
      float y1 = (dragPoints[3]-y)/scale;
      float gain;
      if (x0 == x1 && y0 == y1) {
        gain = sqrt(x0*x0 + y0*y0);
      } else if (x0 == 0 && y0 == 0) {
        gain = exp(-1)*sqrt(x1*x1 + y1*y1);
      } else if (x1 == 0 && y1 == 0) {
        gain = exp(-1)*sqrt(x0*x0 + y0*y0);
      } else {
        gain = exp((x0*y1 - x1*y0)*atan2(x0*y1 - x1*y0, x0*x1 + y0*y1)/(sq(x0 - x1) + sq(y0 - y1)) - 1)*pow(x0*x0 + y0*y0, (x1*(x0 - x1) + y0*(y0 - y1) + sq(x0 - x1))/(2*(sq(x0 - x1) + sq(y0 - y1))))*pow(x1*x1 + y1*y1, (x0*(x1 - x0) + y1*(y1 - y0) + sq(x1 - x0))/(2*(sq(x1 - x0) + sq(y1 - y0))));
      }
      int intensity10 = round(log(gain)/log(10)*0x200)&0xff;
      int intensity1 = round(log(gain)/log(10)*(0x200*10))&0xff;
      bg.pixels[y*appW + x] = color(intensity10, 0xff, intensity1);
    }
  }
  image(bg, 0, 0);
  noFill();
  stroke(0, 0, 255);
  strokeWeight(1);
  line(dragPoints[0], dragPoints[1], dragPoints[2], dragPoints[3]);  

  //ellipse(originX, originY, scale, scale);  
  if (!activated) {
    textFont(fnt,16);
    fill(0, 0, 0);
    text("Click to activate",10,20);
    for (int x = 0; x < appW; x++) {
      color c = color(110*x/appW+128, 110*x/appW+128, 110*x/appW+128);
      set(x, 0, c);  
    }
    for (int y = 0; y < appH; y++) {
      color c = color(110*y/appH+128, 110*y/appH+128, 110*y/appH+128);
      set(0, y, c);  
    }
  }

  for (int u = 0; u < numDragPoints; u++) {
    stroke(0, 0, 255);
    if (dragPoint == u) {
      if (dragging) {
        fill(0, 0, 255);
        strokeWeight(3);
        ellipse(dragPoints[u*2], dragPoints[u*2+1], 5, 5);
      } else {
        noFill();
        strokeWeight(3);
        ellipse(dragPoints[u*2], dragPoints[u*2+1], 6, 6);
      }
    } else {
      //noFill();
      //strokeWeight(1);
      //ellipse(dragPoints[u*2], dragPoints[u*2+1], 6, 6);
    }
  }
  noLoop();
}