漫谈FIR(上)
前言
在广大的工程师看来,FIR是一个滤波器工具,对它的设计、仿真、分析,matlab等软件提供了堪称完美的工具。对于工程师而言,基本上已经足够了,既然如此,还有什么可以说的呢?可以说的只是个人的、零散的、浅薄的知识碎碎念罢了。
本文自推公式较多,也不是学数字信号处理的,如有错漏,还望包涵。
变换
工程上,力学的、电子的系统(当然还有其他的,比如种群数量的增长等)往往采用线性微分方程去抽象,对方程的求解,我能掌握的方法有三种:待定系数法、微分算子法和Laplace变换法。
在远去的大学生涯中,求解线性电路的响应的时候,我就喜欢用Laplace变换的方法,一个微分方程,经过变换,成了简单的代数方程。特别的,当初始条件为零,s域中电容的阻抗是1sC\frac{1}{\text{sC}}sC1,电感的阻抗是sL\text{sL}sL,复杂的微分方程可以用阻抗的串并联来等效。即使simulink等工具可以完美的仿真,但是我还是喜欢拿起A4纸,写一写、画一画。这一爱好在往后的工作中得以延续。
相较于Laplace变换,Fourier变换也是工作中的常客。
周期为N的数列的Fourier级数对,常见的写法如下:
x[n]=∑k=1Nakejkω0n=∑k=1Nakejk2πNn(1.1)x[n]=\sum\limits_{\text{k}=1}^{\text{N}}{{{a}_{k}}{{e}^{jk{{\omega }_{0}}n}}}=\sum\limits_{\text{k}=1}^{\text{N}}{{{a}_{k}}{{e}^{jk\frac{2\pi }{N}n}}}\tag{1.1}x[n]=k=1∑Nakejkω0n=k=1∑NakejkN2πn(1.1)
ak=1N∑n=1Nx[n]e−jkω0n=1N∑n=1Nx[n]e−jk2πNn(1.2){{a}_{k}}=\frac{1}{N}\sum\limits_{\text{n}=1}^{\text{N}}{x[n]{{e}^{-jk{{\omega }_{0}}n}}}=\frac{1}{N}\sum\limits_{\text{n}=1}^{\text{N}}{x[n]{{e}^{-jk\frac{2\pi }{\text{N}}n}}}\tag{1.2}ak=N1n=1∑Nx[n]e−jkω0n=N1n=1∑Nx[n]e−jkN2πn(1.2)
式(1.1)表示将x[n]x[n]x[n]分解为复指数ejkω0n{{e}^{jk{{\omega }_{0}}n}}ejkω0n的线性组合,式(1.2)是一个内积。但是由于式(1.2)当中突兀的 1N\frac{1}{N}N1,显得两公式并不对称,同样的,函数的级数展开也有一个12π\frac{1}{2\pi }2π1的问题:
x(t)=∑k=−∞+∞akejkt(1.3)x(t)=\sum\limits_{k=-\infty }^{+\infty }{{{a}_{k}}{{e}^{jkt}}}\tag{1.3}x(t)=k=−∞∑+∞akejkt(1.3)
ak=12π∫−ππx(t)e−jktdt(1.4){{a}_{k}}=\frac{1}{2\pi }\int_{-\pi }^{\pi }{x(t){{e}^{-jkt}}}dt\tag{1.4}ak=2π1∫−ππx(t)e−jktdt(1.4)
相信这是令广大工程师困惑的地方,数字信号处理的书籍,也读过几本,对于这点都没有过解释,类似的困惑相信还有很多。或许由于这些作者都是工科领域的,公式仅仅是工具,抛出来了,记住了就行了,洋洋洒洒的笔墨更多的落在应用上面。
在数学家们写的分析学作品中,Fourier级数是L2[ab]{{\text{L}}^{2}}[\begin{matrix}
a & b \\\end{matrix}]L2[ab]空间上的一个普通的正交分解,一个普通的有界线性算子,尽管着墨并不多,却能让我们对它的理解远比工科类书籍的叙述来得深刻得多。有的书籍(比如Rudin的《数学分析原理》、Albert Boggess的《小波与傅里叶分析基础》)会提到一个标准正交系:
{...,cos(2x)π,cos(x)π,12π,sin(x)π,sin(2x)π,...}(1.5)\left\{ ...,\frac{\cos (2x)}{\sqrt{\pi }} \right.,\frac{\cos (x)}{\sqrt{\pi }},\frac{1}{\sqrt{2\pi }},\frac{\sin (x)}{\sqrt{\pi }},\frac{\sin (2x)}{\sqrt{\pi }},\left. ... \right\}\tag{1.5}{...,πcos(2x),πcos(x),2π1,πsin(x),πsin(2x),...}(1.5)
12π\frac{1}{2\pi }2π1的存在与否便取决于正交系是否是标准的。式(1.5)经过Euler公式变换之后带入式(1.4),12π\frac{1}{2\pi }2π1便消失了。序列的级数,经过类似的操作,1N\frac{1}{N}N1也会消失。
N维复空间CN{{\mathbb{C}}^{\text{N}}}CN可以由向量
Vk=[ejkω01,ejkω02,...,ejkω0(N−1),ejkω0N]T(k=1,2,3,...,N)(1.6){{\mathbf{V}}_{\text{k}}}\text{=}{{\left[ {{e}^{jk{{\omega }_{0}}1}},{{e}^{jk{{\omega }_{0}}2}},...,{{e}^{jk{{\omega }_{0}}(N-1)}},{{e}^{jk{{\omega }_{0}}N}} \right]}^{\text{T}}}(k=1,2,3,...,N)\tag{1.6}Vk=[ejkω01,ejkω02,...,ejkω0(N−1),ejkω0N]T(k=1,2,3,...,N)(1.6)
张成,并且Vk{{\mathbf{V}}_{\text{k}}}Vk和Vr{{\mathbf{V}}_{\text{r}}}Vr的内积
⟨Vk,Vr⟩=VkTVr∗=∑n=1Nejkω0ne−jrω0n=∑n=1Nej(k−r)ω0n={Nk=r0k≠r(1.7)\left\langle {{\mathbf{V}}_{\text{k}}},{{\mathbf{V}}_{r}} \right\rangle \text{=}\mathbf{V}_{k}^{T}\mathbf{V}_{\text{r}}^{*}=\sum\limits_{n=1}^{N}{{{e}^{jk{{\omega }_{0}}\text{n}}}{{e}^{-jr{{\omega }_{0}}\text{n}}}}=\sum\limits_{n=1}^{N}{{{e}^{j(k-r){{\omega }_{0}}\text{n}}}}=\left\{ \begin{matrix}
N\begin{matrix}
{} & k=r \\
\end{matrix} \\
0\begin{matrix}
{} & k\ne r \\
\end{matrix} \\
\end{matrix} \right.\tag{1.7}
⟨Vk,Vr⟩=VkTVr∗=n=1∑Nejkω0ne−jrω0n=n=1∑Nej(k−r)ω0n={Nk=r0k=r(1.7)
这组基是正交的,但不是标准的,我们进一步对它进行标准化:对Vk{{\mathbf{V}}_{\text{k}}}Vk乘以系数1N\frac{1}{\sqrt{\text{N}}}N1。
V^k=1NVk=[1Nejkω01,1Nejkω02,...,1Nejkω0(N−1),1Nejkω0N]T(1.8){{\mathbf{\hat{V}}}_{\text{k}}}=\frac{1}{\sqrt{\text{N}}}{{\mathbf{V}}_{\text{k}}}={{\left[ \frac{1}{\sqrt{\text{N}}}{{e}^{jk{{\omega }_{0}}1}},\frac{1}{\sqrt{\text{N}}}{{e}^{jk{{\omega }_{0}}2}},...,\frac{1}{\sqrt{\text{N}}}{{e}^{jk{{\omega }_{0}}(N-1)}},\frac{1}{\sqrt{\text{N}}}{{e}^{jk{{\omega }_{0}}N}} \right]}^{\text{T}}}\tag{1.8}V^k=N1Vk=[N1ejkω01,N1ejkω02,...,N1ejkω0(N−1),N1ejkω0N]T(1.8)
即
⟨V^k,V^r⟩={1k=r0k≠r(1.9)\left\langle {{{\mathbf{\hat{V}}}}_{\text{k}}},{{{\mathbf{\hat{V}}}}_{r}} \right\rangle \text{=}\left\{ \begin{matrix}
1\begin{matrix}
{} & k=r \\
\end{matrix} \\
0\begin{matrix}
{} & k\ne r \\
\end{matrix} \\
\end{matrix} \right.\tag{1.9}
⟨V^k,V^r⟩={1k=r0k=r(1.9)
接下来我们看一下X=[x1,x2,...,xN−1,xN]T\mathbf{X}={{\left[ {{x}_{1}},{{x}_{2}},...,{{x}_{N-1}},{{x}_{N}} \right]}^{\text{T}}}X=[x1,x2,...,xN−1,xN]T在每一个V^k{{\mathbf{\hat{V}}}_{\text{k}}}V^k上的投影:
qk=⟨X,V^k⟩=∑n=1Nxn1Ne−jkω0n(1.10){{\text{q}}_{\text{k}}}=\left\langle \mathbf{X},{{{\mathbf{\hat{V}}}}_{\text{k}}} \right\rangle \text{=}\sum\limits_{n=1}^{N}{{{x}_{n}}\frac{1}{\sqrt{\text{N}}}{{e}^{-jk{{\omega }_{0}}n}}}\tag{1.10}qk=⟨X,V^k⟩=n=1∑NxnN1e−jkω0n(1.10)
在这里,我们把“1Ne−jkω0n\frac{1}{\sqrt{\text{N}}}{{e}^{-jk{{\omega }_{0}}n}}N1e−jkω0n”看做一个整体,于是式(1.2)当中的1N\frac{1}{N}N1便消失了。而且可得
qk=Nak(1.11){{\text{q}}_{\text{k}}}=\sqrt{\text{N}}{{a}_{k}}\tag{1.11}qk=Nak(1.11)
qk{{\text{q}}_{\text{k}}}qk组合成的向量
q=[q1,q2,...,qN−1,qN]T\mathbf{q}={{[{{\text{q}}_{\text{1}}},{{\text{q}}_{2}},...,{{\text{q}}_{N-1}},{{\text{q}}_{N}}]}^{\text{T}}}q=[q1,q2,...,qN−1,qN]T
是一个N维的向量。V^k{{\mathbf{\hat{V}}}_{\text{k}}}V^k组合成的矩阵
A=[V^1,V^2,...,V^N-1,V^N]\mathbf{A}=\left[ {{{\mathbf{\hat{V}}}}_{\text{1}}},{{{\mathbf{\hat{V}}}}_{2}},...,{{{\mathbf{\hat{V}}}}_{\text{N-1}}},{{{\mathbf{\hat{V}}}}_{N}} \right]A=[V^1,V^2,...,V^N-1,V^N]
是一个N阶方阵。不难看出
AHX=q(1.12){{\mathbf{A}}^{\text{H}}}\mathbf{X}=\mathbf{q}\tag{1.12}AHX=q(1.12)
式中AH{{\mathbf{A}}^{\text{H}}}AH是A\mathbf{A}A的Hermitian转置。因此我喜欢把q\mathbf{q}q看成X\mathbf{X}X经过AH{{\mathbf{A}}^{\text{H}}}AH的线性变换后的输出
X→AH q\mathbf{X}\overset{{{\mathbf{A}}^{\text{H}}}}{\mathop{\to }}\,\mathbf{q}X→AHq
显然的,AH{{\mathbf{A}}^{\text{H}}}AH是一个映射:X×X→X\mathbf{X}\times \mathbf{X}\to \mathbf{X}X×X→X。我们将式(1.12)展开
[1Ne−j12πN11Ne−j12πN2...1Ne−j12πN(N−1)1N1Ne−j22πN11Ne−j22πN2...1Ne−j22πN(N−1)1N...............1Ne−j(N−1)2πN11Ne−j(N−1)2πN2...1Ne−j(N−1)2πN(N−1)1N1N1N...1N1N][x1x2...xN−1xN]=[q1q2...qN−1qN](1.13)\left[ \begin{matrix}
\frac{1}{\sqrt{\text{N}}}{{e}^{-j1\frac{2\pi }{\text{N}}1}} & \frac{1}{\sqrt{\text{N}}}{{e}^{-j1\frac{2\pi }{\text{N}}2}} & ... & \frac{1}{\sqrt{\text{N}}}{{e}^{-j1\frac{2\pi }{\text{N}}(N-1)}} & \frac{1}{\sqrt{\text{N}}} \\
\frac{1}{\sqrt{\text{N}}}{{e}^{-j2\frac{2\pi }{\text{N}}1}} & \frac{1}{\sqrt{\text{N}}}{{e}^{-j2\frac{2\pi }{\text{N}}2}} & ... & \frac{1}{\sqrt{\text{N}}}{{e}^{-j2\frac{2\pi }{\text{N}}(N-1)}} & \frac{1}{\sqrt{\text{N}}} \\
... & ... & ... & ... & ... \\
\frac{1}{\sqrt{\text{N}}}{{e}^{-j(N-1)\frac{2\pi }{\text{N}}1}} & \frac{1}{\sqrt{\text{N}}}{{e}^{-j(N-1)\frac{2\pi }{\text{N}}2}} & ... & \frac{1}{\sqrt{\text{N}}}{{e}^{-j(N-1)\frac{2\pi }{\text{N}}(N-1)}} & \frac{1}{\sqrt{\text{N}}} \\
\frac{1}{\sqrt{\text{N}}} & \frac{1}{\sqrt{\text{N}}} & ... & \frac{1}{\sqrt{\text{N}}} & \frac{1}{\sqrt{\text{N}}} \\
\end{matrix} \right]\left[ \begin{matrix}
{{x}_{1}} \\
{{x}_{2}} \\
... \\
{{x}_{N-1}} \\
{{x}_{N}} \\
\end{matrix} \right]=\left[ \begin{matrix}
{{q}_{1}} \\
{{q}_{2}} \\
... \\
{{q}_{N-1}} \\
{{q}_{N}} \\
\end{matrix} \right]
\tag{1.13}N1e−j1N2π1N1e−j2N2π1...N1e−j(N−1)N2π1N1N1e−j1N2π2N1e−j2N2π2...N1e−j(N−1)N2π2N1...............N1e−j1N2π(N−1)N1e−j2N2π(N−1)...N1e−j(N−1)N2π(N−1)N1N1N1...N1N1x1x2...xN−1xN=q1q2...qN−1qN(1.13)
将变换用矩阵写出,它的线性性质便是一目了然的了。熟悉线性代数的工程师们可以看出,(AH)−1=A=AT{{({{\mathbf{A}}^{\text{H}}})}^{-1}}=\mathbf{A}={{\mathbf{A}}^{\text{T}}}(AH)−1=A=AT,是unitary矩阵(数学家翻译成酉矩阵,物理学家翻译成幺正矩阵,真是搞事情),它包含了本节的全部内容,我们要做的只是对它进行阐述。对式(1.13)左乘A\mathbf{A}A可得
[x1x2...xN−1xN]=[1Nej12πN11Nej12πN2...1Nej12πN(N−1)1N1Nej22πN11Nej22πN2...1Nej22πN(N−1)1N...............1Nej(N−1)2πN11Nej(N−1)2πN2...1Nej(N−1)2πN(N−1)1N1N1N...1N1N][q1q2...qN−1qN](1.14)\left[ \begin{matrix}
{{x}_{1}} \\
{{x}_{2}} \\
... \\
{{x}_{N-1}} \\
{{x}_{N}} \\
\end{matrix} \right]=\left[ \begin{matrix}
\frac{1}{\sqrt{\text{N}}}{{e}^{j1\frac{2\pi }{\text{N}}1}} & \frac{1}{\sqrt{\text{N}}}{{e}^{j1\frac{2\pi }{\text{N}}2}} & ... & \frac{1}{\sqrt{\text{N}}}{{e}^{j1\frac{2\pi }{\text{N}}(N-1)}} & \frac{1}{\sqrt{\text{N}}} \\
\frac{1}{\sqrt{\text{N}}}{{e}^{j2\frac{2\pi }{\text{N}}1}} & \frac{1}{\sqrt{\text{N}}}{{e}^{j2\frac{2\pi }{\text{N}}2}} & ... & \frac{1}{\sqrt{\text{N}}}{{e}^{j2\frac{2\pi }{\text{N}}(N-1)}} & \frac{1}{\sqrt{\text{N}}} \\
... & ... & ... & ... & ... \\
\frac{1}{\sqrt{\text{N}}}{{e}^{j(N-1)\frac{2\pi }{\text{N}}1}} & \frac{1}{\sqrt{\text{N}}}{{e}^{j(N-1)\frac{2\pi }{\text{N}}2}} & ... & \frac{1}{\sqrt{\text{N}}}{{e}^{j(N-1)\frac{2\pi }{\text{N}}(N-1)}} & \frac{1}{\sqrt{\text{N}}} \\
\frac{1}{\sqrt{\text{N}}} & \frac{1}{\sqrt{\text{N}}} & ... & \frac{1}{\sqrt{\text{N}}} & \frac{1}{\sqrt{\text{N}}} \\
\end{matrix} \right]\left[ \begin{matrix}
{{q}_{1}} \\
{{q}_{2}} \\
... \\
{{q}_{N-1}} \\
{{q}_{N}} \\
\end{matrix} \right]
\tag{1.14}x1x2...xN−1xN=N1ej1N2π1N1ej2N2π1...N1ej(N−1)N2π1N1N1ej1N2π2N1ej2N2π2...N1ej(N−1)N2π2N1...............N1ej1N2π(N−1)N1ej2N2π(N−1)...N1ej(N−1)N2π(N−1)N1N1N1...N1N1q1q2...qN−1qN(1.14)
即
X=Aq(1.15)\mathbf{X}=\mathbf{Aq}\tag{1.15}X=Aq(1.15)
即
xn=∑k=1Nqk1Nejkω0n(1.16){{x}_{n}}=\sum\limits_{\text{k}=1}^{\text{N}}{{{q}_{k}}\frac{1}{\sqrt{N}}{{e}^{jk{{\omega }_{0}}n}}}\tag{1.16}xn=k=1∑NqkN1ejkω0n(1.16)
X=∑k=1NqkV^k(1.17)\mathbf{X}=\sum\limits_{k=1}^{N}{{{\text{q}}_{\text{k}}}{{{\mathbf{\hat{V}}}}_{k}}}\tag{1.17}X=k=1∑NqkV^k(1.17)
式(1.10)和式(1.16)组成一个Fourier级数对。它和式(1.1)和式(1.2)组成的变换对是等价的,我更喜欢前者,式(1.10)表示X\mathbf{X}X在各个正交的单位向量上的投影,式(1.17)表示X\mathbf{X}X是V^k{{\mathbf{\hat{V}}}_{k}}V^k.的线性组合,式(1.16)是式(1.17)的再次分解——这一系列的变换对,都体现出了一个对称性,这点非常好(对称性与守恒律总是形影相伴,在下面的叙述中,相信读者会慢慢体会到)。在一些量子力学著作(比如Cohen的),也有类似的推导过程,只不过没展开写,并且用的是狄拉克符号,没这么啰嗦。
事实上,根据线性代数的正交分解的性质,式(1.17)可以由式(1.10)直接导出,还可以直接导出
∥X∥2=∑k=1N∣qk∣2(1.18){{\left\| \mathbf{X} \right\|}^{2}}=\sum\limits_{\text{k}=1}^{N}{{{\left| {{\text{q}}_{\text{k}}} \right|}^{2}}}\tag{1.18}∥X∥2=k=1∑N∣qk∣2(1.18)
∑n=1N∣xn∣2=∑k=1N∣qk∣2(1.19)\sum\limits_{\text{n}=1}^{N}{{{\left| {{x}_{n}} \right|}^{2}}}=\sum\limits_{\text{k}=1}^{N}{{{\left| {{\text{q}}_{\text{k}}} \right|}^{2}}}\tag{1.19}n=1∑N∣xn∣2=k=1∑N∣qk∣2(1.19)
也就是数字信号处理中所说的Parseval定理,只不过在这些书籍中写作
1N∑n=1N∣xn∣2=∑k=1N∣ak∣2(1.20)\frac{1}{\text{N}}\sum\limits_{\text{n}=1}^{N}{{{\left| {{x}_{n}} \right|}^{2}}}=\sum\limits_{\text{k}=1}^{N}{{{\left| {{\text{a}}_{\text{k}}} \right|}^{2}}}\tag{1.20}N1n=1∑N∣xn∣2=k=1∑N∣ak∣2(1.20)
我更喜欢式(1.18)的形式,它具备了几何上的直观,可以看成是勾股定理的推广(连续的情形也一样,推导过程可以参看Cohen的书的第一篇附录),也是Bessel不等式的极限情况。ejθ{{e}^{j\theta }}ejθ也可以看出是两个正交分量的组合:ejθ=cosθ+jsinθ{{e}^{j\theta }}=\cos \theta +\text{j}\sin \thetaejθ=cosθ+jsinθ,它的范数的平方就是勾股定理:∥ejθ∥2=cos2θ+sin2θ=1{{\left\| {{e}^{j\theta }} \right\|}^{2}}={{\cos }^{2}}\theta +{{\sin }^{2}}\theta =1ejθ2=cos2θ+sin2θ=1,我们熟知的弹簧振子模型,模型的机械能E=12kx2+12mv2=(12kx)2+(12mv)2=a2+b2\text{E}=\frac{1}{2}k{{x}^{2}}+\frac{1}{2}m{{v}^{2}}={{\left( \sqrt{\frac{1}{2}k}x \right)}^{2}}+{{\left( \sqrt{\frac{1}{2}m}v \right)}^{2}}={{a}^{2}}+{{b}^{2}}E=21kx2+21mv2=(21kx)2+(21mv)2=a2+b2也是勾股定理的形式。Euler公式给我们启示,a2+b2{{a}^{2}}+{{b}^{2}}a2+b2可以写成(a+jb)(a−jb)\left( a+\text{j}b \right)\left( a-\text{j}b \right)(a+jb)(a−jb)即⟨a+jb,a+jb⟩\left\langle a+\text{j}b,a+\text{j}b \right\rangle⟨a+jb,a+jb⟩。弹簧振子模型在复平面上可以用向量(a+jb)\left( a+\text{j}b \right)(a+jb)描述,实轴是位置,虚轴是动量,位置a(振子的x)和速度b(振子的v)是相应的坐标值,它们的范的平方就是势能和动能,向量的范的平方就是机械能。振子的运动,就是向量绕原点做匀速圆周运动(二次型能量的系统,角频率ω=km\omega \text{=}\sqrt{\frac{k}{m}}ω=mk),动能势能此消彼长,而机械能守恒(在实数域,该守恒律的学习需要几堂课程,而在复数域却是一目了然的事情,当我们拥有了更高级的数学工具,可以统一处理很多以前费脑筋的事情)。也就是说a和b一对正交的同频的三角函数,就是Euler公式,是Helmholtz方程的解。兜兜转转又回来了。cosθ\cos \thetacosθ和sinθ\sin \thetasinθ是对偶,勾和股是对偶,位置和动量是对偶,势能和动能是对偶,对偶总是和共轭与正交联系在一起,时域和频域也可以看作是对偶,这些对于和振动打交道的我们来说是再熟悉不过的事情了。此外,倘若弹簧振子加上阻尼器,就形成一个弹簧阻尼模型,离散化得到AR模型,便可以在计算机上求解它的共振频率了,我曾经在提取某个动部件的共振频率的项目中用的就是这一方法。
利用(AH)−1=A{{({{\mathbf{A}}^{\text{H}}})}^{-1}}=\mathbf{A}(AH)−1=A的性质,Parseval定理还可以通过如下的数学游戏得出。
∥X∥2=XHX(1.21){{\left\| \mathbf{X} \right\|}^{2}}={{\mathbf{X}}^{\text{H}}}\mathbf{X}\tag{1.21}∥X∥2=XHX(1.21)
∥X∥2=qHAHAq(1.22){{\left\| \mathbf{X} \right\|}^{2}}={{\mathbf{q}}^{\text{H}}}{{\mathbf{A}}^{\text{H}}}\mathbf{Aq}\tag{1.22}∥X∥2=qHAHAq(1.22)
∥X∥2=qHq(1.23){{\left\| \mathbf{X} \right\|}^{2}}={{\mathbf{q}}^{\text{H}}}\mathbf{q}\tag{1.23}∥X∥2=qHq(1.23)
∥X∥2=∥q∥2(1.24){{\left\| \mathbf{X} \right\|}^{2}}={{\left\| \mathbf{q} \right\|}^{2}}\tag{1.24}∥X∥2=∥q∥2(1.24)
也可以用在中间插入AAH\mathbf{A}{{\mathbf{A}}^{\text{H}}}AAH
∥X∥2=XHAAHX(1.25){{\left\| \mathbf{X} \right\|}^{2}}={{\mathbf{X}}^{\text{H}}}\mathbf{A}{{\mathbf{A}}^{\text{H}}}\mathbf{X}\tag{1.25}∥X∥2=XHAAHX(1.25)
∥X∥2=(qHAHA)(AHAq)(1.26){{\left\| \mathbf{X} \right\|}^{2}}=\left( {{\mathbf{q}}^{\text{H}}}{{\mathbf{A}}^{\text{H}}}\mathbf{A} \right)\left( {{\mathbf{A}}^{\text{H}}}\mathbf{Aq} \right)\tag{1.26}∥X∥2=(qHAHA)(AHAq)(1.26)
∥X∥2=qHq{{\left\| \mathbf{X} \right\|}^{2}}={{\mathbf{q}}^{\text{H}}}\mathbf{q}∥X∥2=qHq
∥X∥2=∥q∥2{{\left\| \mathbf{X} \right\|}^{2}}={{\left\| \mathbf{q} \right\|}^{2}}∥X∥2=∥q∥2
无聊的时候,这些游戏常常帮我消磨时光,熬过无数孤独的日子。
让我们把Fourier变换对式(1.12)和式(1.15)放一起,
{AHX=qX=Aq(1.27)\left\{ \begin{matrix}
{{\mathbf{A}}^{\text{H}}}\mathbf{X}=\mathbf{q} \\
\mathbf{X}=\mathbf{Aq} \\
\end{matrix} \right.\tag{1.27}{AHX=qX=Aq(1.27)
应该很容易看出,这是一个简单的坐标变换,我们可以从几何上去理解。空间中的一个点,怎样去描述它呢?首先要选定一个原点,从原点到该点,一个向量就确定下来了。怎样去描述一个向量呢?要指定一组坐标系,用在各个坐标轴上投影值去描述。为了方便计算,坐标系往往选择正交的和单位长度的。非正交系当然也能用,但是有一堆“协变”和“逆变”的麻烦事。如果是正交的,那么合成法则就是简单得多的矩形了,正因为如此,我们往往要把非正交系正交化,所用的工具便是Schmidt正交化。——这一点是受到力学书籍的启发,我手头上的舒幼生的《力学》、丹尼尔《力学概论》,运动学部分基本上是在讲力在直角坐标系下如何如何、在极坐标系下如何让如何、在柱坐标系下又如何如何。
从坐标变换的角度看在,时域的X\mathbf{X}X是在坐标系E\mathbf{E}E(单位矩阵)下的坐标,频域的q\mathbf{q}q是在坐标系A\mathbf{A}A下的坐标,它两描述的是同一个量,是等价的,或者说是同一个量在不同坐标系下的描述(坐标系的选择不会对物理规律产生影响,简单来说,就是复数对加法构成群,非零复数对乘法构成群,保证了不同坐标系之间总是可以相互转换。这也是我们在力学的解题中总是可以根据难易程度随意选择坐标系的原因),因此我们上文说时域和频域是一对对偶(容易证明,X\mathbf{X}X、q\mathbf{q}q各自所属的空间都是完备且线性的,它们互为对偶空间)。在线性代数中,用“基”来定量表示一个坐标系,在量子力学中,把“基”称为表象,所谓的粒子的波粒二象性,是指同一个态矢量在不同的表象(坐标表象和动量表象)下的描述,上文说位置和动量是一对对偶,在这里也一样。套用量子力学常用的术语,EH{{\mathbf{E}}^{\text{H}}}EH是时域算符,AH{{\mathbf{A}}^{\text{H}}}AH是频域算符,于是有
{q=AHXX=EHX(1.28)\left\{ \begin{matrix}
\mathbf{q}\text{=}{{\mathbf{A}}^{\text{H}}}\mathbf{X} \\
\mathbf{X}\text{=}{{\mathbf{E}}^{\text{H}}}\mathbf{X} \\
\end{matrix} \right.\tag{1.28}{q=AHXX=EHX(1.28)
两个算符是不对易的,带来了不确定性的问题(海森堡不确定原理背后的原因也即在于此),导致我们在做STFT时,时域和频域的分辨率存在矛盾。
X→q\mathbf{X}\to \mathbf{q}X→q的过程是整个空间的一系列的坐标旋转和伸缩变换的级联。q→X\mathbf{q}\to \mathbf{X}q→X则是其逆过程,旋转和伸缩反方操作一遍就行了,就像电影的倒放一样。
让我们再看一眼式(1.19)或者式(1.24)——从坐标变换的角度去看,该式表示出了某种变换不变性(一个系统,知道哪些东西是不变的,比知道哪些东西是变化的更为关键)。在抽象的数学公式上表现出的是长度不变性,而在具象的物理世界比如信号处理领域表现的是能量不变性。从具体的事物中抽象出规律,在抽象中看到现实的身影,这种意识和能力是工程师需要留意的——所谓色不异空、空不异色,讲得就是这个道理。
用矩阵的工具去讲述Fourier变换,大多数工程师可能会不太适应——这无疑是和在学生时代我们所受到的常规训练有关,在学校的课程当中,线性代数不仅不是重点,而且非常突兀。事实上我本人也没有见过这样的讲述——一定有的,只是我的阅读范围极其有限——之所有选用矩阵,一开始是由于个人对矩阵较为熟悉,用起来顺手罢了,到后来还发现,新的工具带给我们新的启示。
比如,Fourier变换是无数个坐标变换中的一种。
比如,式(1.)和式(1.2),时域是N维,频域也是N维,为什么呢?如果从信号与系统或者分析学的角度,只能说公式推导结果就是这样,没什么好解释的,但是从矩阵的角度看,理解起来简直是理所当然又一目了然的:Fourier变换对作为“一对”,能变换过来又能变换回去,必然是N→N,没有升维,也没有降维,用集合的术语来讲,叫等势,用映射的术语来讲,叫双射。
比如,我们始终在强调对称性,细心的工程师会注意到,式(1.10):xn→qk{{x}_{n}}\to {{q}_{k}}xn→qk是分解,式(1.16):qk→xn{{q}_{k}}\to {{x}_{n}}qk→xn是合成。在分析学上是理所当然的,但是在坐标变换角度来看却要回答一个问题:既然我们说xn{{x}_{n}}xn是时域坐标,qk{{q}_{k}}qk是频域坐标,是平等的,等价的,两者之间的变换为什么不是相同的,或者为什么不是反过来的。
以前学习二维平面坐标旋转矩阵[cosθsinθ-sinθcosθ]\left[ \begin{matrix}
\cos \theta & \sin \theta \\
\text{-}\sin \theta & \cos \theta \\
\end{matrix} \right][cosθ-sinθsinθcosθ]的时候也向自己提出过类似的疑问。实际上其原因在于旋转后的坐标轴的单位向量(cosθsinθ)\left( \begin{matrix}
\cos \theta & \sin \theta \\
\end{matrix} \right)(cosθsinθ)和(-sinθcosθ)\left( \begin{matrix}
\text{-}\sin \theta & \cos \theta \\
\end{matrix} \right)(-sinθcosθ)是由原坐标轴的单位向量(10)\left( \begin{matrix}
1 & 0 \\
\end{matrix} \right)(10)和(01)\left( \begin{matrix}
0 & 1 \\
\end{matrix} \right)(01)所表示的。如果把表示关系倒过来,疑问便会得到解答。
Fourier变换对的问题是二维平面旋转矩阵问题的升级版,顺着这个思路,将会发现正变换和逆变换是相同的。具体地说,式(1.8)的基V^k{{\mathbf{\hat{V}}}_{\text{k}}}V^k改为
V^k=δ[n−k](n=1,2,3,...,N;k=1,2,3,...,N)(1.29){{\mathbf{\hat{V}}}_{\text{k}}}\text{=}\delta \left[ n-k \right]\left( n=1,2,3,...,N;k=1,2,3,...,N \right)\tag{1.29}V^k=δ[n−k](n=1,2,3,...,N;k=1,2,3,...,N)(1.29)
时域的基则变成
E^k=[1Ne−jkω01,1Ne−jkω02,...,1Ne−jkω0(N−1),1Ne−jkω0N]T(k=1,2,3,...,N)(1.30){{\mathbf{\hat{E}}}_{k}}={{\left[ \frac{1}{\sqrt{\text{N}}}{{e}^{-jk{{\omega }_{0}}1}},\frac{1}{\sqrt{\text{N}}}{{e}^{-jk{{\omega }_{0}}2}},...,\frac{1}{\sqrt{\text{N}}}{{e}^{-jk{{\omega }_{0}}\left( N-1 \right)}},\frac{1}{\sqrt{\text{N}}}{{e}^{-jk{{\omega }_{0}}N}} \right]}^{T}}\left( k=1,2,3,...,N \right)\tag{1.30}E^k=[N1e−jkω01,N1e−jkω02,...,N1e−jkω0(N−1),N1e−jkω0N]T(k=1,2,3,...,N)(1.30)
于是乎式(1.10)和式(1.16)的关系便倒过来了,完美地向我们展示出了它的对称性。
比如,变换是线性的,原点不变,单位向量模长不变(式(1.9)的结论),Parseval定理又是显然的了。至于式(1.20)出现的系数1N\frac{1}{\text{N}}N1,当然是单位向量模长拉伸了N\sqrt{\text{N}}N倍(式(1.7)的结论),从这个意义上讲式(1.20)应该写成
∑n=1N∣xn∣2=∑k=1N∣Nak∣2(1.31)\sum\limits_{n=1}^{N}{{{\left| {{x}_{n}} \right|}^{2}}}=\sum\limits_{k=1}^{N}{{{\left| \sqrt{N}{{a}_{k}} \right|}^{2}}}\tag{1.31}n=1∑N∣xn∣2=k=1∑NNak2(1.31)
的形式。人们用“保范性”来称呼这一性质,同时我们可以进一步追问,它具有“保角性”吗?不妨简单推导一下。
<q1,q2>=q1Tq2∗=(AHX1)T(AHX2)∗<{{\mathbf{q}}_{1}},{{\mathbf{q}}_{2}}>=\mathbf{q}_{1}^{T}\mathbf{q}_{2}^{*}={{({{\mathbf{A}}^{H}}{{\mathbf{X}}_{1}})}^{T}}{{({{\mathbf{A}}^{H}}{{\mathbf{X}}_{2}})}^{*}}<q1,q2>=q1Tq2∗=(AHX1)T(AHX2)∗
<q1,q2>=X1T(AH)T(AH)∗X2∗<{{\mathbf{q}}_{1}},{{\mathbf{q}}_{2}}>=\mathbf{X}_{1}^{T}{{({{\mathbf{A}}^{H}})}^{T}}{{({{\mathbf{A}}^{H}})}^{*}}\mathbf{X}_{2}^{*}<q1,q2>=X1T(AH)T(AH)∗X2∗
<q1,q2>=X1TAHAX2∗<{{\mathbf{q}}_{1}},{{\mathbf{q}}_{2}}>=\mathbf{X}_{1}^{T}{{\mathbf{A}}^{H}}\mathbf{AX}_{2}^{*}<q1,q2>=X1TAHAX2∗
<q1,q2>=X1TX2∗<{{\mathbf{q}}_{1}},{{\mathbf{q}}_{2}}>=\mathbf{X}_{1}^{T}\mathbf{X}_{2}^{*}<q1,q2>=X1TX2∗
<q1,q2>=<X1,X2>(1.32)<{{\mathbf{q}}_{1}},{{\mathbf{q}}_{2}}>=<{{\mathbf{X}}_{1}},{{\mathbf{X}}_{2}}>\tag{1.32}<q1,q2>=<X1,X2>(1.32)
由此我们知道答案是肯定的。看看式(1.2),我们应该能理解它确实属于保角变换。“保范性”决定了滤波器的幅频特性,“保角性”决定了滤波器的相频特性,并且可以得出,逆变换是正变换的伴随算子。
又比如,它提供了一个具备几何直观的理解不确定性的方法。当我们从数列的角度去看待x[n]x[n]x[n]时,便可以引入均值和方差的概念,方差越小即意味着在均值位置能量越集中。当我们切换到向量视角,均值位置对应的取样时刻n可以看做是一个坐标轴,方差越小意味着向量越靠近均值坐标轴,即与该坐标轴的内积越大。当方差小于某个临界值,在小下去,可以想象出在频域各坐标轴上的投影越分散,即方差越大。反之也一样。
这些都是车轱辘话,念叨过来又念叨过去,说得透透的,甚至有点烦了。
在R^N上的分解
FIR在我看来,是一个乘累加运算,是对输入x[n]x[n]x[n]的加权平均,是一个N×N→K\text{N}\times \text{N}\to \text{K}N×N→K的映射,是一个内积运算,是一个相关运算。而它的系数是一个有限长数列,是一个Euclidean空间RN{{\mathbb{R}}^{\text{N}}}RN上的向量。FIR的公式
y[n]=∑k=0Nbkx[n−k](2.1)y[n]=\sum\limits_{k=0}^{N}{{{b}_{k}}x[n-k]}\tag{2.1}y[n]=k=0∑Nbkx[n−k](2.1)
让我们忽略x[n]x[n]x[n]的动态输入,将式(2.1)改写成
y=∑k=0Nbkxk(2.2)y=\sum\limits_{k=0}^{N}{{{b}_{k}}{{x}_{\text{k}}}}\tag{2.2}y=k=0∑Nbkxk(2.2)
y=⟨X,B⟩=XTB∗(2.3)y=\left\langle \mathbf{X},\mathbf{B} \right\rangle ={{\mathbf{X}}^{\text{T}}}{{\mathbf{B}}^{*}}\tag{2.3}y=⟨X,B⟩=XTB∗(2.3)
之所以说它是一个相关运算(让我们忽略“时间反转”,那只是个先后顺序的对齐问题),从公式上看,当输入X\mathbf{X}X被缩放成单位长度的时候X~=X∥X∥\mathbf{\tilde{X}}=\frac{\mathbf{X}}{\sqrt{\left\| \mathbf{X} \right\|}}X~=∥X∥X,y的大小取决于X~\mathbf{\tilde{X}}X~长得像不像B\mathbf{B}B。两个物体,a和b,长得像不像,应该怎么去判断呢?张三有张一寸照片,人们说,照片像张三。可是照片那么小,张三脑袋那么大,怎么就像了呢?猪脑袋和张三的脑袋一样大,为什么没有人说张三长得像猪呢?在人们看来,这是不言自明的事情,也许我比较笨,被这个问题困扰了许久。现在,我尝试着给出一个答案,a和b长得像不像,取决于它们之间的夹角θ\thetaθ,和它们的大小无关。
X\mathbf{X}X和B\mathbf{B}B的夹角cosθ=⟨X,B⟩∥X∥∥B∥\cos \theta =\frac{\left\langle \mathbf{X},\mathbf{B} \right\rangle }{\sqrt{\left\| \mathbf{X} \right\|}\sqrt{\left\| \mathbf{B} \right\|}}cosθ=∥X∥∥B∥⟨X,B⟩,概率论当中的相关系数ρ(X,Y)=cov(X,Y)var(X)var(Y)\rho (X,Y)=\frac{\operatorname{cov}(X,Y)}{\sqrt{\operatorname{var}(X)}\sqrt{\operatorname{var}(Y)}}ρ(X,Y)=var(X)var(Y)cov(X,Y)就是cosθ\cos \thetacosθ,Schwarz不等式添加上∣cosθ∣\left| \cos \theta \right|∣cosθ∣就是等式。cosθ\cos \thetacosθ越大,X\mathbf{X}X和B\mathbf{B}B越像,当我们用X~\mathbf{\tilde{X}}X~取代X\mathbf{X}X的时候,cosθ\cos \thetacosθ和⟨X,B⟩\left\langle \mathbf{X},\mathbf{B} \right\rangle⟨X,B⟩相差的也就一个系数而已,可以不做区分。我曾经用求解相关函数的方法解决过一个轮盘的齿圈识别问题,也算是学有所用。
从这里我们可以理解为什么式(1.2)的e−jkω0n{{e}^{-jk{{\omega }_{0}}n}}e−jkω0n为什么会有个负号。在分析学中,我们知道这是为了保证向量自己和自己内积的正性,即⟨X,X⟩≥0\left\langle \mathbf{X},\mathbf{X} \right\rangle \ge 0⟨X,X⟩≥0。现在我们知道,求两向量的角度差当然要相减。在对代数特别是群论的学习中,我们会知道实数加群R\mathbb{R}R、复数乘法群{eiθ∣θ∈R}\left\{ {{e}^{i\theta }}|\theta \in \mathbb{R} \right\}{eiθ∣θ∈R}存在同构的关系,<x,ejkω0n><x,{{e}^{jk{{\omega }_{0}}n}}><x,ejkω0n>运算蕴含了复数xxx和ejkω0n{{e}^{jk{{\omega }_{0}}n}}ejkω0n的夹角,即实数∠x\angle x∠x和实数kω0n{k{{\omega }_{0}}n}kω0n的差,也就是在R\mathbb{R}R上是(∠x)+(−kω0n)(\angle x)+(-k{{\omega }_{0}}n)(∠x)+(−kω0n),对应在{eiθ∣θ∈R}\left\{ {{e}^{i\theta }}|\theta \in \mathbb{R} \right\}{eiθ∣θ∈R}上是ei((∠x)+(−kω0n)){{e}^{i((\angle x)+(-k{{\omega }_{0}}n))}}ei((∠x)+(−kω0n))。
讨论完X\mathbf{X}X在B\mathbf{B}B上的投影⟨X,B⟩\left\langle \mathbf{X},\mathbf{B} \right\rangle⟨X,B⟩,下面我们讨论一下X\mathbf{X}X在B\mathbf{B}B的各正交分量的上的投影。FIR的系数B\mathbf{B}B可以依照式(1.17)进行分解。
B=∑k=1NpkV^k(2.4)\mathbf{B}=\sum\limits_{k=1}^{N}{{{\text{p}}_{\text{k}}}{{{\mathbf{\hat{V}}}}_{k}}}\tag{2.4}B=k=1∑NpkV^k(2.4)
带入式(2.3)可得
y=XT∑k=1N(pkV^k)∗y={{\mathbf{X}}^{\text{T}}}\sum\limits_{k=1}^{N}{{{\left( {{\text{p}}_{\text{k}}}{{{\mathbf{\hat{V}}}}_{k}} \right)}^{*}}}y=XTk=1∑N(pkV^k)∗
y=∑k=1NXT(pkV^k)∗(2.5)y=\sum\limits_{k=1}^{N}{{{\mathbf{X}}^{\text{T}}}{{\left( {{\text{p}}_{\text{k}}}{{{\mathbf{\hat{V}}}}_{k}} \right)}^{*}}}\tag{2.5}y=k=1∑NXT(pkV^k)∗(2.5)
令Bk=pkV^k{{\mathbf{B}}_{k}}={{\text{p}}_{\text{k}}}{{\mathbf{\hat{V}}}_{k}}Bk=pkV^k,则
y=∑k=1NXTBk∗(2.6)y=\sum\limits_{k=1}^{N}{{{\mathbf{X}}^{\text{T}}}\mathbf{B}_{k}^{*}}\tag{2.6}y=k=1∑NXTBk∗(2.6)
y=∑k=1N⟨X,Bk⟩(2.7)y=\sum\limits_{k=1}^{N}{\left\langle \mathbf{X},{{\mathbf{B}}_{k}} \right\rangle }\tag{2.7}y=k=1∑N⟨X,Bk⟩(2.7)
也就是说,X\mathbf{X}X在B\mathbf{B}B上的投影等效于X\mathbf{X}X在B\mathbf{B}B的正交分量上的投影的和。
由式(2.5)还可得
y=∑k=1N⟨X,pkV^k⟩y=\sum\limits_{k=1}^{N}{\left\langle \mathbf{X},{{\text{p}}_{\text{k}}}{{{\mathbf{\hat{V}}}}_{k}} \right\rangle }y=k=1∑N⟨X,pkV^k⟩
y=∑k=1Npk∗⟨X,V^k⟩y=\sum\limits_{k=1}^{N}{\text{p}_{k}^{*}\left\langle \mathbf{X},{{{\mathbf{\hat{V}}}}_{k}} \right\rangle }y=k=1∑Npk∗⟨X,V^k⟩
y=∑k=1Npk∗qk(2.8)y=\sum\limits_{k=1}^{N}{\text{p}_{k}^{*}{{\text{q}}_{k}}}\tag{2.8}y=k=1∑Npk∗qk(2.8)
通俗来讲,FIR滤波器,描述的就是什么信号能通过,什么信号不能通过,就是看长得像不像(相关运算,或者说求夹角),就是把B\mathbf{B}B在频域上分解,从频域上观察pk{{\text{p}}_{\text{k}}}pk。其实我们还可以把各个的Bk=pkV^k{{\mathbf{B}}_{k}}={{\text{p}}_{\text{k}}}{{\mathbf{\hat{V}}}_{k}}Bk=pkV^k画出来,这样频域又转回时域了,由式(2.7)可知,还是一个长得像不像的问题。人们常说,学会从时域和频域的角度看问题,式(2.4)告诉我们,要从时域看到频域的信息,而式(2.7)告诉我们要从频域看到时域的图景。
事实上,X\mathbf{X}X和B\mathbf{B}B可以同时做正交分解,一条公式中就出现两个式(1.14)所示的大块头矩阵,哪位工程师有兴趣,可以试着推导一下。
当一个FIR系统运行起来,X\mathbf{X}X的值在不断地更新,如果X\mathbf{X}X的周期也是N,那么qk{{\text{q}}_{k}}qk的周期也是N,输出y的周期也是N。在某些情况下,X\mathbf{X}X的周期是可知的,曾经我在对市电做FFT的项目中就利用到了这一点。市电的频率大约50HZ,其谐波是50HZ的倍数,尽管存在其他频率,占比是小到可以忽略不计的,那么在对市电进行采样的时候,样本点绘制的波形就应该是市电波形的整数倍,规避掉频谱泄漏的问题(如果采样频率和市电频率都是理想的)。
更多的情况,X\mathbf{X}X的周期并不是N的倍数,我们需要知道B\mathbf{B}B对所有频率输入的响应。为了达到这一目的,我们需要在l2{{l}^{2}}l2空间讨论FIR。
在l^2上的分解
l2{{l}^{2}}l2上的向量是可数无穷的,计算机系统是当然不能存储无穷的数据。但是有限维的B\mathbf{B}B却可以等价地表示一个无穷维的向量,只要B\mathbf{B}B后续的元素都是0。我把FIR当成IIR的一个特例——N之后的元素都是0。符号太多了,让我们还是继续用B\mathbf{B}B代表这一无穷维向量吧!
在有限维上使用的是Fourier级数,可数无穷维上使用的是离散时间Fourier变换(DTFT)
X(ω)=∑n=−∞+∞x[n]e−jωn(3.1)X\left( \omega \right)=\sum\limits_{n=-\infty }^{+\infty }{x\left[ n \right]{{e}^{-j\omega n}}}\tag{3.1}X(ω)=n=−∞∑+∞x[n]e−jωn(3.1)
B\mathbf{B}B可以在单位向量的移位向量δ[n−k]\delta \left[ n-k \right]δ[n−k]组成的基上展开,容易看出,这组基也是正交的,它的运算满足内积的条件
(1)⟨X,X⟩≥0\left\langle \mathbf{X},\mathbf{X} \right\rangle \ge 0⟨X,X⟩≥0,⟨X,X⟩=0\left\langle \mathbf{X},\mathbf{X} \right\rangle =0⟨X,X⟩=0当且仅当X=0\mathbf{X}=\mathbf{0}X=0
(2)⟨X,Y⟩=⟨Y,X⟩∗\left\langle \mathbf{X},\mathbf{Y} \right\rangle ={{\left\langle \mathbf{Y},\mathbf{X} \right\rangle }^{*}}⟨X,Y⟩=⟨Y,X⟩∗
(3)⟨aX,Y⟩=a⟨X,Y⟩\left\langle a\mathbf{X},\mathbf{Y} \right\rangle =a\left\langle \mathbf{X},\mathbf{Y} \right\rangle⟨aX,Y⟩=a⟨X,Y⟩
(4)⟨X+Y,Z⟩=⟨X,Z⟩+⟨Y,Z⟩\left\langle \mathbf{X}+\mathbf{Y},\mathbf{Z} \right\rangle =\left\langle \mathbf{X},\mathbf{Z} \right\rangle +\left\langle \mathbf{Y},\mathbf{Z} \right\rangle⟨X+Y,Z⟩=⟨X,Z⟩+⟨Y,Z⟩
并且还是标准的
⟨δ[n−k],δ[n−r]⟩={1k=r0k≠r(3.2)\left\langle \mathbf{\delta }\left[ n-k \right],\mathbf{\delta }\left[ n-r \right] \right\rangle =\left\{ \begin{matrix}
1\begin{matrix}
{} & k=r \\
\end{matrix} \\
0\begin{matrix}
{} & k\ne r \\
\end{matrix} \\
\end{matrix} \right.\tag{3.2}
⟨δ[n−k],δ[n−r]⟩={1k=r0k=r(3.2)
那么在上面的展开便是
B=[b0,b1,...,bN−1]T=∑k=0N-1bkδ[n−k]\mathbf{B}={{\left[ {{b}_{0}},{{b}_{1}},...,{{b}_{N-1}} \right]}^{T}}=\sum\limits_{k=0}^{\text{N-1}}{{{b}_{k}}\mathbf{\delta }[n-k]}B=[b0,b1,...,bN−1]T=k=0∑N-1bkδ[n−k]
bkδ[n−k]{{b}_{k}}\mathbf{\delta }\left[ n-k \right]bkδ[n−k]的DTFT很容易得出:
Xk(ω)=bke−jωk(3.3){{X}_{k}}\left( \omega \right)={{b}_{k}}{{e}^{-j\omega k}}\tag{3.3}Xk(ω)=bke−jωk(3.3)
根据内积空间的线性性质,B\mathbf{B}B的DTFT是
XB(ω)=∑k=0N−1bke−jωk(3.4){{X}_{B}}\left( \omega \right)=\sum\limits_{k=0}^{N-1}{{{b}_{k}}{{e}^{-j\omega k}}}\tag{3.4}XB(ω)=k=0∑N−1bke−jωk(3.4)
上面一系列叙述,都默认了一个前提,系统有一个好的特性,在信号处理中,称作时不变(我把它看成一个物理系统的时间平移不变性)。这当然是一个理想化的特性,大千世界,有形的物体,总熬不过时间,无可奈何的宿命。幸运的是,很多情况下在允许的时间内,这种不变性可以认为是存在的,不幸的是,总有些情况系统是时变的,我在做共振频率提取的课题时,就必须考虑到随着时间的推移,机械件的磨损带来的一系列问题。
IIR
FIR和IIR往往出现在离散系统的讨论中,我却常常将这一概念延伸到连续系统当中去(连续系统的响应都是无限的吧?也就不会出现FIR和IIR的划分吧?)。一个连续系统,如前所述,我们常常用微分方程去抽象,它的冲激响应几乎都是无限长的,离散化之后的数列B=b[n](n∈N)\mathbf{B}=b\left[ n \right]\left( n\in N \right)B=b[n](n∈N)长度也是无限的,也就是实际上的IIR。工程师常常面临这样的一个话题,如何在一个不具备无限存储能力的计算机上实现IIR。
有两个手段可以去实现,一是做一个阉割版的IIR——用足够长的FIR去近似。一个稳定的因果的系统的δ(t)\delta \left( \text{t} \right)δ(t)输入的响应看成是L2[0∞]{{\text{L}}^{2}}[\begin{matrix}
0 & \infty \\
\end{matrix}]L2[0∞]上的函数h(t)\text{h}\left( t \right)h(t),显然的,limt→∞ h(t)=0\underset{t\to \infty }{\mathop{\lim }}\,\text{h(t)=0}t→∞limh(t)=0,即b[n]→0b\left[ n \right]\to 0b[n]→0,这个数列是收敛到0的:
对任意的ε>0\varepsilon >0ε>0,存在N,使得n>N\text{n}>\text{N}n>N,∣b[n]−0∣<ε\left| \text{b}\left[ \text{n} \right]-0 \right|<\varepsilon∣b[n]−0∣<ε
根据这一性质,选择一个足够小的ε\varepsilonε以及合适的长度N,使得在工程上b[N]\text{b}\left[ \text{N} \right]b[N]后面的元素可忽略不计。在一个主动降噪的系统辨识项目中我们便是这么做的。
第二个手段是更加常用的,即引入反馈,或者依据系统微分方程做离散化,即
∑r=0Mary[n−r]=∑k=0Nbkx[n−k](4.1)\sum\limits_{\text{r}=0}^{M}{{{a}_{r}}y[n-r]}=\sum\limits_{k=0}^{N}{{{b}_{k}}x[n-k]}\tag{4.1}r=0∑Mary[n−r]=k=0∑Nbkx[n−k](4.1)
如前文所述,我将FIR看成是特殊的IIR,其实,我也将IIR看成是特殊的FIR,y[n]y[n]y[n]可以递归展开成x[n−r]x[n-r]x[n−r]和y[k](k<0)y[k](k<0)y[k](k<0)的组合,生成新的bk{{b}_{k}}bk,系统初始值一般都是0,即y[k]=0(k<0)y[k]=0(k<0)y[k]=0(k<0),那么y[n]y[n]y[n]便可以由x[n−r]x[n-r]x[n−r]和bk{{b}_{k}}bk唯一地决定了。
2727

被折叠的 条评论
为什么被折叠?



