数字信号频域变换分析

2022-02-02

学习提升 / 信号处理与雷达 / 信号与系统

字数统计: 13k | 阅读时长≈ 52 分钟

一、数字信号处理中常用矩函数计算和其物理意义

(一) 常用矩函数

1. 均值

均值表示信号中直流分量的大小，用$E(x)$表示。对于高斯白噪声信号而言，它的均值为0，所以它只有交流分量。

$E(x) = \bar x = \frac{x_1 + x_2 + \cdots +x_n}{n} = \frac{\sum_{i=1}^{n}x_i}{n}$

2. 均值的平方

均值的平方，用${E^2(x)}$表示，它表示的是信号中直流分量的功率。

3. 均方值

均方值表示信号平方后的均值，用$E(x^2)$表示。均方值表示信号的平均功率。信号的平均功率 = 信号交流分量功率信号直流分量功率。

4. 均方根值

用RMS（root meansquare），即均方值的开根号。

$X_{rms} = \sqrt{\frac{\sum_{i=1}^{n}X_i^2}{N}}$

5. 均方差

均方差（mean squareerror），用MSE表示。均方差是各数据偏离真实值的距离平方和的平均数，也即误差平方和的平均数，计算公式形式上接近方差，它的开方叫均方根误差，均方根误差才和标准差形式上接近。均方差有时候被认为等同于方差。

$X_{MSE} = \frac{1}{n}$

6. 均方根误差

均方根误差用RMSE（root mean squareerror）表示。它是观测值与真值偏差的平方和观测次数n比值的平方根，在实际测量中，观测次数n总是有限的，真值只能用最可信赖（最佳）值来代替。方根误差对一组测量中的特大或特小误差反映非常敏感，所以均方根误差能够很好地反映出测量的精密度。均方根误差有时被认为是标准差。

7. 方差

方差用variance或deviation 或Var表示。方差描述信号的波动范围，表示信号中交流分量的强弱，即交流信号的平均功率。

$D(X) = E[(X-E(X))^2]$ $s^2 = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar x)^2 ——或者另一种形式——s^2 = \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar x)^2$

注意上面除以的是n-1，只有这样由样本值估计出的方差才是无偏的，即上面式子的期望才是X的方差。但是有的地方也有用除以n来表示方差，只不过这样求出的结果不是方差的无偏估计，计算结果的数学期望并不是X的方差，而是X方差的倍。

8. 标准差

标准差（Standard Deviation）用σ表示，有的时候标准差又可以被称为均方根误差RMSE。标准差是各数据偏离平均数的距离的平均数，它是离均差平方和平均后的方根，用σ表示，标准差能反映一个数据集的离散程度。

标准差σ，反映了测量数据偏离真实值的程度，σ越小，表示测量精度越高，因此可用σ作为评定这一测量过程精度的标准。

$S = \sqrt{\frac{\sum_{i=1}^{n}(x_i-\bar x)^2}{n}}$

有了方差为什么要使用标准差？标准差比方差有什么优势？

因为方差与我们要处理的数据的量纲是不一致的，虽然能很好的描述数据与均值的偏离程度，但是处理结果是不符合我们的直观思维的。

(二) 总结

（1）总的来说，均方差，均方根误差和方差，标准差是不能够等同的，尽管它们的公式相似。我们需要从真实值和均值之间的关系来区分它们

（2）对于方差和标准差而言，它们反映的是数据序列与均值的关系。

（3）对于均方差和均方根误差而言，它们反映的是数据序列与真实值之间的关系。

说明

文章转自：新浪博客/未注明作者，感谢原作者的辛勤付出，如涉及版权，请联系我删除！ http://www.360doc.com/content/18/0715/14/908538_770556276.shtml

二、数字域频率与模拟频率

(一) 二者关系

重点：数字域频率是模拟角频率相对于采样频率的归一化。

1. 第一种理解：

从”连续时间正弦信号与离散时间正弦信号“的关系来理解。

敲黑板，重点是

2. 第二种理解

从”采样信号的傅里叶变换与离散时间信号的傅里叶变换“的关系来理解。从时域采样出发，先看时域上离散时间信号$x(n)$的自变量$n$与连续时间信号$x(t)$的自变量$t$的关系，如下图。

在频域上模拟信号、采样信号、和采样序列(也就是离散时间信号)三者的频谱的关系，如下图所示。看看横轴，时域上，横轴由$t$变成$n$，是按照$n=t/T$的规则；而频域上，横轴由$\Omega$(模拟角频率)变成了$\omega$(数字域频率)，按照$\omega = \Omega × T$的规则。

数字域频率$\omega$的最高频率是$\pi$，这一点可以结合采样定理来理解。采样间隔$T$，意味着采样角频率为$2\pi/T$，根据采样定理，最高的频率是采样频率的一半，也就是$\pi/T$，再将它转换为数字域频率（(乘以$T$)，不就是$\pi$吗?

3. 第三种理解

从z平面与s平面的映射关系上来理解。

(二) 总结

不管从哪种角度、哪种方式，殊途同归，结论是一样的，那就是——数字域频率是模拟角频率相对于采样频率的归一化。这一结论，是我们能够正确解读采样后的信号它的频谱分析结果的前提。
模拟角频率和数字域频率是数字信号处理中非常重要的两个概念，因为经常需要将模拟信号离散化后再进行频域分析，那么，得到的是数字域频率，必须正确转换为模拟角频率后，才能得到分析对象——模拟信号的频域信息。

说明

文章转自：个人图书馆/未注明作者，感谢原作者的辛勤付出，如涉及版权，请联系我删除！ http://www.360doc.com/content/19/0611/18/908538_841808354.shtml

三、离散时间傅里叶变换DTFT

首先要需要理解英文缩写DTFT的含义：离散时间DT、傅里叶变换FT，连起来就是DTFT。注意第二个T (也就是Time)，至关重要。它是“离散D”这个特性所描述的主体：也就是，时间是离散的。

(一) 从z变换到DTFT

重点1：从“单位圆上的z变换“这个角度来理解DTFT正变换的定义。
重点2：理解DTFT的周期性——离散时间信号的傅里叶变换都是以$2\pi$为周期的。

1. DTFT的正变换：

$\text{DTFT}[x(n)] = X(e^{j\omega}) = \sum_{n = -\infty}^{\infty}x(n)e^{-j\omega n}$

2. DTFT的基本性质：

$X(e^{j\omega})$是以$\omega$为自变量的连续函数；

$X(e^{j\omega})$是以$2\pi$为周期的周期函数；

3. DTFT存在的充分条件：

$x(n)绝对可和：\sum_{n = -\infty}^{\infty}|x(n)| < \infty$

(二) DTFT变换对及物理含义

重点1：理解离散信号频谱的物理概念及特点。
重点2：会求几种常用信号的DTFT(单位样值信号、矩形脉冲、指数衰减信号、离散Sa函数等)。

1. DTFT公式整理

2. 物理含义

(1) $x(n)$可以表示成复指数信号的线性组合(只需一个$2\pi$区间内的频率)；

(2) $X(e^{j\omega})$表示了$x(n)$中各个频率分量的相对大小及位置，称为$x(n)$的频谱(以$2\pi$为周期)。

$X(e^{j \omega})=|X(e^{j \omega})| e^{j \angle X(e^{j \omega})}$

说明

文章转自：个人图书馆/未注明作者，感谢原作者的辛勤付出，如涉及版权，请联系我删除！ http://www.360doc.com/content/19/0611/18/908538_841808425.shtml

四、离散傅里叶级数DFS

(一) 离散时间周期信号的DFS

重点1：对照连续时间周期信号的FS的思想，理解离散时间周期信号的FS，二者的相同(离散谱)和不同之处。
重点2：理解离散时间周期信号的FS展开式为何只有N项，理解离散时间周期信号频谱的周期性。

需要说明的是，DFS有两种表示形式，如下面两个式子，这两种方式的区别在于，1/N的系数在正变换(FS系数求解式)中，还是反变换(级数展开式)中。没有实质区别，只是差一个常数N。这里均采用第二种表示形式，即将1/N的系数放在反变换中。

$第一种：\tilde{X} _ k = \frac{1}{N}\sum _ {n = 0}^{N-1}\tilde{x} _ ne^{-j\frac{2\pi}{N}kn}\text{ } - - - - \text{ }\tilde{x}(n) = \sum _ {k=0}^{N-1}\tilde{X} _ ke^{j\frac{2\pi}{N}kn}$ $第二种：\tilde{X}(k) =\sum _ {n = 0}^{N-1}\tilde{x} _ ne^{-j\frac{2\pi}{N}kn} \text{ } - - - -\text{ }\tilde{x}(n) = \frac{1}{N}\sum _ {k=0}^{N-1}\tilde{X}(k)e^{j\frac{2\pi}{N}kn}$

(二) 离散时间周期信号的频谱

重点1：掌握离散时间周期信号频谱的特点——离散性、谐波性、周期性；
重点2：会求常用离散时间周期信号的DFS(见例题1和例题2)。

1. 两类典型的求解DFS的题目：

第一类：周期信号直接以正弦、余弦之和形式给出(如例1)
方法：与标准形式的FS对照，直接得出FS的系数（即频谱)。
第二类：一般的周期信号(如例2)
方法：用FS的系数求解公式。

2. DFS与DTFT的关系

重点：理解DFS与DTFT的关系——周期序列的离散谱是其主值序列连续谱的离散抽样。

与上一讲中离散矩形脉冲的DTFT进行比较：可见，二者在时域上的关系是周期延拓，频域上的关系是离散抽样。如下图所示。

3. 四种傅里叶变换的关系

重点：一个域周期，对应另一个域离散。

说明

文章转自：个人图书馆/未注明作者，感谢原作者的辛勤付出，如涉及版权，请联系我删除！ http://www.360doc.com/content/19/0611/18/908538_841808371.shtml

五、离散傅里叶变换DFT

(一) 背景介绍

为什么需要定义一种“新”的变换?

计算机处理的两个基本条件：
第一，只能处理离散的数据(时域和频域都要离散)；
第二，要有限长。

DTFT，时域上离散，但频域是连续的；DFS，时域频域都是离散的，但同时又都是周期的，周期序列长度为无限长。但同时我们也注意到，周期序列实际上只有有限个序列值有意义，因而它的离散傅里叶级数也适用于有限长序列，这就得到有限长序列的离散傅里叶变换(DFT)。

所以，DFT并不是一种新的变换。它只是将DFS时域和频域上都取一个周期而已。DFT正反变换的定义式与DFS相同，只是加了一个取值范围的限定而已。换汤不换药。

(二) DFT的定义及物理含义

重点：DFT与DFS、DTFT的关系

1. 变换公式

$正变换：X(k)=\sum_{n=0}^{N-1} x(n) e^{-j \frac{2 \pi}{N} k n}, k=0,1, \ldots, N-1$ $逆变换：x(n)=\frac{1}{N} \sum_{k=0}^{N-1} X(k) e^{j \frac{2 \pi}{N} k n}, n=0,1, \ldots, N-1$

定义旋转因子$W_N=e^{-j\frac{2\pi}{N}}$，则上式转变为：

$正变换：X(k)=\sum_{n=0}^{N-1} x(n) W_{N}^{k n}, k=0,1, \ldots, N-1$ $逆变换：x(n)=\frac{1}{N} \sum_{k=0}^{N-1} X(k) W_{N}^{-k n}, n=0,1, \ldots, N-1$

总上，这对变换记为：

$x(n)\overset{N点\text{DFT}}{\longleftrightarrow}X(k)$

DFT不是序列x(n)的真正的频谱。x(n)的真正的频谱是DTFT，DFT只是对其真正频谱的一个周期上的离散抽样值。

(三) DFT的计算

重点：DFT的两种计算方法。

1. 方法一：利用DFT定义式

$X(k) = \sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}{N}kn},k=0,1,\cdots,N-1$

2. 方法二：先求DTFT再抽样

$X(k)=X(e^{j \omega})| _ {\omega=\frac{2 \pi}{N} k} R _ {N}(k)$

3. 方法三：DFT矩阵

设有一序列$x(n), n = 0, 1, \cdots, N-1$，求其$N$点DFT有：

$X(k) = \mathrm{DFT}[x(n)] = \sum_{n=0}^{N-1} x(n) \exp\left[{-\mathrm j}\dfrac{2\pi}{N}nk\right] = \sum_{n=0}^{N-1} x(n) W_N^{nk}$

其中，$k = 0,1, \cdots, N-1$。

令 $\boldsymbol x = \left[x(0), x(1),\cdots, x(N-1)\right]^{\mathrm T}$，则DFT计算可以看做向量$\boldsymbol x$和一个矩阵的乘积，我们一般讲这个矩阵称作DFT矩阵：

$\begin{aligned} \boldsymbol{F_1} &= \left[\begin{array}{ccccc} \mathrm e^{-\mathrm j \frac{2\pi}{N} 0 \cdot 0} &\mathrm e^{-\mathrm j \frac{2\pi}{N} 0 \cdot 1} & \mathrm e^{-\mathrm j \frac{2\pi}{N} 0 \cdot 2} & \cdots & \mathrm e^{-\mathrm j \frac{2\pi}{N} 0 \cdot (N-1)} \\ \mathrm e^{-\mathrm j \frac{2\pi}{N} 1 \cdot 0} &\mathrm e^{-\mathrm j \frac{2\pi}{N} 1 \cdot 1} & \mathrm e^{-\mathrm j \frac{2\pi}{N} 1 \cdot 2} & \cdots & \mathrm e^{-\mathrm j \frac{2\pi}{N} 1 \cdot (N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathrm e^{-\mathrm j \frac{2\pi}{N} (N-1) \cdot 0} &\mathrm e^{-\mathrm j \frac{2\pi}{N} (N-1) \cdot 1} & \mathrm e^{-\mathrm j \frac{2\pi}{N} (N-1) \cdot 2} & \cdots & \mathrm e^{-\mathrm j \frac{2\pi}{N} (N-1) \cdot (N-1)} \end{array}\right] \\ \\ &= \left[\begin{array}{ccccc} W_N^{0 0} & W_N^{0 1} & W_N^{0 2} & \cdots & W_N^{0 (N-1)} \\ W_N^{1 0} & W_N^{1 1} & W_N^{1 2} & \cdots & W_N^{1 (N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ W_N^{(N-1) 0} &W_N^{(N-1) 1} & W_N^{(N-1) 2} & \cdots & W_N^{(N-1)(N-1)} \end{array}\right] \end{aligned}$

此时，序列$x(n)$的DFT为：

$\boldsymbol X = \left[\begin{array}{c} X(0) \\ X(1) \\ \vdots \\ X(k) \\ \vdots \\ X(N-1) \end{array}\right] = \boldsymbol{F_1 x} = \boldsymbol{F_1}\left[\begin{array}{c} x(0) \\ x(1) \\ \vdots \\ x(n) \\ \vdots \\ x(N-1) \end{array}\right]$

或者，若将DFT矩阵定义为：

$\begin{aligned} \boldsymbol{F_2} = \boldsymbol{F_1}^{\mathrm H} &= \left[\begin{array}{ccccc} \mathrm e^{\mathrm j \frac{2\pi}{N} 0 \cdot 0} &\mathrm e^{\mathrm j \frac{2\pi}{N} 0 \cdot 1} & \mathrm e^{\mathrm j \frac{2\pi}{N} 0 \cdot 2} & \cdots & \mathrm e^{\mathrm j \frac{2\pi}{N} 0 \cdot (N-1)} \\ \mathrm e^{\mathrm j \frac{2\pi}{N} 1 \cdot 0} &\mathrm e^{\mathrm j \frac{2\pi}{N} 1 \cdot 1} & \mathrm e^{\mathrm j \frac{2\pi}{N} 1 \cdot 2} & \cdots & \mathrm e^{\mathrm j \frac{2\pi}{N} 1 \cdot (N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathrm e^{\mathrm j \frac{2\pi}{N} (N-1) \cdot 0} &\mathrm e^{\mathrm j \frac{2\pi}{N} (N-1) \cdot 1} & \mathrm e^{\mathrm j \frac{2\pi}{N} (N-1) \cdot 2} & \cdots & \mathrm e^{\mathrm j \frac{2\pi}{N} (N-1) \cdot (N-1)} \end{array}\right] \\ \\ &= \left[\begin{array}{ccccc} W_N^{-0 0} & W_N^{-0 1} & W_N^{-0 2} & \cdots & W_N^{-0 (N-1)} \\ W_N^{-1 0} & W_N^{-1 1} & W_N^{-1 2} & \cdots & W_N^{-1 (N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ W_N^{-(N-1) 0} &W_N^{-(N-1) 1} & W_N^{-(N-1) 2} & \cdots & W_N^{-(N-1)(N-1)} \end{array}\right] \end{aligned}$

此时，序列$x(n)$的DFT变换为：

$\boldsymbol X = \left[\begin{array}{c} X(0) \\ X(1) \\ \vdots \\ X(k) \\ \vdots \\ X(N-1) \end{array}\right] = \boldsymbol{F_2}^{\mathrm H} \boldsymbol x = \boldsymbol{F_2}^{\mathrm H}\left[\begin{array}{c} x(0) \\ x(1) \\ \vdots \\ x(n) \\ \vdots \\ x(N-1) \end{array}\right]$

容易发现，DFT矩阵具有如下性质：

① DFT矩阵是对称矩阵，即$\boldsymbol{F_1} = \boldsymbol{F_1}^{\mathrm T} \qquad \boldsymbol{F_2} = \boldsymbol{F_2}^{\mathrm T}$；
② $\boldsymbol{F_2} = \boldsymbol{F_1}^{\mathrm H}$；

参考链接1：第三章离散傅里叶变换（DFT）及其快速算法（FFT）
参考链接2：DFT和IDFT分析 - CSDN

4. 例题

下面的例题，分别用这两种方法进行求解。

求5点矩形窗函数$x(n) = R_5(n)$的5点DFT$X_1(k)$和10点DFT$X_2(x)$

对于本题来说，方法二特别直观，便于理解DFT与DTFT的关系。

可见，同一个序列的不同点数的DFT，得到的结果不同。DFT的点数$N$越大，$X(k)$越能反映连续频谱的形状。
当DFT的点数$N$>序列的长度$N_0$时，相当于在序列后面补上$N-N_0$个零，故称为”补零DFT“。

(四) DFT性质

1. 圆周移位(循环移位)

首先看，为什么要定义一种新的移位?这种新的移位为什么称为“圆周移位”？它与普通的移位有何异同？见下图。

可见，圆周移位与普通移位的区别在于：多一个周期延拓、再取主值区间的过程。定义式如下：

$圆周移位(循环移位)——x(n+m)_NR_N(n) = \tilde{x}(n+m)R_N(n)$

圆周移位的步骤：先周期延拓、再移位、最后取主值区间；或者先移位、再周期延拓、最后取主值区间。前两步谁先谁后都可以。见下图。

可见，这种新的移位，相当于先把序列放在一个圆上，然后转圈圈：左移，沿顺时针方向转(即上图中的情况);右移，沿逆时针方向转。
与其他傅里叶变换一样，DFT依然满足：时域移位，频域线性相移；频域移位，时域调制。需要注意的是，这里的移位，都指的是圆周(循环)移位。

2. 圆周卷积(循环卷积)

(1) 圆周卷积的定义
仿照上面我们对移位操作做的改进——增加了“周期延拓、取主值区间”的过程，将“移位”改造成“圆周移位”一样，我们也对之前很熟悉的卷积和公式做改进，就得到一种新的卷积”圆周卷积“(或称为循环卷积)的定义。

$x_1(n)$的长度为$N_1$，$x_2(n)$的长度为$N_2$，$N \geq \max(N_1,N_2)$，则$N$点圆周卷积为：

$x_1(n) \odot x_2(n) = \sum_{m=0}^{N-1}x_1(m)x_2(n-m)_NR_N(n)$

为了区分，大家在《信号与系统》中所学习的卷积和，称为”线性卷积“。
显然易见，圆周卷积的结果与$N$有关。
(2) DFT的圆周卷积定理
在圆周卷积的定义下，DFT依然满足”一个域卷积，另外一个域相乘“这一更古不变的关系。

时域圆周卷积定理：

$x_1(n) \odot x_2(n) \overset{N点\text{DFT}}{\longleftrightarrow}X_1(k)X_2(k)$

频域圆周卷积定理：

$x_1(n) x_2(n) \overset{N点\text{DFT}}{\longleftrightarrow}\frac{1}{N}X_1(k) \odot X_2(k)$

(3) 圆周卷积的计算
下面重点来看圆周卷积的计算。

(4) 圆周卷积与线性卷积的关系
上面的例题中，两个序列的线性卷积是多少呢?圆周卷积与之相同吗?为什么?
线性卷积，是直接将序列2反转、平移，而没有”周期延拓、再取主值“这两步。可以用图解法或者竖式法得到（信号与系统中已学，此处直接给出结果)，线性卷积的结果如下：

$两序列的5点线性卷积为:x _ 1(n) * x _ 2(n)=[1,3,6,9,12,9,5]$

两种卷积为何不同呢?又有何关系呢?
下图从公式上推导圆周卷积与线性卷积的关系

我们可以得到以下结论：

圆周卷积是线性卷积周期延拓后的主值序列：

$y_c(n) = \sum_{r=-\infty}^{\infty}y_l(n+rN)R_N(n)$

当$N\geq N_1+N_2-1$时，圆周卷积和线性卷积结果相同。

根据以上结论，我们可以根据线性卷积的结果，直接得到N点圆周卷积的结果。如下图：

3. 圆周共轭对称性

这里不讲证明(教材上都有)，重点讲怎么理解教材上让人眼花缭乱的公式。我们把“公式”翻译成“人话”。
首先说明，文中所说的N点长序列，都指的是自变量取值范围为$[0,N-1]$，除此之外的区间，序列值为0。
先看第一个。
(1) 共轭序列的DFT
时频域有这样一个基本对应关系——时域取共轭，对应频域自变量取负然后函数取共轭。具体到DFT呢？“自变量取负”也就是“反转”，而“DFT隐含着周期性”，所以这里的“反转”要加上“周期延拓，再取主值区间”，所以，公式及证明过程如下

时域取共轭，对应DFT是：先周期延拓，再反转，再取主值区间，最后取共扼。当然，第一步与第二步可以交换次序，取共轮可以放在任意步骤上。关键是理解这个操作用公式的三种描述方式(上图中画红线处)：

第①种：$X^*((-k))_N R_N(k)$，是最直观地展现上述过程的；
第②种：$X^*((N-k))_N R_N(k)$，可认为用其周期性(周期延拓得到的当然是以$N$为周期啦)，把$-k$换成$N-k$；
第③种：去掉了双括号，也去掉了$R_N(k)$，好像看不出“周期延拓”和“取主值区间”的操作了。大家会心存疑虑，这个等号成立吗?

我们用下图的例子来说明一下这个等号成立，为了画图的方便，我们用函数值为实数的情况，图中是以$n$为自变量，换作$k$当然也是一样的。

$x(N-n)$可以看作简写形式，优点在于形式简洁明了，缺点在于掩盖了周期延拓再取主值的过程。用这种简写形式，要注意一点，$N$点长序列$x(n)$，$n$的取值范围为$0<n<N-1$，也即本应该$x(N)=0$，但此处当$n=0$时，$x(N-n)=x(N)$，不能认为$x(N)=0$，而要认为$x(N)=x(0)$。也就是说，要把$x(n)$的这$N$个点，认为是周期序列的主值区间，那么$x(N)$就是下一个周期的第一个点，所以$x(N)=x(0)$。

用这种简写形式来描述这个性质，就是：时域取共轭，对应的DFT，相当于把序号$k$与序号$N-k$做一个互换，然后取共轭。

下面看这个性质的两个推论：

第一个推论：实序列的DFT是圆周共轭对称序列。

对于实序列$x(n)$，其DFT即$X(k)$满足： $X(k) = X^*(N-k)_NR_N(k) = X^*(N-k)$ 实序列的DFT $X(k)$是圆周共轭对称序列。即：实部圆周偶对称，虚部圆周奇对称；模圆周偶对称，相角圆周奇对称。

“圆周共轭对称”是个什么鬼？我们按照以下几步来解释一下：

第一步：从“偶对称、奇对称”到“共轭对称/共轭反对称”

偶对称/奇对称地球人都知道吧。共轭对称/反对称就不是地球人都知道了，大学生才知道。

对于实函数$x(n)$，如果$x(-n)=x(n)$，称之为偶对称，$x(-n)=-x(n)$称之为奇对称。

扩展到复函数$x(n)$，如果$x^*(-n)=x(n)$，称之为共轭对称，与之相对应的还有共轭反对称。
第二步：从“共轭对称”到“圆周共轭对称”

圆周共轭对称定义：对于$N$点长序列$x(n)$，若$x(n)=x((N-n))_NR_N(n)$，或者用简写形式：$x(n)=x(N-n)$，那么称之为“圆周共轭对称”。

可以理解为：把$x(n)$放在一个圆周的$N$个等分点上，或者说把横轴掰弯成一个圆($n=N-1$与原点重合)，则这$N$个序列值关于原点对称，或者说关于$\dfrac{N}{2}$也对称，如下图所示。

我们前面求解过的例题：5点矩形脉冲的DFT，如下图，也体现出圆周偶对称的特点。

第二个推论：实部/虚部与圆周共轭对称/反对称分量的关系

首先解释一下什么叫圆周共轭对称分量和圆周共轭反对称分量。需要经过以下几步循序渐进的理解。

第一步：实函数可以分解为偶分量+奇分量

$x(n) = x_e(n)+x_o(n)$

其中：偶分量 —— $x_e = \dfrac{x(n)+x(-n)}{2}$，奇分量 —— $x_o = \dfrac{x(n)-x(-n)}{2}$

第二步：从“偶分量/奇分量”到复函数可以分解为“共轭对称分量+共轭反对称分量”，把上式中的$x(-n)$改为$x^*(-n)$即可

$x(n) = x_e(n)+x_o(n)$

其中：共轭对称分量 —— $x_e(n) = \dfrac{x(n)+\overline{x(-n)}}{2}$，共轭反对称分量——$x_o = \dfrac{x(n)-\overline {x^*(-n)}}{2}$

以上两式，无论是对无限长序列，还是有限长序列，都是适用的。如果$x(n)$为$N$点长，并且$0<n<N-1$，那么$x_e(n)$和$x_o(n)$是$2N-1$点长，并且$-(N-1) \leq n \leq N-1$。

第三步：改造成适合DFT的。凡是涉及到自变量取负(也就是反转)的，都加上“周期延拓，再取主值区间”的操作。也就是把上式中的$x(-n)$改为$x((N-n))_NR_N(n)$，用简写形式表示就是$x(N-n)$。

因此，得到圆周共轭对称分量和圆周共轭反对称分量的定义：

$N$点长的序列$x(n)$，可以分解为圆周共轭对称分量+圆周共轭反对称分量

$x(n) = x_{ep}(n) + x_{op}(n)$

其中：圆周共轭对称分量——$x_{ep}(n) = \dfrac{x(n)+\overline{x(N-n)}}{2}$，圆周共轭反对称分量——$x_{op}(n) = \dfrac{x(n)-\overline{x(N-n)}}{2}$

注意，前提是$x(n)$为$N$点长序列，并且n的范围是$0≤n≤N-1$，圆周共轭对称/反对称分量的长度仍是$N$，$n$的范围也不变。而且如前所述，$n=0$时，$x(N-0)=x(N)=x(0)$。

上面，是以$x(n)$为例，同样，对于DFT $X(k)$，也可以定义圆周共扼对称/反对称分量，不再赘述。

解释完这些，我们的核心公式就出来啦（证明过程省略，直接看结论)。序列$x(n)$及其DFT的实部/虚部与圆周共轭对称/反对称分量之间的关系，见下图：

(此处省略若干公式)，翻译成人话(绕口令)就是：
序列实部的DFT是序列DFT的共轭对称分量
序列虚部×j的DFT是序列DFT的共轭反对称分量
序列共轭对称分量的DFT是序列DFT的实部
序列共轭反对称分量的DFT是序列DFT的虚部巧j

是不是像绕口令，但总比公式强多了。这—切，意义何在?

第一，从图形上可以淋漓尽致地体现DFT隐含的周期性。
第二，为DFT的简化运算提供了思路。

4. Parseval定理

有限长序列的能量：

$\sum\limits_{n=0}^{N-1}|x(n)|^{2}=\dfrac{1}{N} \sum\limits_{k=0}^{N-1}|X(k)|^{2}$

(五) 频域抽样

实际上，DFT就是频域抽样。包括三个问题，这三个问题环环相扣、层层推进。

1. DFT与DTFT、z变换的关系

先从公式上看三个变换的关系，再结合z平面的单位圆的概念，从图形上理解。如下图:

毫无疑问，DFT的自变量$k$为离散的，而DTFT的自变量$\omega$、以及z变换的自变量$z$都是连续变量。DFT是另外两种变换的离散采样值。因为这种采样是在频域，所以称为”频域采样“。

那么问题来啦：能否由离散频谱值$X(k)$得到$X(z)$和$X(e^{j\omega})$？

不管在那个域进行抽样，其数学本质都是用一些离散的数值代替原来连续变化的函数，或者说用一些离散的点代表原来连续的曲线。能不能代表取决于两个因素：一是这些离散的点的间隔，即抽样间隔；二是原来那条连续曲线的变化起伏程度，(频域抽样定理)。

2. 频域抽样定理

傅里叶分析方法的好处在于，建立起时域和频域的一种重要的对应关系：一个域离散抽样，另外一个域周期延拓。所以，研究时域抽样时，把问题对应到频域上去研究；那么现在研究频域抽样时，又要把问题对应到时域上去研究。毫无疑问，时域上会周期延拓。如下图：

既然是以$N$为周期延拓，条件自然而然就出来了：

对于$M$点长的序列$x(n)$，频域抽样不失真的条件是： 一个周期内的频域抽样点数$N \geq M$ ，此时有：

$x_N(n) = \sum_{r = -\infty}^{\infty}x(n+rN)R_N(n) = x(n)$

即满足该条件时，$N$个频率抽样值$X(k)$（即$x(n)$的$N$点DFT)能够完全代表$x(e^{j\omega})$及$X(z)$。

问题又来了，怎么表示？这就是第三个问题：频域的插值恢复。

3. 频域的插值恢复

与时域抽样的恢复完全相同的思路，用离散的样本值乘以一个插值函数，得到一个连续的函数，只不过这里的插值函数是关于$\omega$或$z$的函数。下面的任务就是找这个函数$\varphi(w)$或$\varphi(z)$。

z变换的形式更为简洁，因此首先解决由$X(k)$得到$X(z)$的问题。

以下推导过程的大致思路：把z变换定义式中的$x(n)$用IDFT的公式替换，然后交换求和次序，再利用旋转因子的性质，即可得到。如下图：

解决了由$X(k)$得到$X(z)$的问题，将$z$换成$e^{j\omega}$，自然就得到了$X(e^{j\omega})$。如下图：

把内插公式和内插函数总结如下图，这个内插函数的幅度部分的图形我们可以画出来，我们发现，它在一些固定的位置($\dfrac{2\pi}{N}$的整数倍处）是零，而$\dfrac{2\pi}{N}$恰好是频域抽样时的间隔，这是巧合吗？显然不是，这是必然的。

我们把内插公式展开来看，如下图所示。也就是说，把各个频域抽样值$X(k)$与做相应平移后的内插函数（平移$\dfrac{2\pi}N$的$k$倍）相乘，再相加，就得到连续的频谱函数$X(e^{j\omega})$。与第$k$个抽样值相乘的内插函数，在所有其他抽样点处刚好是零点，只有在第$k$个抽样点处的值不为零(值为1)。所以，重建后的这个连续函数，在每个抽样位置（也就是$\dfrac{2}{N}$的整数倍)上的值，就等于$X(k)$这一点的值，不需要任何其他抽样值参与;而在两个抽样点之间的值（没抽到的地方)，需要所有抽样值来参与共同构成。

这个问题的理解，与“时域抽样后信号的重建”问题是一样的。但有的人可能会说，时域抽样后信号的重建，我记得是通过理想低通滤波器来推导出重建的内插公式，这里怎么不是呢?
如果你能提出这样的疑问，要表扬，说明你“信号与系统”学的不错。其实，图7和图8中的内插公式，完全可以用“时域抽样信号重建连续时间信号的内插公式”相同的推导方法推导出来。感兴趣的同学可以试一下哦。

说明

文章转自：个人图书馆/未注明作者，感谢原作者的辛勤付出，如涉及版权，请联系我删除！

http://www.360doc.com/content/19/0611/18/908538_841808337.shtml

http://www.360doc.com/content/19/0611/18/908538_841808319.shtml

http://www.360doc.com/content/19/0611/18/908538_841808298.shtml

http://www.360doc.com/content/19/0611/18/908538_841808274.shtml

六、快速傅里叶变换FFT

(一) 改进DFT计算的方法

1. DFT计算量分析

对于$N$点长的序列$x(n)$其DFT变换为

$正变换：X(k) = \text{DFT}[x(n)] = \sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}{N}kn} = \sum_{n=0}^{N-1}x(n)W_N^{nk}$ $逆变换：x(n) = \text{IDFT}[X(k)] = \frac{1}{N}\sum_{k=0}^{N-1}X(k)e^{j\frac{2\pi}{N}kn} = \sum_{k=0}^{N-1}X(k)W_N^{-nk}$

DFT正反变换的运算方式和运算量是相同的，后面的分析以DFT正变换为例。

观察正变换和反变换的公式可知，二者的运算方式和运算量是完全相同的。下面的分析均以DFT正变换为例。(顺便说一句，大家要像熟悉自己的手机一样熟悉旋转因子，闭着眼睛都知道它)
观察DFT正变换的公式，容易看出：每计算一个点的数据，需要$N$次复数乘法、$N-1$次复数加法，所以$N$点DFT，需要$N$的平方次复数乘法，$N(N-1)$次复数加法。我们知道，DFT的点数，至少要取成序列的长度，当序列长度很长时，运算量巨大！如下表所示。

序列长度	复数乘法$N^2$	复数加法$N(N-1)$
2	4	2
8	64	56
64	4096	4032
1024	1048576	1047552
2048	4194304	4192256

以1024点为例，复数乘法的次数100万次之多。

1965年，库利(J.w.Cooley)和图基(J.W.Tukey)在《Mathmatics of Computation》上发表了《AnAlgorithm for the Machine Calculation of Complex Fourier Series》，提出一种DFT的快速算法，后人称为快速傅里叶变换(Fast Fourier Transform ——FFT)。

2. 改进DFT计算效率的基本途径

$N$点DFT，直接计算，需要$N$的平方次乘法；分成2个$\dfrac{N}{2}$点DFT分别计算，乘法的次数减少了一半；分成4个$\dfrac{N}{4}$点DFT，乘法的次数又减少了一半。如果能够继续下去，前景很让人向往。

为了能够一直分下去，我们限定$N$为2的整数次幂，即：$N=2^M$，称为基2FFT。由此可见，FFT的基本思想是：把长序列分成几个较短的序列。

但怎么分？不能随便分，基本原则：要保证这几个短序列的DFT组合起来后能够很方便地构成原来长序列的DFT。所以DFT快速算法要解决的两个核心问题是：怎么分？怎么合？

根据分与合的方式不同，有两种基2FFT算法，分别称为:

按时间抽取的FFT算法———Decimation-in-Time，简称DIT-FFT。
按频率抽取的FFT算法———Decimation-in-Frequency，简称DIF-FFT。

下面我们分别来归纳总结两种基2FFT算法。

(二) 两种基2FFT算法

1. 按时间抽取DIT-FFT算法

以第一次分解($N$点分为2个$\dfrac{N}{2}$点)为例来说明算法原理，首先解决怎么分的问题。

通俗地说，就是大家列队、报数（从0开始），报偶数的一组，奇数的一组。

然后解决怎么合的问题，我们略过看似艰苦卓绝实际很简单的推导过程，直接上结论：

公式不好看，有人画了一幅图，并且起了个好听的名字：蝶形运算符号。下面的动图演示了蝶形运算的过程:

以8点长序列为例，我们来看分解为2个4点长DFT，是如何通过蝶形运算合成8点DFT的：

经过第一次分解之后，总的运算量=两个N/2点DFT的运算+N/2个蝶形的运算。而每次蝶形运算，只需要1次乘法，2次加法。所以，总的乘法次数为：

$\frac{N^2}{2}+\frac{N}{2} \approx \frac{N^2}{2}$

总加法次数为：

$N(\frac{N}{2}-1)+N \approx \frac{N^2}{2}$

当$N$很大时，运算量减少了近一半。

这就给了我们信心，说明这种分解思路是可以有效减少运算量的。我们继续分解下去，经过M-1次分解，分解为N/2 个 2 点长序列。

而2点DFT也用蝶形运算来表示（因为计算机最擅长一致而重复的东西），就得到下面的流图：

2. 按频率抽取DIF-FFT算法

仍以第一次分解（$N$点分为2个$\dfrac N2$点）为例来说明算法原理。

以8点长序列为例，我们来看分解为2个4点长DFT，是如何通过蝶形运算合成8点DFT的：

注意到，输出的频率数据，序号是按照偶数一组、奇数一组的顺序排列的，所以这种算法称为：按频率抽取。我们继续分解下去，经过$M-1$次分解，分解为$N/2$ 个 2 点长序列，就得到下面的流图：

3. 运算量分析

通过前面的分析可见，两种基2FFT算法，运算量是一样的，N点DFT，就分解成了若干个蝶形的运算而已。多少个蝶形呢？序列长度$N=2^M$，共有 $M$级蝶形，每级$\dfrac N2$个蝶形，共$\dfrac{MN}{2}$个。而每个蝶形是1次复数乘法，2次复数加法。所以总的运算量为：

	FFT	DFT
复数乘法	$\dfrac{NM}{2} = \dfrac{N}{2}\log_2N$	$N^2$
复数加法	$NM = N\log_2N$	$N(N-1)$

频率作为自然界的一个基本物理量，是很多领域研究的重要内容。人们很早就认识到，用DFT的方法可以有效进行信号的频率分析。但是因为DFT算法运算量很大，在数字计算机发明以前，运算效率普遍很低的情况下，DFT也更多的是一种理论分析工具，很难被用于实际的信号处理。

FFT的出现，破解了这一历史性难题，极大地促进了数字信号处理这门学科的应用和发展。有人甚至以FFT算法提出的1965年作为数字信号处理这门学科的诞生之年。

4. 算法特点

在计算机看来，这两种算法是非常相像的。两者互为转置。

	DIF-FFT	DIF-FFT
同址运算	将蝶形运算的结果仍然保存在原输入量的存储单元中	将蝶形运算的结果仍然保存在原输入量的存储单元中
输入/输出顺序	输入倒位序,输出自然顺序	输入自然顺序，输出倒位序
蝶形运算	先乘旋转因子，后加减第$m$级节点间距离:$2^{m-1}$	先加减，后乘旋转因子第$m$级节点间距离:$2^{L-m}$

首先来看第一个特点：同址运算（又称同位运算或原位运算），每完成一个蝶形运算，输入的两个数据就没有用的，这就意味着，不需要再重新开辟新的存储单元来保存输出数据，计算结果仍保留在原输入数据占据的存储单元即可。

再来看第二个特点：输入/输出数据的顺序。这是两种算法的不同之处。以DIT-FFT为例来说明为什么会输入倒位序。

还是以8点长数据为例，输入数据的正常顺序是$x(0)、x(1)、x(2)……x(7)$，我们称之为自然顺序。按照序号的奇偶分为两组，第一组是$x(0)、x(2)、x(4)、x(6)$，第二组是$x(1)、x(3)、x(5)、x(7)$。每个新的组再重新排队报数，按奇偶分，第一组又分成两个组，分别是$x(0)、x(4)$和$x(2)、x(6)$，第二组分成两个组，分别是$x(1)、x(5)$和$x(3)、x(7)$。

也就是说，8点长序列的DIT-FFT，输入数据的顺序是：$x(0)、x(4)、x(2)、x(6)、x(1)、x(5)、x(3)、x(7)$。这个序号的顺序乍看杂乱无章，其实有规律性。0、1、2、3、4、5、6、7的顺序与0、4、2、6、1、3、5、7有何关系的呢？用二进制来写一目了然，看下面的动图：

倒位序，是将二进制数的最高有效位到最低有效位的位序进行颠倒排列而得到的二进制数。

DIT-FFT算法中，对时域序列按照序号的奇偶进行分解，造成输入序列的序号按照倒位序排列。

最后再说一说蝶形运算的规律。两种FFT算法，最终都是转换成了M列、每列N/2个、一共MN/2个蝶形运算。但二者蝶形运算的规律有差异。

第一个差异：基本蝶形不同。DIT是先乘旋转因子，再加或减；而DIF是先加或减，再乘旋转因子。
第二个差异：两种算法，蝴蝶翅膀的距离（即节点间的距离）和旋转因子的数目恰好相反。

仔细观察两种算法的流图，我们会发现，二者互为转置。

(三) 其他FFT算法简介

1. 混合基FFT

2. Chirp-z变换

实际应用中，有时只对信号的某一频段感兴趣，即只需要计算单位圆上某一段的频谱值，而不需要计算[0，Π]区间的所有频谱采样值。此时，可以用”Chirp-z变换“（CZT）。

适用场合：窄带信号的DFT。

3. Goertzel算法

在某些应用场合，只需计算很少几个频率点的频谱值。例如，双音多频信号（DTMF）的检测。此时可以采用卡泽尔（Goertzel）算法。

除此之外，傅里叶变换的快速算法还有很多种。不过应用最广泛的依然能是基2FFT算法，它是数字信号处理最经典算法之一，几乎各种主流的计算机编程语言都有现成的函数可以调用。不同型号的芯片，硬件开发商也都会给出优化后的FFT算法源代码，一般情况下直接调用就可以。

说明

文章转自：个人图书馆/未注明作者，感谢原作者的辛勤付出，如涉及版权，请联系我删除！ http://www.360doc.com/content/19/0611/18/908538_841808253.shtml

七、FFT算法的应用

说明

文章转自：个人图书馆/未注明作者，感谢原作者的辛勤付出，如涉及版权，请联系我删除！

http://www.360doc.com/content/19/0611/18/908538_841808231.shtml

http://www.360doc.com/content/19/0611/18/908538_841808103.shtml

八、Matlab中的fftshift函数

一般情况下，在fft()之间先对输入数据进行fftshift()处理是没有必要的。没有什么必然性。fftshift()这个函数的存在的目的并不是为了这个，单纯地就是上面所说的为了让频谱观测显得更自然一些而已。

这里先给出一个简单的解释（有时间再来补详细的解释）。这是从根本上来说是因为fft()处理的离散数据，进行的是离散傅里叶变换（DFT）。如果所要解析的数据本身就是一个离散周期性信号，那fft()给出的结果就反映了真实情况。但是现实应用中，所分析的信号并不是离散周期性信号，而我们所能做的又只能是调用fft()执行离散傅里叶变换进行分析，这种情况下fft()结果并没有完全反映真实情况，而是一种近似的或者“变形”的体现。这个时候我们就需要对fft()结果进行合理的解释间接地得到真实的情况，某种意义上fftshift()将频域数据前后半颠倒以使得频域数据显示出来的效果显得更自然也是属于这种情况。

比如(4)中给出的例子是对对称时域信号进行调用fft()处理，然后将其结果与傅里叶变换的理论解析结果进行对比。原始信号的时间区间是$[-\dfrac{T}{2}, \dfrac{T}{2}]$，这个数据采样后直接用fft()进行处理的话，由于fft()自然地认为数据是从t=0开始，所以可以理解为fft()的输入数据其实是将原始信号在时域上右移$\dfrac{T}{2}$所得到的信号，而fft()的变换结果则代表着这个时域上右移$\dfrac{T}{2}$所得到的信号的周期性拓展后的信号的离散时间傅里叶变换（DTFT）的一个周期内的数据。接下来从这个结果进行解析可反推得到原始信号的DFT结果。在fft()之前做fftshift()可以看作是一个小的trick, as an alternative to the post-analysis of the fft() result, without physical significance and necessity.

参考链接8.1：fftshift讲解 - 信号处理小王子的文章 - 知乎

参考链接8.2：关于MATLAB中fft，fftshift，fftshift(fft(fftshift(x)))，FFT要乘以dx 等问题理解 - 博客园

参考链接8.3：Matlab fftshift and ifftshift and some confusions - 笨牛慢耕 - CSDN

参考链接8.4：fftshift函数详解 - 博客园

九 Matlab实现部分简单信号处理算法

（一）傅里叶变换

$F(\omega) = \int_{-\infty}^{+\infty} f(t) \exp({-\mathrm j \omega t}) \text{ d}t$ $F(k) = \sum_{n=0}^{N-1} f(n) \exp\left(-\mathrm j \dfrac{2\pi}{N}nk\right)$

1. 实现方法1：matlab内置函数fft

%% matlab内置fft函数
clc; clear; close all;
% 采样频率
fs = 100;
% 采样点数
N = 128;
% 时间轴
n = 0:N-1;
t = n/fs;
% 正弦信号频率
f0 = 20;
% 正弦信号
x = sin(2*pi*f0*t);

% 作正弦信号的时域波形
figure;
subplot(1, 2, 1);
plot(t, x);
xlabel('t'); ylabel('y');
title('正弦信号时域波形');
grid on;

% 进行FFT变换并做频谱图
y = fft(x ,N);
% 进行对应的频率转换
f = (0:length(y)-1)'*fs/length(y);

% 做频谱图
subplot(1, 2, 2);
plot(f, abs(y));
xlabel('频率(Hz)'); ylabel('幅值');
title("正弦信号幅频谱图");
grid on;

输出结果

2. 实现方法2：自己编写函数

%% 自己编写dft函数  参考链接 —— https://www.jianshu.com/p/06c65b2a0408
clc; clear; close all;
% 采样频率
fs = 100;
% 采样点数
N = 128;
% 时间轴
n = 0:N-1;
t = n/fs;
% 正弦信号频率
f0 = 20;
% 正弦信号
data = sin(2*pi*f0*t);
% 时域信号转化为列向量
xk = data.';

% 作正弦信号的时域波形
figure;
subplot(1, 2, 1);
plot(t, xk);
xlabel('t'); ylabel('data');
title('正弦信号时域波形');
grid on;

% DFT需要的相关矩阵:
WN = exp(-1i*2*pi/N);
WN_nk = zeros(N, N)+WN;
% 辅助的E(WN_kn的幂,单独拿出来算)
E = zeros(N, N);   

% 构造DFT矩阵
for row = 0:N-1
    for cow = 0:N-1
        E(row+1, cow+1) = row*cow;
        WN_nk(row+1, cow+1) = WN_nk(row+1, cow+1)^E(row+1, cow+1);
    end
end

% 利用DFT矩阵计算频谱
Xk = WN_nk * xk;
% 进行对应的频率转换
f = (0:length(Xk)-1)'*fs/length(Xk);

% 做频谱图
subplot(1, 2, 2);
plot(f, abs(Xk));
xlabel('频率(Hz)'); ylabel('幅值');
title("正弦信号幅频谱图");
grid on;

输出结果

3. 实现方法3：轴采样实际数据，而不是离散化标号

%% 自己编写函数，轴采样实际数据
clc; clear; close all;
% 采样频率
fs = 100;
% 采样点数
N = 128;
% 时间轴
n = 0:N-1;
t = n/fs;
% 正弦信号频率
f0 = 20;
% 正弦信号
data = sin(2*pi*f0*t);
xk = data.';

% 作正弦信号的时域波形
figure;
subplot(1, 3, 1);
plot(t, xk);
xlabel('t'); ylabel('data');
title('正弦信号时域波形');
grid on;

% 频率轴设置方式1
CR_axis = (0:N-1)*fs/N;
% 傅里叶矩阵
F_matrix = exp(-1i*2*pi*CR_axis.'*t);
Xk = F_matrix*xk;

% 做频谱图
subplot(1, 3, 2);
plot(CR_axis, abs(Xk));
xlabel('频率(Hz)'); ylabel('幅值');
title("正弦信号幅频谱图(频率轴1)");
grid on;

% 频率轴设置方式2
CR_axis = (0:N-1)*fs/N-fs/2;
% 傅里叶矩阵
F_matrix = exp(-1i*2*pi*CR_axis.'*t);
Xk = F_matrix*xk;
% 做频谱图
subplot(1, 3, 3);
plot(CR_axis, abs(Xk));
xlabel('频率(Hz)'); ylabel('幅值');
title("正弦信号幅频谱图(频率轴2)");
grid on;

输出结果

（二）卷积、相关与类相关计算

1. 卷积

$y(t) = s(t) \otimes h(t) = \int_{-\infty}^{+\infty} s(u) h(t-u) \text{ d}u = \int_{-\infty}^{+\infty}s(t-u)h(u) \text{ d}u$

设滤波器的时间长度为$T$，以$t=0$为中心（非因果滤波器），则：

$y(t) = s(t) \otimes h(t) = \int_{t-T/2}^{t+T/2} s(u) h(t-u) \text{ d}u = \int_{-T/2}^{+T/2}s(t-u)h(u) \text{ d}u$

对于实际中的因果滤波器，定义域为$[0, T]$，则：

$y(t) = s(t) \otimes h(t) = \int_{t}^{t+T} s(u) h(t-u) \text{ d}u = \int_{0}^{T}s(t-u)h(u) \text{ d}u$

2. 自相关

$R(\tau) = \int_{-\infty}^{+\infty} s(t)s^*(t-\tau) \text{ d}\tau = \int_{-\infty}^{+\infty} s(t)s^*(t+\tau) \text{ d}\tau$ $R(m) = \sum_{n=-\infty}^{+\infty} s(n)s^*(n-m) = \sum_{n=-\infty}^{+\infty} s(n)s^*(n+m)$

实际中自相关函数$R[m]$，对于一个离散信号$x[n]$的定义为：

$R(m) = \sum_{n=-0}^{N-1} s(n)s^*(n-m) = \sum_{n=0}^{N-1} s(n)s^*(n+m)$

其中，$m$是滞后值，范围是$[−N+1,N−1]$， $N$是信号的长度。

实现方法1：matlab内置xcorr函数

%% 内置xcorr函数
clc; clear; close all;
% 计算并绘制信号的估计自相关。在零滞后时（此时 x 与自身完全匹配），出现最大峰值。
% 创建一个离散信号示例
x = [1+1i, 2-1i, 3+2i, 4-2i, 5+3i, 6-2i, 3+4i, 5+5i];
% 计算自相关函数
[r, lags] = xcorr(x);

% 绘制自相关函数
figure;
stem(lags, abs(r));
title('自相关函数');
xlabel('Lag'); ylabel('自相关');
grid on;

实现方法2：根据s(n)*s(n+m)

%% 自己编写，根据s(n)*s(n+m)
clc; clear; close all;
x = [1+1i, 2-1i, 3+2i, 4-2i, 5+3i, 6-2i, 3+4i, 5+5i];

R = my_xcorr_complex(x);
% 绘制自相关函数
lags = -length(x)+1:length(x)-1;
figure;
stem(lags, abs(R));
title('复数信号的自相关函数');
xlabel('滞后'); ylabel('自相关');
grid on;

function R = my_xcorr_complex(x)
    N = length(x);          % 信号的长度
    R = zeros(1, 2*N-1);    % 初始化自相关函数结果向量
    for k = -N+1:N-1        % 滞后值范围从 -N+1 到 N-1
        for n = 1:N         % 遍历信号的每个样本点
            if (n+k > 0) && (n+k <= N)  % 检查索引是否在有效范围内
                R(k+N) = R(k+N) + x(n) * conj(x(n+k));  % 累加乘积，包含共轭
            end
        end
    end
end

实现方法3：根据s(n)*s(n-m)

%% 自己编写，根据s(n)*s(n-m)
clc; clear; close all;
x = [1+1i, 2-1i, 3+2i, 4-2i, 5+3i, 6-2i, 3+4i, 5+5i];
R = my_xcorr_complex(x);

% 绘制自相关函数
lags = -length(x)+1:length(x)-1;
figure;
stem(lags, abs(R));
title('复数信号的自相关函数');
xlabel('滞后'); ylabel('自相关');
grid on;

function R = my_xcorr_complex(x)
    N = length(x);          % 信号的长度
    R = zeros(1, 2*N-1);    % 初始化自相关函数结果向量
    for m = -N+1:N-1        % 滞后值范围从 -N+1 到 N-1
        for n = 1:N         % 遍历信号的每个样本点
            if (n-m > 0) && (n-m <= N)  % 检查索引是否在有效范围内
                R(m+N) = R(m+N) + x(n) * conj(x(n-m));  % 累加乘积，包含共轭
            end
        end
    end
end

3. 类相关

这主要是是来自Wigner-Ville变换公式：

$W_x(t, \Omega) = \int_{-\infty}^{+\infty} x\left(t+\dfrac{\tau}{2}\right) x^*\left(t-\dfrac{\tau}{2}\right) \exp(-\mathrm j \Omega \tau) \text{ d}\tau$

离散化后：

$W_x(n,l) = 2 \sum_k r_x(n, k) \exp\left(-\mathrm j \dfrac{4\pi}{N}kl\right)$

注意，上面例子中，当$n=3$时候，我在计算的时候出现了下面可笑的结果：

$n=3时，r_x(3, k) = \{x_0x_6^*, x_1x_5^*, x_2x_4^*, x_3x_3^*, x_4x_2^*, x_5x_1^*, x_6x_0^* \}$

而实际上，这里是没有$x_6$这一项的，故$n=3$时只有5项。

（三）连续信号变换的离散化

这里简单用几句话总结一下目前的理解：
- ① 首先离散化积分变量(这里积分变量在信号处理中一般是时间变量或时延变量或频率变量)；
- ② 其次，当没有积分项后，其实剩下的连续变量主要是看该连续变量的取值范围，然后直接离散化就行。
两点注意：
- ì. 当考察到变量存在周期性时候，只需要在变量的一个周期内离散采样即可。
- ìì. 注意变量的取值范围，确定不了就试，一般就试$(0 \sim N-1)$或$(-N/2 \sim N/2-1)$范围。