我已阅读主题:移动平均线的截止频率
我在我的算法中使用第二个答案来计算我的滤波器的 3dB 截止频率,效果很好,因为我的滤波器长度通常在 300 以上。我用阶跃响应验证了它。
但我想知道这个公式的来源或推导。
我在第二个和第三个学期后手动尝试了泰勒级数。我接近但不完全符合公式,mapple 给了我一个有效但极其复杂的结果。
希望大家能帮忙。
亲切的问候
Slev1n
我已阅读主题:移动平均线的截止频率
我在我的算法中使用第二个答案来计算我的滤波器的 3dB 截止频率,效果很好,因为我的滤波器长度通常在 300 以上。我用阶跃响应验证了它。
但我想知道这个公式的来源或推导。
我在第二个和第三个学期后手动尝试了泰勒级数。我接近但不完全符合公式,mapple 给了我一个有效但极其复杂的结果。
希望大家能帮忙。
亲切的问候
Slev1n
考虑长度为 $N$ 的零相位移动平均线:
$$\text{y}[n] = \begin{cases} \displaystyle\frac{\text{x}[n] + \displaystyle\sum_{k=1}^{\frac{N-1}{2 }}\left(\text{x}[n+k] + \text{x}[nk]\right)}{N},&n\in\mathbb{Z}&\text{for }N\text{奇数}\\ \displaystyle\frac{\displaystyle\sum_{k=1}^{\frac{N}{2}}\left(\text{x}[n+(k-\frac{1}{2} )] + \text{x}[n-(k-\frac{1}{2})]\right)}{N},&n+\frac{1}{2}\in\mathbb{Z}&\ text{for }N\text{ even} \end{cases}$$
在具有整数时间索引的离散序列上操作的偶数长度滤波器不能是零相位。我们通过使输出时间索引始终具有 $\frac{1}{2}$ 的小数部分来规避这一点,即使是 $N$。作为一个真实世界的例子,如果输入在每个午夜进行采样,则将在每个中午计算偶数长度的零相位移动平均值。这种不寻常的索引方便地给出了相同的零相位形式的频率响应 $F_N(\omega)$ 用于 $N$ 奇数和 $N$ 偶数:
$$F_N(\omega) = \begin{cases} \displaystyle\frac{\displaystyle1 + \displaystyle\sum_{k=1}^{\frac{N-1}{2}}\left(e^{ik \omega}+e^{-ik\omega}\right)}{N}&\text{对于 }N\text{ 奇数}\\ \displaystyle\frac{\displaystyle\sum_{k=1}^{\ frac{N}{2}}\left(e^{i(k-\frac{1}{2})\omega}+e^{-i(k-\frac{1}{2})\omega }\right)}{N}&\text{for }N\text{ even} \end{cases}\\ = \frac{\sin(\frac{N\omega}{2})}{N\sin (\frac{\omega}{2})}$$
不幸的是,频率响应对于 -3 dB 截止频率 $\omega_c$ 没有符号解,例如:
$$F_N(\omega_c)=\sqrt{\frac{1}{2}}.$$
严格来说 $\sqrt{\frac{1}{2}}$ 约为 -3.01 dB,但我认为这就是人们所说的 -3 dB 的意思,否则它只是一个任意数字。一个近似的频率响应 $\hat{F}_N(\omega)$ 使用积分而不是总和:
$$\hat{F}_N(\omega)=\frac{1}{N}\int_{-\frac{N}{2}}^\frac{N}{2}e^{ik\omega} = \frac{2\sin(\frac{N\omega}{2})}{N\omega}$$
真实(总和)和近似(积分)频率响应的主瓣在大 $N$ 处收敛:
我们可以通过引入函数 $G_N(\chi) = F_N(\omega)$ 和 $\hat{G}_N(\chi) = \hat{F}_N(\omega)$ 来证明收敛性,参数标准化为$\omega = \frac{2\pi\chi}{N}$,将两个函数的第一个零带到 $\chi = 1$:
$$G_N(\chi) = \frac{\sin(\pi\chi)}{N\sin\left(\frac{\pi\chi}{N}\right)}\\ \hat{G}_N (\chi) = \frac{\sin(\pi\chi)}{\pi\chi}\\ \lim_{N\rightarrow \infty}G_N(\chi) = \frac{\sin(\pi\chi )}{\pi\chi}$$
$G_N(\chi)$ 被称为$N$-周期性带限脉冲序列。它的极限 $N$ 和函数 $\hat{G}_N(\chi)$ 都是$\text{sinc}$ 函数。不幸的是,-3 dB 截止频率在近似 $\hat{F}_N(\omega)$ 中也没有符号解。对于不同的$N$,近似值与$N = 1$ 的近似值仅通过映射$\omega \rightarrow \omega N$ 不同,因此求解近似的-3 dB 截止频率$\hat\omega_c(N )$ 数字表示 $N = 1$:
$$\frac{2\sin(\frac{\hat\omega_c(1)}{2})}{\hat\omega_c(1)} = \sqrt{\frac{1}{2}}\\ \ Rightarrow\hat\omega_c(1) = 2.78311475650302030063992,$$
给出任意 $N$ 的近似截止频率:
$$\hat\omega_c(N) = \frac{\hat\omega_c(1)}{N} $$
这似乎是另一个比 Massimo 更简单的近似。对于您的 $N > 300$,使用它应该没有问题。Massimo 和这个答案的常数与以下因素相关:
$$\frac{\hat\omega_c(1)}{2\pi} = 0.442946470689452340308369.$$
我进一步看了一下,发现 Massimo 用 $\hat{F}_M(\omega)$ 逼近 $F_N(\omega)$,选择 $M$ 使得频率响应的(的)二阶导数和$\omega = 0$ 处的近似匹配:
$$F''_N(\omega)=\frac{\sin(\frac{N \omega}{2}) \left(2 \sin(\omega) \cos(\frac{\omega}{2} ) + (N^2 - 1) \sin(\frac{\omega}{2}) (\cos(\omega) - 1)\right)}{2 N (\cos(\omega) - 1)^ 2} - \frac{2 N \sin(\omega) \sin(\omega/2) \cos(\frac{N \omega}{2})}{2 N (\cos(\omega) - 1) ^2}\\ \hat{F}''_M(\omega)=\left(\frac{4}{M \omega^3} - \frac{M}{2 \omega}\right) \sin\左(\frac{M \omega}{2}\right) - \frac{2 \cos(\frac{M \omega}{2})}{\omega^2}\\ \lim_{\omega\rightarrow 0}F''_N(\omega) = \frac{1-N^2}{12}\\\lim_{\omega\rightarrow 0}\hat{F}''_M(\omega) = -\frac {M}{12}\\ \lim_{\omega\rightarrow 0}F''_N(\omega) = \lim_{\omega\rightarrow 0}\hat{F}''_M(\omega)\\ \右箭头 M = \sqrt{N^2 - 1}\\ \右箭头 \hat\omega_c(M) = \frac{\hat\omega_c(1)}{\sqrt{N^2 - 1}} $$
这改进了小 $\omega$ 处的近似值,其中包括 -3dB 截止点,尤其是在小 $N$ 处:
Massimo 的近似值总是高估截止频率(参见误差比较),通过改变常数 $1$ 来改进它。$N = 2$ 时误差最大。如果它的误差被限制为等于 $N = 3$ 处的(当前第二大)误差,我们会得到一个更好但同样便宜的近似值:
$$\hat\omega_c(N) = \frac{2.78311475650302030063992}{\sqrt{N^2 - 0.8603235992780290790596}}$$
常量的这个和其他调整,比如Matt 的常量 $0.863031932778066$,对于大的 $N$ 工作得非常好(参见错误比较)。对于较大的 $N$,$N$ 每增加 10 倍,误差就会下降 1000 倍。对这些事情的解释是,作为 $N$ 函数的真实截止频率具有Laurent 级数:
$$\omega_c(N) = \displaystyle\sum_{k=0}^{\infty}\frac{a_k}{N^k},\ a_k = 0\text{ if }k\text{ even},$ $
近似值及其 Laurent 级数为:
$$\hat\omega_c(N) = \frac{a}{\sqrt{N^2+c}}\\= \frac{a}{N} - \frac{ac}{2N^3} + \ frac{3ac^2}{8N^5} - \frac{5ac^3}{16N^7} + \frac{35ac^4}{128N^9} - \ldots,$$
这样: $$a_1 = a = 2.78311475650302030063992\\ a_3 \approx -\frac{ac}{2}$$
如果 $N^{-3}$-term 中的近似匹配是精确的,那么近似误差应该减少 $10^5$,因为大的 $N$ 增加了 $10$。函数 $f(x)$ 的 Laurent 级数 $f(x) = \sum_{k=0}^{\infty}\frac{a_k}{x^k}$ 的系数 $a_k$ as $x \rightarrow\infty$ 可以通过以下方式迭代找到:
$$f_0(x) = f(x)\\ a_k = \lim_{x\rightarrow\infty}f_k(x)\\ f_{k+1}(x) = \left(f_{k}(x) -a_k\right)x$$
当我们没有符号形式的 $f(x)$,但可以对非常大的 $x$ 以任意精度数值求解时,我们可以在数值上执行上述过程的等价。以下 Python 脚本使用SymPy
并将mpmath
计算给定数量(此处为 10)的第一个系数 $a_k$,以达到真实截止频率的 Laurent 系列所需的精度:
from sympy import *
from mpmath import *
num_terms = 10
num_decimals = 24
num_use_decimals = num_decimals + 5 #Ad hoc headroom
def y(omega):
return sin(N*omega/2)/(N*sin(omega/2))-sqrt(0.5)
a = []
h = mpf('1e'+str(num_decimals))
for i in range(0, num_terms):
mp.dps = 2*2**(num_terms - i)*num_use_decimals*(i + 1) #Ad hoc headroom
N = mpf('1e'+str(2**(num_terms - i)*num_use_decimals))
r = findroot(y, [2.7/N, pi/N]) #Safe search range
for j in range(0, i):
r = (r - a[j])*N
a.append(r)
mp.dps = num_decimals
print 'a_'+str(i)+' = '+str(nint(r*h)/h)
在我的电脑上,程序运行了大约 7 分钟。它打印以下内容,表明 Laurent 级数仅包含奇数的负幂:
a_0 = 0.0
a_1 = 2.78311475650302030063992
a_2 = 0.0
a_3 = 1.20103120331802052770586
a_4 = 0.0
a_5 = 0.767601805674195789947052
a_6 = 0.0
a_7 = 0.537174869947196574599614
a_8 = 0.0
a_9 = 0.387986305863698148870773
这些数字显示到小数点后 24 位,并不是洛朗级数唯一意义上的近似值。没有其他 Laurent 级数等于 $\omega_c(N)$。仅使用 $a_1$ 和 $a_3$,可以构造一个简单的两项截断 Laurent 级数近似:
$$\hat\omega_c(N) = \frac{2.78311475650302030063992}{N} + \frac{1.20103120331802052770586}{N^3},$$
并通过 $c=-\frac{2a_3}{a_1}$ 的近似值:
$$\hat\omega_c(N) = \frac{2.78311475650302030063992}{\sqrt{N^2-0.863084212040982824494534}}$$
两者在大 $N$ 时都有 $1/N^5$ 误差衰减,分别参见误差比较,列 h) 和 i)。截断较长的 Laurent 级数,具有来自脚本输出的更多项,衰减更快,在误差比较中,对于列 j) 处的 5 项近似值,$1/N^{11}$。
让我们比较截止频率的不同近似值的实际数值误差。表中给出的误差是通过从近似值中减去真实的(数值求解的)-3 dB 截止频率 $\omega_c$ 计算得出的。
$$\begin{array}{rl} \text{a})&\frac{2.78311475650302030063992}{N}\\ \text{b})&\frac{2.78311475650302030063992}{\sqrt{N^2-1}} \\ \text{c})&\frac{2.78311475650302030063992}{\sqrt{N^2 - 0.8603235992780290790596}}\\ \text{d})&\sqrt{\frac{12}{2 N^2 - 1} }\\ \text{e})&\sqrt{ \frac{15(2N^2-1) - 6\sqrt{5 \left(N^4 - 5N^2 - \tfrac{13}{4} \右)}}{2N^4 - 1} }\\ \text{f})&\frac{2.78311475650302030063992}{\sqrt{N^2 - 0.863031932778066}}\\ \text{g})&\sum_{i =1}^{5}\frac{a_{2i-1}}{N^{2i-1}}\\\text{h})&\frac{2.78311475650302030063992}{N} + \frac{1.20103120331802052770586}{ N^3}\\ \text{i})&\frac{2.78311475650302030063992}{\sqrt{N^2-0.863084212040982824494534}}\\ \text{j})&\sum_{i=1}^{5}\ frac{a_{2i-1}}{N^{2i-1}}\\ \end{array}$$ 近似 g) 的系数 $a_{2i-1}$ 可以在这个答案中找到作者:Matt L。近似 j) 的系数 $a_{2i-1}$ 来自 $\omega_c(N)$ 的 Laurent 系列,可以在这个答案中找到由奥利。$$\small\begin{array}{l} \begin{array}{r|r|l} N&\omega_c&\text{a)}&\text{b)}&\text{c)}&\text {d)}&\text{e)}\\\hline 2 & 1.57079632679491894158970 & \text{-1.8E-01} & \text{3.6E-02} & \text{-1.1E-04 } & \text {-2.6E-01} & \text{n/a}\\ 3 & 0.97561347715844261315892 & \text{-4.8E-02} & \text{8.4E-03} & \text{-1.1E-04 } & \text{-1.4E-01} & \text{7.7E-02}\\ 4 & 0.71532874990708873792278 & \text{-2.0E-02} & \text{3.3E-03} & \text{-5.4E- 05 } & \text{-9.3E-02} & \text{3.6E-02}\\ 5 & 0.56648391394608301620215 & \text{-9.9E-03} & \text{1.6E-03} & \text{- 2.9E-05 } & \text{-7.2E-02} & \text{2.4E-02}\\ 6 & 0.46951346150003474555592 & \text{-5.7E-03} & \text{9.2E-04} & \文本{-1.7E-05 } & \文本{-5.8E-02} &
注:近似 e) 不允许 $N=2$。一些误差被列为 0,但这仅意味着它们的大小小于 1E-17 左右。这和其他可能的不准确性是由于在计算近似值和误差时使用了双精度浮点运算。
随意编辑/添加另一个近似值。
来自我的向上箭头,奥利。
但出于某种原因,我认为答案要简单得多。通常我喜欢设计非因果对称 FIR 滤波器,因为它们是零相位,但通常我将自己限制为奇数个非零抽头。为了更普遍地做到这一点,我可能会坚持因果 FIR 移动平均线。
假设抽头数为 $N$。
$$ y[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} x[nk] $$
应用 $\mathcal{Z}$-transform(和几何求和公式):
$$ \begin{align} Y(z) & = \frac{1}{N} \sum\limits_{k=0}^{N-1} X(z) z^{-k} \\ & = X(z) \frac{1}{N} \sum\limits_{k=0}^{N-1} z^{-k} \\ & = X(z) \frac{1}{N} \ frac{1 - z^{-N}}{1 - z^{-1}} \\ \end{align} $$
代入 $z \ \leftarrow \ e^{j \omega}$ 得到 DTFT:
$$ \begin{align} Y(e^{j \omega}) & = X(e^{j \omega}) \ \frac{1}{N} \frac{1 - (e^{j \omega })^{-N}}{1 - (e^{j \omega})^{-1}} \\ & = X(e^{j \omega}) \ \frac{1}{N} \ frac{1 - e^{-j \omega N}}{1 - e^{-j \omega}} \\ & = X(e^{j \omega}) \ \frac{e^{-j \欧米茄 N/2}}{N \ e^{-j \omega/2}} \ \frac{e^{j \omega N/2} - e^{-j \omega N/2}}{e^ {j \omega/2} - e^{-j \omega/2}} \\ & = X(e^{j \omega}) \ \frac{e^{-j \omega (N-1)/ 2}}{N} \ \frac{e^{j \omega N/2} - e^{-j \omega N/2}}{e^{j \omega/2} - e^{-j \ omega/2}} \\ & = X(e^{j \omega}) \ \frac{e^{-j \omega (N-1)/2}}{N} \ \frac{ \sin(\欧米茄 N/2)}{ \sin(\omega/2) } \\ \end{align} $$
通常我们将乘以 $X(z)$ 的东西称为“传递函数”
$$ H(z) = \frac{1}{N} \frac{1 - z^{-N}}{1 - z^{-1}} $$
以及乘以 $X(e^{j \omega})$ 的东西,即“频率响应”
$$ H(e^{j \omega}) = e^{-j \omega (N-1)/2} \ \frac{\sin(\omega N/2)}{N \ \sin(\omega /2)} $$
$e^{-j \omega (N-1)/2}$ 因子表示 $\frac{N-1}{2}$ 个样本的线性相位、恒定延迟。它不会改变增益。
$\frac{\sin(\omega N/2)}{N \ \sin(\omega/2)}$ 因子是增益因子。“-3 dB 频率”,$\omega_c$,(通常我们指的是 -3.0103 dB 频率,因为它对应于“半功率”频率)是这样的
$$ \左| H(e^{j \omega_c}) \right|^2 = \frac{1}{2} $$
要么
$$ \left( \frac{\sin(\omega_c N/2)}{N \ \sin(\omega_c/2)} \right)^2 = \frac{1}{2} $$
要么
$$ 2 \ \sin^2(\omega_c N/2) = N^2 \ \sin^2(\omega_c/2) \ .$$
所以给定抽头数$N$,你必须求解$\omega_c$。对于封闭式表格,这可能并不容易做到,但是您可以挖出计算器并插上插头,直到得到足够精确的答案。或者你可以让 MATLAB 来做。
通过使用三角恒等式(我在摆弄双线性变换时常用的一个)和$\cos( )$。
$$ \begin{align} \sin^2(\theta) & = \frac{1}{2}(1 \ - \ \cos(2 \theta)) \\ & \approx \frac{1}{2 }\left(1 \ - \ \left(1 - \frac{(2\theta)^2}{2!} + \frac{(2\theta)^4}{4!} \right) \right) \\ & = \frac{1}{2}\left( \frac{(2\theta)^2}{2!} - \frac{(2\theta)^4}{4!} \right) \ \ & = \theta^2 \left(1 - \frac{\theta^2}{3} \right) \\ \end{align} $$
如果您在前面的等式中插入 $\sin^2()$ 的近似值并求解...(跳过很多步骤,因为我懒得 $\LaTeX$ 出来...)
你得到
$$ \omega_c \approx \sqrt{\frac{12}{2 N^2 - 1}} \xrightarrow{ \quad N \to \infty \quad} \frac{\sqrt{6}}{N} $$
奥利,这与你的结果相比如何?
用另一个近似 $\sin^2()$ 的术语来做这个更好,是可行的,只需要 $\omega_0^2$ 的二次解。使用的近似值(保留 $\cos()$ 扩展的前四项)是:
$$ \begin{align} \sin^2(\theta) & = \frac{1}{2}(1 \ - \ \cos(2 \theta)) \\ & \approx \frac{1}{2 }\left(1 \ - \ \left(1 - \frac{(2\theta)^2}{2!} + \frac{(2\theta)^4}{4!} - \frac{(2 \theta)^6}{6!} \right) \right) \\ & = \frac{1}{2}\left( \frac{(2\theta)^2}{2!} - \frac{ (2\theta)^4}{4!} + \frac{(2\theta)^6}{6!} \right) \\ & = \theta^2 \left(1 - \frac{1}{ 3}\theta^2 + \frac{2}{45}\theta^4 \right) \\ \end{align} $$
尝试该近似值并求解 $\omega_c^2$。
我得到的最一致的答案是
$$ \begin{align} \omega_c^2 & = (15)\frac{2N^2-1}{2N^4-1} \pm \sqrt{ \left((15)\frac{2N^2- 1}{2N^4-1} \right)^2 - \frac{360}{2N^4-1} } \\ & = \frac{15(2N^2-1) \pm 6\sqrt{5 \left(N^4 - 5N^2 - \tfrac{13}{4} \right)}}{2N^4 - 1} \\ \\ \omega_c^2 & \xrightarrow{ \quad N \to \infty \quad} \frac{15 \pm 3\sqrt{5}}{N^2} \end{align} $$
使用 "$+$" 选项,它看起来像
$$ \omega_c \xrightarrow{ \quad N \to \infty \quad} \frac{\sqrt{21.7}}{N} $$
并且使用“$-$”选项看起来像
$$ \omega_c \xrightarrow{ \quad N \to \infty \quad} \frac{\sqrt{8.3}}{N} $$
这更接近我上面所做的“一阶”近似。所以我想我会选择“$-$”选项。
所以,即使我不能分析性地说明为什么“$+$”选项应该被拒绝,我想我最准确的答案是
$$ \omega_c = \sqrt{ \frac{15(2N^2-1) - 6\sqrt{5 \left(N^4 - 5N^2 - \tfrac{13}{4} \right)}}{ 2N^4 - 1} } $$
对于大的 $N$,它有限制,如上所示。
还有其他人有更好的方法来查看一个好的近似封闭形式的解决方案吗?
好的,这很有趣。我将添加自己的想法和近似值,其中第一个结果与 Massimo 在此答案中给出的相同,而 Olli 在此线程中得出的相同。我仍然将它包括在这里,因为它的推导不同。然后我将展示一个更好的近似值,对于 $N=2$,它的最大相对误差为 $0.002$(对于这种情况,我们当然有一个精确截止频率的解析解:$\omega_c=\pi/ 2$),并且对于 $N\ge 10$,相对误差小于 $1.2\cdot 10^{-4}$。
众所周知,Olli 和 Robert 在他们的回答中表明,长度为 $N$ 的移动平均滤波器的实值幅度函数由下式给出
$$H_N(\omega)=\frac{\sin\left(\frac{N\omega}{2}\right)}{N\sin\left(\frac{\omega}{2}\right)} \标签{1}$$
$3$ dB 截止频率 $\omega_c$ 满足
$$H_N(\omega_c)=\frac{1}{\sqrt{2}}\tag{2}$$
据我所知,方程没有相当简单的解析解。$(2)$。这里提出的近似值的关键是 - 毫不奇怪 - 泰勒近似值。罗伯特答案中使用的泰勒级数之间的区别在于,我没有单独逼近正弦函数(或罗伯特答案中的平方值),而是直接逼近 $(1)$ 中给出的完整幅度函数。逼近 $\sin(N\omega/2)$(或其平方值)将导致比逼近完整函数时更大的误差,因为参数 $N\omega/2$ 永远不会接近零,即使对于较大的值$N$。只逼近分母 $\sin(\omega/2)$(或其平方值)是可以的,因为它的参数 $\omega=\omega_c$ 对于大的 $N$ 确实接近于零。无论如何,我不会使用这两个近似值,但我将使用 $H_N(\omega)$ 的泰勒级数。对于更简单的表示法,我将使用 $x=\omega/2$ 和 $F_N(x)=H_N(\omega)$。$x_0=0$ 附近的 $F_N(x)$ 的泰勒级数由下式给出
$$F_N(x)\约 1-\frac{N^2-1}{6}x^2+\frac{3N^4-10N^2+7}{360}x^4\tag{3} $$
对于较大的 $N$ 值,这种近似是合理的,因为截止频率 $\omega_c$ 倾向于较小的值。
对于第一个近似值,我只使用 $(3)$ 中的前两项:
$$1-\frac{N^2-1}{6}x_c^2=\frac{1}{\sqrt{2}}\tag{4}$$
求解 $(4)$ 给出第一个近似解:
$$\omega_{c1}=2x_c=\frac{2\sqrt{6(1-\frac{1}{\sqrt{2}})}}{\sqrt{N^2-1}}=\frac {2.65130859228473}{\sqrt{N^2-1}}\tag{5}$$
这个解决方案的问题在于它是有偏差的,这意味着它的误差对于大的 $N$ 不会收敛到零。然而,事实证明,通过 $(5)$ 的简单缩放,可以消除这种偏差。对于零偏差,我们需要
$$\lim_{N\rightarrow \infty}H_N(\omega_{c1}(N))=\frac{1}{\sqrt{2}}\tag{6}$$
我使用符号 $\omega_{c1}(N)$ 来强调它对 $N$ 的依赖。用一般表达式求解 $(6)$
$$\omega_{c1}=\frac{a}{\sqrt{N^2-1}}\tag{7}$$
将我们引向等式
$$\frac{\sin(a/2)}{a/2}=\frac{1}{\sqrt{2}}\tag{8}$$
对于现在著名的解决方案,必须对其进行数值求解
$$a=2.78311475650302\标签{9}$$
由 $(9)$ 给出的 $(7)$ 与 $a$ 的近似值与 Massimo 的公式相同(您必须除以 $2\pi$ 才能得到他的魔常数),它也与由奥利在这个线程中以不同的方式。我们看到泰勒近似为我们提供了方程的正确形式,但常数必须通过极限过程来确定,才能得到零偏差的公式。对于大多数实际目的,此公式足够准确,对于 $N\ge 10$ 的最大相对误差为 $6.9\cdot 10^{-4}$。
$$\omega_{c2}(N)=2\sqrt{6}\sqrt{\frac{5(N^2-1)-\sqrt{5}\sqrt{(3\sqrt{2}-1) N^4+10(1-\sqrt{2})N^2+7\sqrt{2}-9}}{3N^4-10N^2+7}}\tag{10}$$
请注意,在四次方程的四个解中,我们需要选择两个正解中较小的一个,因为这是泰勒级数非常接近 $F_N(x)$ 的值。另一个积极的解决方案是泰勒级数偏离 $F_N(x)$ 的范围内的伪影。近似值 $(10)$ 与 $(5)$ 给出的先前近似值的第一个版本有相同的小问题,因为它有一个小的偏差。可以通过考虑限制 $(6)$ 以与以前完全相同的方式消除这种偏差,这次是 $\omega_{c2}(N)$。我基于 $(10)$ 但零偏差的最终近似值由下式给出
$$\omega_{c3}(N)=b\omega_{c2}(N)\tag{11}$$
其中 $b$ 也可以通过求解类似于 $(8)$ 的方程得到。它实际上可以写成由 $(9)$ 给出的 $a$:
$$b=\frac{a}{2\sqrt{2}\sqrt{5-\sqrt{5(3\sqrt{2}-1)}}}=0.997314251642175\tag{12}$$
$$\epsilon_i=\frac{\omega_c-\omega_{ci}}{\omega_c}\tag{13}$$
这允许比较不同的近似值$\omega_{ci}$。我将只讨论零偏差的近似值:$\omega_{c1}$ 由 $(7)$ 给出,$a$ 由 $(9)$ 给出,$\omega_{c3}$ 由 $(11 )$(和 $(10)$),其中 $b$ 由 $(12)$ 给出。下图显示了由 $(13)$ 定义为 $N$ 的函数的相对误差。红色曲线是近似$(7)$的相对误差,绿色曲线是近似$(11)$的误差。两种近似都具有零偏差(它们收敛到大 $N$ 的精确值),但绿色曲线的收敛速度明显更快。
上面显示的零偏置公式是实际截止频率的近似值,但更好的公式(公式 $(10,11,12)$)非常尴尬。奥利有一个好主意,可以在简单的公式 $(7)$ 中调整分母常数。只要我们使用 $(9)$ 给出的 $a$ 的最优值,我们就可以改变分母常数而不会失去零偏差属性。所以我们得到了一个新的公式
$$\omega_{c4}(N)=\frac{a}{\sqrt{N^2-c}}\tag{14}$$
用常数 $c$ 进行优化。如果我理解正确的话,Olli 对 $c$ 的优化是基于 $N=2$ 的错误值。但是,我认为 $N=2$ 的值不是很相关,因为对于 $N=2$,$\omega_c$ 可以解析计算:$\omega_c(2)=\pi/2$。所以我们不一定需要为 $N=2$ 的情况优化公式 $(14)$,如果它是以牺牲 $N$ 较大值的近似为代价的。我通过以下方式优化了 $(14)$ 中的常量 $c$。如果 $\omega_c(N)$ 是给定一组滤波器长度 $N$ 的精确截止频率,我们有一个超定方程组:
$$\omega_c(N)\stackrel{!}{=}\frac{a}{\sqrt{N^2-c}},\qquad N=2,3,\ldots\tag{15}$$
我们可以为 $N$ 选择任何合理的值集。重新排列 $(15)$ 给出另一组方程,这一次在未知 $c$ 中是线性的:
$$N^2-c\stackrel{!}{=}\frac{a^2}{\omega_c^2(N)}\tag{16}$$
$(16)$ 的最优最小二乘解是
$$c_0=\frac{1}{L}\sum_N\left(N^2-\frac{a^2}{\omega_c(N)}\right)\tag{17}$$
其中 $L$ 是总和中使用的 $N$ 的不同值的数量。如果您使用 $[2,100]$ 范围内的 $N$ 的所有整数值,您将得到
$$c_0=0.863031932778066\标签{18}$$
这接近 Olli 的值,但它为所有 $N\ge 3$ 提供了更好的近似值。我将错误值添加到该表的 f) 列。
在他的回答中,罗伯特想知道为什么在对 $\sin^2(x)$ 使用四阶泰勒级数时,他必须放弃 $\omega_c$ 的第二个(更大的)正解。下图说明了原因。原始平方幅度函数显示为蓝色($N=10$)。3dB 线为红色。绿色函数是泰勒近似,两次穿过红线。这是$\omega_c$ 的两个正解。由于函数是偶数,我们也得到相同的两个带负号的解,这使它成为四个,就像四阶多项式的情况一样。然而,很明显,两个正解中较大的一个是人为的,因为较大参数的泰勒近似存在差异。因此,只有较小的解决方案才有意义,而另一个则没有。