使用矩阵运算替代 for 循环实现信号的DTFT1. 问题的引入2. 使用矩阵运算实现DTFT3. 利用MATLAB绘制DTFT的幅频特性曲线4. 小结

225 阅读 0 评论 149 点赞

我是靠谱客的博主无聊服饰，最近开发中收集的这篇文章主要介绍使用矩阵运算替代 for 循环实现信号的DTFT1. 问题的引入2. 使用矩阵运算实现DTFT3. 利用MATLAB绘制DTFT的幅频特性曲线4. 小结，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

1. 问题的引入

最近的数字信号处理实验中有这样的一道题

对 $x_a(t)=e^{-1000|t|}$ 信号进行采样分析，讨论信号的抽样的最低频率，验证奈奎斯特采样定理。要求对信号 $x_a(t)$ 分析，确定采样频率，求出采样后的序列 $x(n)$ 的频谱，比较不同采样频率情况下（满足和不满足采样定理） FT频谱 $x(jOmega)$ 和 DTFT频谱 $x(e^{jomega})$ 关系。

实验思路如下。首先对信号进行截断，因为计算机只能分析有限长的序列，所以需要截断，截断的依据是当幅值衰减到小于某一阈值（如 $e^{-5}$ ）时截断；其次对截断后的连续非周期信号进行傅里叶变换（FT），这里通过 MATLAB 的 Fourier 命令可以直接得到；得到原信号的FT之后根据频谱判断截止频率（最高频率）大概在什么位置，从而确定采样信号的频率（根据香农采样定理，采样频率是最高频率的二倍）；最后对采样后的信号进行DTFT。因为在 MATLAB 中没有直接求取 DTFT 的函数，所以需要自己编写。而本文的主要灵感也来自于自己编写DTFT的过程。

假设我们现在已经得到了经过抽样之后的离散非周期序列 $x(n)$ ，现在要对其进行 DTFT 变换。由于 DTFT 的定义式如下

X(e^{jomega})=Sigma_{n=-infty}^{infty}x(n)e^{-jomega n}tag{1}label{1}

$X(e^{jomega})=Sigma_{n=-infty}^{infty}x(n)e^{-jomega n}tag{1}label{1}$
所以如果从定义式下手，需要对上述的

x(n)e−jω $x(n)e^{-jomega}$ 在使得

x(n) $x(n)$ 有意义上的区间内求和，也就是说是有限项累加求和的形式，这种形式的计算往往是通过 for 循环实现的。但是无论是在MATLAB 还是在 python 的 Numpy 计算库上，矩阵相乘的计算速度又要高于 for 循环的计算速度，所以将 for 循环改写成矩阵相乘的形式也是一件很应该做的事情，在机器学习中这个常常叫做“向量化”。

2. 使用矩阵运算实现DTFT

假设 $x(n)$ 是一个5 点的序列，即

[x (0) x (1) x (2) x (3) x (4)]

$begin{bmatrix} x(0) & x(1) & x(2) & x(3) & x(4) end{bmatrix}$
根据公式 (1) 可知，

x(n) $x(n)$ 的 DTFT 为

X (e j ω) = x (0) e - j ω 0 + x (1) e - j ω 1 + x (2) e - j ω 2 + x (3) e - j ω 3 + x (4) e - j ω 4

$X(e^{jomega}) =x(0)e^{-jomega 0} +x(1)e^{-jomega 1}+x(2)e^{-jomega 2}+x(3)e^{-jomega 3}+x(4)e^{-jomega 4}$
上式虽然是 5 项相加，但是实际上一个关于

ω $omega$ 的函数。这 5 项求和的过程可以使用 for 进行求和得到，当然也可以通过矩阵相乘的形式获得

X (e j ω) = [x (0) x (1) x (2) x (3) x (4)] * ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ e - j ω 0 e - j ω 1 e - j ω 2 e - j ω 3 e - j ω 4 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$X(e^{jomega}) = begin{bmatrix} x(0) & x(1) & x(2) & x(3) & x(4) end{bmatrix} * begin{bmatrix} e^{-jomega 0}\ e^{-jomega 1}\ e^{-jomega 2}\ e^{-jomega 3}\ e^{-jomega 4} end{bmatrix}$
通过上式的矩阵相乘就可以替代了一个 for 循环。上面的5点序列是一个行信号，在很多时候我们处理的是咧信号，如

[x (0) x (1) x (2) x (3) x (4)] T

$begin{bmatrix} x(0) & x(1) & x(2) & x(3) & x(4) end{bmatrix}^T$
因此

X (e j ω) = [e - j ω \times 0 e - j ω \times 1 e - j ω \times 2 e - j ω \times 3 e - j ω \times 4] * ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ x (0) x (1) x (2) x (3) x (4) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

$X(e^{jomega}) = begin{bmatrix} e^{-jomega ×0} & e^{-jomega ×1} & e^{-jomega ×2} & e^{-jomega ×3} & e^{-jomega ×4} end{bmatrix} * begin{bmatrix} x(0)\ x(1)\ x(2)\ x(3)\ x(4) end{bmatrix}$
上述的过程是显而易见的，因为

(AB)T=BTAT $(AB)^T=B^TA^T$

3. 利用MATLAB绘制DTFT的幅频特性曲线

很明显通过 DTFT 得到的表达式是一个关于 $omega$ 以 $2pi$ 为周期的周期连续函数， $omega$ 的取值范围为 $[-infty, +infty]$ ，因此在利用 MATLAB 在绘制DTFT的幅频特性曲线时，可以首先去一个或两个周期进行绘图（如取2个周期的话，应该绘制[ $-2pi,2pi$ ]区间内的图像即可）。但实际上 MATLAB 无法直接绘制出连续函数的图像，是通过求取一系列的函数值，在把这些值连接起来得到的函数曲线。所以我们要绘制出 $X(e^{jomega})$ 的图像，还需要求取 $X(e^{jomega})$ 若干个点的值。理论上当然是取的点数越多绘制的图像越接近于真实图像，但是为了书写方便，在这里我们以 5 点为例进行说明。

假设我们求取一个周期上的 5 个点 $omega_0，omega_1, omega_2,omega_3, omega_4$ 。则最后 MATLAB 计算得到的如下

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ x (0) e - j ω 0 \times 0 + x (1) e - j ω 0 \times 1 + x (2) e - j ω 0 \times 2 + x (3) e - j ω 0 \times 3 + x (4) e - j ω 0 4 x (0) e - j ω 1 \times 0 + x (1) e - j ω 1 \times 1 + x (2) e - j ω 1 \times 2 + x (3) e - j ω 1 \times 3 + x (4) e - j ω 1 4 x (0) e - j ω 2 \times 0 + x (1) e - j ω 2 \times 1 + x (2) e - j ω 2 \times 2 + x (3) e - j ω 2 \times 3 + x (4) e - j ω 2 4 x (0) e - j ω 3 \times 0 + x (1) e - j ω 3 \times 1 + x (2) e - j ω 3 \times 2 + x (3) e - j ω 3 \times 3 + x (4) e - j ω 3 4 x (0) e - j ω 4 \times 0 + x (1) e - j ω 4 \times 1 + x (2) e - j ω 4 \times 2 + x (3) e - j ω 4 \times 3 + x (4) e - j ω 4 4 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ T

$begin{bmatrix} x(0)e^{-jomega_0× 0} +x(1)e^{-jomega_0× 1}+x(2)e^{-jomega_0× 2}+x(3)e^{-jomega_0× 3}+x(4)e^{-jomega_0 4} \ x(0)e^{-jomega_1× 0} +x(1)e^{-jomega_1× 1}+x(2)e^{-jomega_1× 2}+x(3)e^{-jomega_1× 3}+x(4)e^{-jomega_1 4} \ x(0)e^{-jomega_2× 0} +x(1)e^{-jomega_2× 1}+x(2)e^{-jomega_2× 2}+x(3)e^{-jomega_2× 3}+x(4)e^{-jomega_2 4} \ x(0)e^{-jomega_3× 0} +x(1)e^{-jomega_3× 1}+x(2)e^{-jomega_3× 2}+x(3)e^{-jomega_3× 3}+x(4)e^{-jomega_3 4} \ x(0)e^{-jomega_4× 0} +x(1)e^{-jomega_4× 1}+x(2)e^{-jomega_4× 2}+x(3)e^{-jomega_4× 3}+x(4)e^{-jomega_4 4} end{bmatrix}^T$

注意上面是矩阵的转置，所以实际上是一个 1×5 的矩阵，它的每一列都包含 $begin{bmatrix} x(0) & x(1) & x(2) & x(3) & x(4) end{bmatrix}$ ，上述矩阵等价于如下形式。

[x (0) x (1) x (2) x (3) x (4)] \times e x p （ - j ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 01234 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ [ω 0 ω 1 ω 2 ω 3 ω 4] ）

$begin{bmatrix} x(0) & x(1) & x(2) & x(3) & x(4) end{bmatrix}×exp（{-j} begin{bmatrix} 0\1\2\3\4 end{bmatrix} begin{bmatrix} omega_0&omega_1&omega_2&omega_3&omega_4 end{bmatrix}）$
这里有三点需要注意的地方，首先上面的

ω $omega$ 我们取了5个，实际上是在所讨论的周期内等间距取的5个

ω $omega$ 的值；其次这里取的

ω $omega$ 的数量与原信号序列长度相同只是一个巧合，实际上

ω $omega$ 的数量是任意的，只要能满足取这么多的点可以画出图来就好；最后是当信号是列向量的时候，同样是根据

(AB)T=BTAT $(AB)^T=B^TA^T$ 计算。

4. 小结

在上式的计算过程中总共涉及了3个向量，即 $x(n)$ 、 $n$ 和 $omega$ 。其中 $x(n)$ 是一个 1×n 的向量，求和的点数 n 也是一个 1×n 的向量，参数 $omega$ 一共取 b 个点，所以是一个 1 × b 的向量，最后我们要得到的是一个 1 × b 的向量。有多少个 $omega$ ，就相当于有多少列，所以明显参数 $omega$ 应该放在最后面，而 n 与 $omega$ 都是位于指数项的参数，所以它们应该放在一起相乘，得到的矩阵再 × $e^{-j}$ 得到一个矩阵（该矩阵明显是一个 n × b 的矩阵），将 $x(n)$ 右乘该矩阵便可以得到最终的结果。这样就将一个双层循环问题转化为三个向量的相乘的问题。