标量/向量/矩阵的相关运算

矩阵乘法

\[\boldsymbol x^T \boldsymbol \Lambda \boldsymbol x = (x_1, \cdots, x_n) \left( {\begin{array}{*{20}{c}} {{\lambda _1}}&{}&{}\\ {}& \ddots &{}\\ {}&{}&{{\lambda _n}} \end{array}} \right) \left( {\begin{array}{*{20}{c}} {{x_1}}\\ \vdots \\ {{x_n}} \end{array}} \right) = \sum_{i=1}^n \lambda_i x_i^2 \]

知乎-矩阵求导术：上篇和下篇。本文基本上是这两篇文章内容的重新整理。
刘建平Pinard系列博客，这个博客主要用于查缺补漏
教材：《矩阵分析与应用》，作者张贤达
博客：矩阵求导总结
查询手册：The Matrix Cookbook

矩阵行列式求导以及矩阵的逆的求导 - CSDN

超级强大的在线矩阵求导工具：Matrix Calculus

2 信号处理中数学公式Matlab编程实现

2.1 Wigner分布的实现

2.1.1 理论介绍

首先，下面这部分内容来自书籍《现代信号处理教程（第二版）》，该书籍的相关信息请参考下面这两个链接：

现代信号处理教程（第二版）- 百度百科

现代信号处理教程（第二版）- 豆瓣

信号\(x(t)\)的自Wigner分布为：

\[W_x(t, \Omega) = \int_{-\infty}^{\infty} x(t+\frac{\tau}{2}) x^*(t-\frac{\tau}{2})\mathrm{e}^{-j\Omega \tau} \text{ d}\tau \]

信号处理算法要服务于实际工程，即最终的目的是要将它们应用于科研或工程的实际。因此，上面这种连续域无限长的理论公式必须要进行离散化和有限化才能在我们这个非连续有限的物理世界中实现，这时所遇到的问题就是信号的离散化及数据的有限长问题。

• 离散化

要离散化，就要对信号进行采样(把时间离散化)，对于信号\(x(t)\)，我们令采样间隔为\(T_s\)，即\(t = nT_s\)，则有：

\[W_x(t, \Omega) = \int_{-\infty}^{\infty} x(nT_s+\frac{\tau}{2}) x^*(nT_s-\frac{\tau}{2})\mathrm{e}^{-j\Omega \tau} \text{ d}\tau \]

进一步，要用计算机计算，连续的积分变量也要离散化，因此我们令\(\dfrac{\tau}{2} = kT_s\)，即\(\tau = 2kT_s\)，故而上式对\(\tau\)的积分变成对k的求和：

\[W_x(t, \Omega) = 2T_s \sum_{k=-\infty}^{\infty} x(nT_s+kT_s) x^*(nT_s-kT_s)\mathrm{e}^{-j2k\Omega T_s} \]

注意：上面这种描述我不知道是否准确，但是在写的过程中，还有如下一种思考：

信号\(x(t)\)的自Wigner分布为\(W_x(t, \Omega)\)，其中\(t\)和\(\Omega\)是连续变量，但是这2个连续变量并不影响在计算机上进行计算，当真正编程实现的时候，其实会发现原理论公式里面的连续积分变量才是真正的阻碍，我们可以随意给定\(t\)和\(\Omega\)的值，然后计算给定值下的具体值，比如令\(t = 1, \Omega = 0.5\)，则有：

\[W_x(1, 0.5) = \int_{-\infty}^{\infty} x(1+\frac{\tau}{2}) x^*(1-\frac{\tau}{2})\mathrm{e}^{-j0.5 \tau} \text{ d}\tau \]

这时候我们只要把\(\tau\)离散化，比如，令\(\dfrac{\tau}{2} = kT_s\)，即\(\tau = 2kT_s\)，则有：

\[W_x(1, 0.5) = 2T_s \sum_{k=-\infty}^{\infty} x(1+kT_s) x^*(1-kT_s)\mathrm{e}^{-j2k × 0.5 T_s} \]

假设这里暂不考虑有限性问题，那么就可以编程实现\((t = 1, \Omega = 0.5)\)这一点的Wigner分布值了。

虽然在\((t, \Omega)\)这个二维平面是随便连续地取值进行计算，但是我们总不能一个点一个点地取吧，要交给计算机有规律地计算这个二维平面上不同\((t, \Omega)\)处的值，我们需要将\((t, \Omega)\)二维屏幕也要有规律地进行离散化，于是就变成了全部离散化公式：

\[W_x(t, \Omega) = 2T_s \sum_{k=-\infty}^{\infty} x(nT_s+kT_s) x^*(nT_s-kT_s)\mathrm{e}^{-j2k\Omega T_s} \]

PS：还有一点要补充地是，为什么\(t\)离散化为\(t=nT_s\)，\(\tau\)离散化为\(\dfrac{\tau}{2} = kT_s\)呢，为什么不是\(\dfrac{\tau}{2} = kT_{\tau}\)呢？

额。。。其实在写这篇博客的时候感觉也不好说，或者自己还没有掌握本质，只能暂时解释为：要考虑\(\tau\)的物理意义——时延，当\(x(t)\)离散化后，只有在离散的\(t\)上可以取到值，如果\(\tau\)的离散化间隔\(T_{\tau} \neq T_s\)，那么就有可能把时延信号取在离散化后的信号\(x(nT_s)\)的间隔之间(连续值是不用考虑这种情况的，毕竟连续值在任何地方都有值)，导致无法取到值，无法进行计算。

离散化之后，信号其实就变成了一系列离散的点了，比如\(x(n) = [0.5, 1, 5, -1, -2, \cdots]\)，计算机在计算的时候只会见到这些数值，这里我就是想表明计算机在带入数值计算时是不知道\(T_s\)等于多少的，因此我们可以人为地把\(T_s\)归一化为1，并考虑到相对离散信号的频率\(\omega = \Omega T_s\)，则上式变为：

\[W_x(n, \omega) = 2 \sum_{k=-\infty}^{\infty} x(n+k) x^*(n-k)\mathrm{e}^{-j2k\omega} \]

注意1：这里\(n=0，1，2，\cdots\)是离散化了，但是\(\omega = \Omega T_s\)还是没有离散化，因为虽然\(T_s\)离散化了，但是\(\Omega\)是连续的，一个连续的值\(\Omega\)乘以一个离散固定值\(T_s\)仍然是连续值。

注意2：只是在对计算机计算时可以人为地把\(T_s\)归一化为1，并不是说我们要忘记\(T_s\)的值，当上述通过上式计算出结果后，我们要把对应\(t\)的轴乘以\(T_s\)，对应\(\Omega\)的轴除以\(T_s\)

在继续讲解有限性之前，还有一个点需要插入，但是这里我就先不讲了，因为我自己也还无法弄懂，简单粘贴一下书中截图，后续如果有机会再补充：

• 有限性

现在还剩下两个问题需要解决，① 一是频率\(\omega\)的离散化；② 二是确定式中\(k\)的求和范围。

现在令：

\[r_x(n, k) = x(n+k)x^*(n-k) \]

并假定\(x(n)\)的长度为\(N\)，即\(n=0,1, \cdots, N-1\)，再分析一下\(r(n,k)\)的取值情况。

对于信号离散化之后的下标问题——即采用正负对称下标还是从0开始的下标，这里暂时我也还无法确定，有可能两种方法都行，也有可能需要按照情况讨论，不管怎么样，按照书中的讲解，这里就以上面的假定为例，即：\(n=0,1, \cdots, N-1\)

这里直接通过具体计算来展示，比如我们令\(N = 6\)，则有：

下标索引\(n\)	\(r_x(n,k)\)的值	\(k\)的取值范围
\(n=0\)	\(r_x(0, k) = [x_0x^*_0]\)	\(k=0\)
\(n=1\)	\(r_x(1, k) = [x_0x^*_2, x_1x^*_1, x_2x^*_0]\)	\(k=-1, 0, 1\)
\(n=2\)	\(r_x(2, k) = [x_0x^*_4, x_1x^*_3, x_2x^*_2, x_3x^*_1, x_4x^*_0]\)	\(k=-2, -1, 0, 1, 2\)
\(n=3\)	\(r_x(3, k) = [x_1x^*_5, x_2x^*_4, x_3x^*_3, x_4x^*_2, x_5x^*_1]\)	\(k=-2, -1, 0, 1, 2\)
\(n=4\)	\(r_x(4, k) = [x_3x^*_5, x_4x^*_4, x_5x^*_3]\)	\(k=-1, 0, 1\)
\(n=5\)	\(r_x(5, k) = [x_5x^*_5]\)	\(k=0\)

根据上表可以看出，由于下标采用的是\(n=0,1, \cdots, N-1\)，因此限制了不同\(n\)取值下\(k\)的取值范围，总之：一条准则，\(k\)的取值范围不能使信号\(x_?\)的下标出现超出\(n=0,1, \cdots, N-1\)的范围。

另外，对于计算机而言，数据长度不一样的时候需要进行0填充（PS：这句话肯定还需进一步修改，这里暂时这样记录），将\(r(n,k)\)都扩充成\(N\)点序列，即在其后补零，那么上表即：

方法1：在序列后面补零：

下标索引\(n\)	\(r_x(n,k)\)的值后填充0	\(k\)的取值范围
\(n=0\)	\(r_x(0, k) = [x_0x^*_0, 0, 0, 0, 0, 0]\)	\(k=0\)
\(n=1\)	\(r_x(1, k) = [x_0x^*_2, x_1x^*_1, x_2x^*_0, 0, 0, 0]\)	\(k=-1, 0, 1\)
\(n=2\)	\(r_x(2, k) = [x_0x^*_4, x_1x^*_3, x_2x^*_2, x_3x^*_1, x_4x^*_0, 0]\)	\(k=-2, -1, 0, 1, 2\)
\(n=3\)	\(r_x(3, k) = [x_1x^*_5, x_2x^*_4, x_3x^*_3, x_4x^*_2, x_5x^*_1, 0]\)	\(k=-2, -1, 0, 1, 2\)
\(n=4\)	\(r_x(4, k) = [x_3x^*_5, x_4x^*_4, x_5x^*_3, 0, 0, 0]\)	\(k=-1, 0, 1\)
\(n=5\)	\(r_x(5, k) = [x_5x^*_5, 0, 0, 0, 0, 0]\)	\(k=0\)

方法2：在序列两侧补零：

下标索引\(n\)	\(r_x(n,k)\)的值	\(k\)的取值范围
\(n=0\)	\(r_x(0, k) = [0, 0, x_0x^*_0, 0, 0, 0]\)	\(k=0\)
\(n=1\)	\(r_x(1, k) = [0, x_0x^*_2, x_1x^*_1, x_2x^*_0, 0, 0]\)	\(k=-1, 0, 1\)
\(n=2\)	\(r_x(2, k) = [x_0x^*_4, x_1x^*_3, x_2x^*_2, x_3x^*_1, x_4x^*_0, 0]\)	\(k=-2, -1, 0, 1, 2\)
\(n=3\)	\(r_x(3, k) = [x_1x^*_5, x_2x^*_4, x_3x^*_3, x_4x^*_2, x_5x^*_1, 0]\)	\(k=-2, -1, 0, 1, 2\)
\(n=4\)	\(r_x(4, k) = [0, x_3x^*_5, x_4x^*_4, x_5x^*_3, 0, 0]\)	\(k=-1, 0, 1\)
\(n=5\)	\(r_x(5, k) = [0, 0, x_5x^*_5, 0, 0, 0]\)	\(k=0\)

参考书籍中的原文说的是选取方法1：在序列后面补零，但是，我觉得可能使用方法2：在序列两侧补零更合适一点，毕竟上面\(k\)变量的取值是正负对称的，当然不管序列后补零还是两侧补零，应该都会涉及fftshit问题，这个我可能穿插在代码部分，这里就先不展开了，这样表达式就变成了：

\[W_x(n, \omega) = 2 \sum_{k=-\lfloor (N+1)/2 \rfloor}^{\lfloor (N+1)/2 \rfloor} r_x(n, k)\mathrm{e}^{-j2k\omega} \]

最后，看一下\(\omega\)的离散化，这里我呢也不深究，只说一下目前的简要理解，后期如果有其他想法再补充，首先我们根据复指数函数的性质可知：

\[W_x(n, \omega + \pi) = 2 \sum_{k}r_x(n, k)\mathrm{e}^{-j2k(\omega + \pi)} = W_x(n, \omega) \]

也即，对于\(\omega\)变量来说，是按照\(\pi\)周期性变化的，因此，我们只需考察\(0 \sim \pi\)内的数值即可，其他部分都是周期重复的，所以现在可把\(\omega=0 \sim \pi\)区间内进行采样离散化，至于采样间隔多大，我也不知道，但是一般和离散信号离散点数一样即可吧，即把\(\omega=0 \sim \pi\)区间划分\(N\)份：

\[\omega = 0, \dfrac{1}{N}\pi, \dfrac{2}{N}\pi, \cdots, \dfrac{N-1}{N}\pi \]

所以，最终公式的完全离散化表达为：

\[W_x(n, l) = 2 \sum_{k=-\lfloor (N+1)/2 \rfloor}^{\lfloor (N+1)/2 \rfloor} r_x(n, k)\mathrm{e}^{-j2k\frac{\pi}{N}l}, \quad l = 0,1, \cdots, N-1 \]

2.1.2 代码实现

既然完成了公式的离散化，下面就可以编写程序实现了，一个示例程序如下所示。

• 首先产生一个线性调频信号

clc; clear; close all;
% 调频率
kr = 4;
% 信号持续时间4s
T = 4;
% 带宽
B = kr*T;
% 采样频率
fs = 4*B;
% 采样间隔
Ts = 1/fs;
% 采样点数
N = T*fs;
% 时间轴 
t = -N/2:N/2-1;
% 解析信号 时间轴 
x = 2*exp(1i*kr*pi*(t*Ts).^2);
X = fftshift(fft(fftshift(x)));

figure;
subplot(2, 2, 1);
plot(t*Ts, real(x));
xlabel('t'); ylabel('x(t)'); legend('x(t)');

subplot(2, 2, 2);
plot(t*fs/N, abs(X));
xlabel('f'); ylabel('|X(f)|'); legend('X(w)');

• 然后，调用一个函数来计算\(r_x(n,l)\)矩阵，也可以自己写，然后通过一个例子来展示经过该函数后的输出

function R = rx(x, N)
# PS：这段代码来自下面的原创力文档链接
R = zeros(N, N);
for n = 0:N-1
    M = min(n, N-1-n);
    for k = 0:M
        R(n+1, k+1) = x(n+k+1)*conj( x(n-k+1) );
    end
    for k = N-1:-1:N-M
        R(n+1, k+1) = conj( R(n+1, N-k+1) );
    end
end
% 上面设计的循环逻辑我个人是人为有点问题，因此按照我自己的理解：下面3行单独把R的第1列挪到最后一列
% PS1：这个问题本质就是fftshit和ifftshift的问题，我现在还没搞懂，因此下面3行可要可不要
% PS2：这个逻辑应该好修改，但是我懒得改了，毕竟不是目前的主要任务
RM = R;
R(:, 1:N-1) = RM(:, 2:N);
R(:, N) = RM(:, 1);
end

示例测试该函数：

clc; clear; close all;
% 偶数情况
x = ones(1, 8);
Rx = rx(x, length(x));
figure(1);
heatmap(Rx)

% 奇数情况
y = ones(1, 9);
Ry = rx(y, length(y));
figure(2);
heatmap(Ry)

输出结果如下：

偶数情况	奇数情况

🚀 有3个点需要注意：
① 第1：偶数情况下，该函数在信号\(x(t)\)对应的\(R_x\)矩阵中间添加了一个全0列矩阵，不这样添加的话，对于偶数信号的\(R_x\)矩阵就不是方阵了。
② 第2：上面的图是保留函数最后3行的结果，自己也可以尝试去掉最后3行运行查看结果，差异就是最后一列在前在后的区别。
③ 第3：通过与我自己用visio绘制的图形可知，该函数的输出结果是采用方式2，并且经过fftshift之后的结果，fftshift就相当于把\(k=[-\lfloor (N+1)/2 \rfloor, \lfloor (N+1)/2 \rfloor]\)变成\(k=[1, 2, ..., N-2, 0]\)

• 进行WVD变换

TF = zeros(N, N);
for n = 0:N-1
    temp = fftshift(fft( R(n+1, :) ));
    TF(n+1, :) = temp;
end
% 频率轴除以2的原因请看下面的注释
fnew = t*fs/2/N;
tnew = t*Ts;
[F, T] = meshgrid(fnew, tnew);
subplot(2, 2, 3);
mesh(F, T, abs(TF));
xlabel('频率');
ylabel('时间');
subplot(2, 2, 4);
contour(F, T, abs(TF));

可以看到，上面的代码利用了fft函数进行快速WVD变换，可以这样做是因为，考察DFT函数：

\[X(l) = \text{DFT}[(x(k)] = \sum_{k=0}^{N-1} x(k) \mathrm{e}^{-jk \frac{2\pi}{N} l}, \quad l = 0, 1, 2, ..., N-1 \]

我们把DFT函数与离散化的WVD表达式对比，就很明显可以得到，把\(n\)视为固定的数，\(r_x(n, k)\)就可视为DFT函数中的\(x(k)\)，所以WVD这不就是变成了一个DFT了吗，最终的输出结果如下：

PS：解释上面代码中频率轴为什么除以2：但是要注意，虽然数值形式上完全一样，但是可以观察到数字2的位置不同，在两个式子中的代表的意义也不同，也就是说，我们只是单纯借助fft来得到WVD的结果，但是在傅里叶变换中，默认频率轴的周期是\(2\pi\)，而WVD中的“频率轴”的周期是\(\pi\)，因此，用fft计算出的结果要除以2。
又或者可以这样认为，直接用fft来计算结果，相当于把\(2\omega\)看成一个整体了（\(2\omega\)这个整体的频率不就是\(2\pi\)了嘛）。

参考链接1：MATLAB辅助现代工程数字信号处理 第二版 作者 李益华 第6-10章_ 第8章.ppt - 原创力文档

2.2 离散傅里叶DFT的实现

2.2.1 连续傅里叶变换的离散化

有了上面的基础，下面我自己想试着写一下DFT的实现。

连续傅里叶变换的公式为：

\[F(\Omega) = \int_{-\infty}^{\infty}f(t)\mathrm{e}^{-\text{j}\Omega t} \text{ d}t \]

首先，对积分变量t进行离散化：\(t = nT_s\)，则\(dt = (n+1)T_s-nT_s=T_s\)，故有：

\[F(\Omega) = \sum_{n=-\infty}^{\infty}f(nT_s)\mathrm{e}^{-\text{j}\Omega nT_s} * T_s \]

对于计算机来说，离散化后的\(f(nT_s)\)和\(f(n)\)没有区别，因为带入计算的都是一串序列数字而已。此外，与WVD类似，考虑到相对离散信号的频率\(\omega = \Omega * T_s\)(Ps：这里的\(T_s\)可不能归一化为1，因为\(T_s\)是有实际的物理意义的，例如\(\sin(t)\)每隔\(T_s=\frac{\pi}{2}\)采样和\(\sin(2t)\)每隔\(T_s=\frac{\pi}{4}\)采样得到的序列是完全一样的，如果不考虑\(T_s\)，那么就无法区分这两个信号了)，注意此时\(\omega\)仍然是连续变量：

\[F(\omega) = T_s * \sum_{n=-\infty}^{\infty}f(n)\mathrm{e}^{-\text{j}n\omega} \]

下面，就要考虑有限化了，这里很简单，信号\(f(t)\)以\(T_s\)采样间隔得到离散化的采样序列，但是实际中不可能无限采样，故将采样序列的索引限制为：\(n = 0, 1, \cdots, N-1\)，或者说不是限制采样序列索引长度，而是我们实际就采样了\(N\)个点：

\[F(\omega) = T_s * \sum_{n=0}^{N-1}f(n)\mathrm{e}^{-\text{j}n\omega} \]

这就是离散时间傅里叶变换DTFT了。理论上，到此，我们就可以计算任意一个\(\omega\)处的\(F(\omega)\)值了，例如计算\(\omega = 0.401\)处的值，则只需把\(0.401\)代入公式然后交给计算机进行计算就行了，但是为了能有规律的计算不同\(\omega\)上的点并绘制出一条线便于观察，我们可以对\(\omega\)也进行离散化，也即：DTFT --> DFT。

观察到，\(F(\omega)\)是周期性的，即：\(F(\omega)\)=\(F(\omega+2\pi)\)，因此，只需要在\([0, 2\pi]\)直接对其离散化就行了，离散多少个点呢，实际上可以随意，即：

\[0, \frac{2\pi}{M}, 2*\frac{2\pi}{M}, \cdots, (M-1)*\frac{2\pi}{M} \]

上面的\(M\)就表示\(\omega\)在一个周期内离散化的点数，一般情况下(比如Matlab官方的fft函数)，都是默认\(M=N\)的，最终得到离散化后的公式：

\[F(k) = T_s * \sum_{n=0}^{N-1}f(n)\mathrm{e}^{-\text{j}n k \frac{2\pi}{M}} \]

至此，算是彻底完成了连续傅里叶变换FT的离散化。下面可以编写程序了。

2.2.2 根据自己离散化的公式编写程序

clc; clear; close all;
% ====================================================
% 随机构建一个信号x(n)
% ====================================================
x = [1 0.6 0.2 -0.2 -0.6 -1 -0.6 -0.2 0.2 0.6 ...
    1 0.6 0.2 -0.2 -0.6 -1 -0.6 -0.2 0.2 0.6 ...
    1 0.6 0.2 -0.2 -0.6 -1 -0.6 -0.2 0.2 0.6 ...
    1 0.6 0.2 -0.2 -0.6 -1 -0.6 -0.2 0.2 0.6 ...
    1 0.6 0.2 -0.2 -0.6 -1 -0.6 -0.2 0.2 0.6 ...
    1 0.6 0.2 -0.2 -0.6 -1 -0.6 -0.2 0.2 0.6 ...
    1 0.6 0.2 -0.2 -0.6 -1 -0.6 -0.2 0.2 0.6];
x = [x, x];

% ====================================================
% 使用fft进行傅里叶变换
% ====================================================
X = fft(x);
figure;
plot(abs(X), 'LineWidth', 2.5);
hold on;

% ====================================================
% 使用自己离散化公式进行DFT
% ====================================================
N = length(x);  % 频域的采样点可以自己调整
n = 0:length(x)-1;
K = 1;
myX = zeros(1, N);
for k = 0:1/N:1-1/N
    temp = exp(-1i*2*pi*k*n).';
    myX(K) = x*temp;
    K = K+1;
end
plot(abs(myX), '.-y')
legend('官方fft','自写DFT')

结果如下：

2.2.4 补充1：波形分辨率与频率分辨率

一般而言，n点信号的离散傅里叶变换 (DFT) 的变换结果（频域）也是n个数据点。但在实际应用中，对实际信号作FFT时，常常涉及到变换前数据需要补零 (Zero padding) 的问题。一些论坛里，曾看到某些专业人士从信息论的角度分析认为：“Zero padding没有增加时域信号的有效信息，因此，不会改变DFT/FFT的分辨率”。

• 什么是补零

简单来说，补零 (Zero Padding) 就是对变换前的时域或空域信号的尾部添加若干个0，以增加数据长度。如下图所示，为含有\(1.00MHz\)和\(1.05MHz\)两个频率成分合成的正弦波实信号。

左图中信号长度为1000个样点，采样频率为\(f_s=100MHz\)时，信号的实际时长则为\(10us\)。在其尾部添加1000个0，即数据增加到了2000个点（时长为\(20us\)），则变为右所示的波形。我们所关心的补零会不会影响计算输出的频率分辨率呢？

• 两种分辨率的辨析

这里涉及到两种意义下的分辨率问题，一种叫“波形频率分辨率”，一种叫“FFT分辨率”。波形频率分辨率是指可以被分辨的2个频率的最小间隔 (Spacing)；而FFT 分辨率则是频谱中的数据点数 (The number of points in the spectrum)，它是与做FFT的点数直接相关的。

因此，波形频率分辨率可定义为：

\[\Delta R_w = \frac{1}{T} \]

其中，\(T\)表示实际有效的采样信号时长，不包括人为补零的时长。

同样，FFT分辨率可以定义为：

\[\Delta R_f = \frac{f_s}{N_{fft}} \]

其中，\(f_s\)为采样频率 (the sampling frequency)，\(N_{fft}\)为FFT的点数。

值得注意的是，可能有很好的FFT分辨率，但却不一定能够很好的把2个频率成分简单的分开。同样，可能有很高的波形分辨率，但波形的能量峰值会通过整个频谱而分散开（这是因为FFT的频率泄漏现象）。可以通过下面的例子解释：

例1：信号时长为1秒，采样率\(1000 Hz\)，FFT点数1000

波形频率分辨率：\(\Delta R_w = \frac{1}{T} = 1/1 = 1 Hz\)（理论上可分辨间隔\(≥1 Hz\)的两个频率）

FFT分辨率：\(\Delta R_f = \frac{f_s}{N_{fft}} = 1000/1000 = 1Hz\)（频谱中每\(1 Hz\)显示一个频点）

例2：信号时长为0.5秒，采样率1000 Hz，FFT点数2000（补零）

波形频率分辨率：\(\Delta R_w = \frac{1}{T} = 1/0.5 = 2 Hz\)（实际只能分辨间隔\(≥2 Hz\)的频率）

FFT分辨率：\(\Delta R_f = \frac{f_s}{N_{fft}} = 1000/2000 = 0.5Hz\)（频谱每\(0.5 Hz\)显示一个频点，但实际无法分辨\(0.5 Hz\)间隔的频率！）

(1) 补零（Zero-Padding）的陷阱

现象：通过补零增加FFT点数（如例2），频谱显示的频点间隔更小（FFT分辨率提高）。

误区：误以为补零能提高物理分辨率。实际上，补零仅增加频谱插值点数，使频谱曲线更平滑，但无法提升真实的频率分辨能力。

验证：若两个频率间隔小于\(\dfrac{1}{T}\)，即使补零到\(N_{fft}→∞\)，频谱主瓣仍会重叠，无法分辨两个信号。

(2) 如何真正提高分辨率？

提升波形频率分辨率：必须增加信号的实际时长T（例如采集更长时间的数据）。

提升FFT分辨率：仅需增加\(N_{fft}\)或降低\(f_s\) ，但这只是优化显示效果，并没有很特别的意义。

关于问题：我不明白波形频率分辨率为什么由信号实际时长T决定？为什么只有增加T才能提高真实分辨率？这个结论是如何推导出来的？

1. 直观解释：时间越长，分辨越清晰
波形频率分辨率\(\dfrac{1}{T}\)表示能够分辨的两个频率之间的最小间隔。为什么它和信号的实际时长有关呢？我们可以打个比方：

想象你在听两个音调，比如钢琴上的两个键，它们的频率很接近（比如 440 Hz 和 441 Hz）。如果只听一瞬间（比如 0.1 秒），你可能听不出差别，因为耳朵需要时间来感知频率的周期性。但如果听的时间延长到 2 秒，你会更容易分辨出这两个音调的微小差异。信号处理也是类似的道理。

2. 数学推导：从傅里叶变换的本质看

对于一个连续时间信号\(x(t)\)，其傅里叶变换为：\(X(f) = \displaystyle\int_{-\infty}^{\infty}x(t)\mathrm{e}^{-\text{j}2\pi f t} \text{ d}t\)，在实际应用中，我们通常处理的是有限长度的信号（即截断信号），假设信号持续时间为\(T\)。此时，信号可以表示为：

\[{x_T}(t) = x(t) \cdot w(t) \qquad w(t)= \left\{ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {1}&{0 \le t \le T} \end{array}}\\ {\begin{array}{*{20}{c}} 0&{{\rm{otherwise}}} \end{array}} \end{array}} \right. \]

对这个有限长度信号进行傅里叶变换，得到的频谱为：\(X_T(f) = \displaystyle\int_{0}^{T}x(t)\mathrm{e}^{-\text{j}2\pi f t} \text{ d}t\)，也就相当于加了窗函数\(w(t)\)。

窗函数的频域影响：矩形窗\(w(t)\)的傅里叶变换是\(\text{sinc}\)函数：\(W(f) = T \cdot \text{sinc}(fT) = T \cdot \dfrac{\sin(\pi fT)}{\pi f}\)，\(\text{sinc}\)函数的主瓣宽度决定了能够分辨的最小频率间隔。

\(\text{sinc}\)函数的主瓣宽度（从零点到零点的距离）是\(\dfrac{2}{T}\)，而主瓣的半宽（中心到第一个零点的距离）是\(\dfrac{1}{T}\)。根据卷积定理，时域乘积对应频域卷积：

\[X_T(f) = X(f)*W(f) \]

这意味着原始信号的频谱\(X(f)\)会被\(\text{sinc}\)函数“模糊”。\(\text{sinc}\)函数的主瓣宽度决定了频率分辨能力。下面是一个例子：

假设信号\(x(t)\)是两个频率接近的正弦波之和：\(x(t) = \cos(2\pi f_1t) + \cos(2\pi f_2t)\)，它的频谱 \(X(f)\) 是两个冲激函数：\(𝛿(𝑓−𝑓_1)+𝛿(𝑓−𝑓_2)\)。经过有限时长\(T\)的截断后，频谱变成：

\[X_T(f) = X(f)*W(f) = [𝛿(𝑓−𝑓_1)+𝛿(𝑓−𝑓_2)]*T \cdot \text{sinc}(fT) \]

结果是两个\(\text{sinc}\)函数，分别以\(f_1\)和\(f_2\)为中心，宽度约为\(\dfrac{1}{T}\)。① 如果\(|f_1 - f_2|<\dfrac{1}{T}\) ：两个\(\text{sinc}\)函数重叠严重，无法区分两个峰值。② 如果\(|f_1 - f_2|>\dfrac{1}{T}\)：两个\(\text{sinc}\)函数的主瓣分开，可以分辨出两个频率。

因此，频率分辨率定义为\(\dfrac{1}{T}\)，它是能够分辨两个频率的最小间隔。

关于问题：补零等价于频域sinc函数内插

进行zero padding只是增加了数据的长度，而不是原信号的长度。就好比本来信号是一个周期的余弦信号，如果又给它补了9个周期长度的0，那么信号并不是10个周期的余弦信号，而是一个周期的余弦加一串0，补的0并没有带来新的信息。其实zero padding等价于频域的sinc函数内插，而这个sinc函数的形状（主瓣宽度）是由补0前的信号长度决定的，补0的作用只是细化了这个sinc函数，并没有改变其主瓣宽度。而频率分辨率的含义是两个频率不同的信号在频率上可分，也就要求它们不能落到一个sinc函数的主瓣上。所以，如果待分析的两个信号频率接近，而时域长度又较短，那么在频域上它们就落在一个sinc主瓣内了，补再多的0也是无济于事的。

1. 问题的核心

补零（zero-padding）是在时域信号末尾添加零，使信号长度从\(𝑁\)增加到\(𝑀\)。我们需要证明：
【1】补零后的频谱是对原始频谱\(X(\omega)\)的插值。
【2】插值函数是\(\text{sinc}\)函数。
这需要从时域和频域的变换关系入手，结合采样定理来推导。

2. 证明：补零与sinc插值的关系

（1）原始信号及其频谱

假设有一个离散信号\(𝑥[𝑛]\)，长度为\(𝑁\)，采样频率为\(𝑓_s\)。其离散时间傅里叶变换（DTFT）为：

\[X(\omega) = \sum_{n=0}^{N-1}x[n]\mathrm{e}^{-j\omega n} \]

\(X(\omega)\)是连续的频域函数，\(\omega\)是归一化频率（范围\([−𝜋,𝜋]\)）。离散傅里叶变换（DFT）是对\(X(\omega)\)在\(𝑁\)个等间隔点采样：

\[X[k] = X(\omega) |_{\omega=\frac{2\pi k}{n}} = \sum_{n=0}^{N-1}x[n]\mathrm{e}^{-j2\pi kn/N}, \qquad k = 0,1,2,\cdots, N-1 \]

频率间隔为\(Δ𝑓=\frac{𝑓_𝑠}/{𝑁}\)。

（2）补零后的信号

补零后，信号变为\(𝑥_{\text{pad}} [𝑛]\)，长度为\(𝑀(M>N)\)：

\[{x_{{\rm{pad}}}}[n] = \left\{ {\begin{array}{*{20}{c}} {x[n],}&{0 \le n < N}\\ {0,}&{N \le n < M} \end{array}} \right. \]

其DFT为：

\[X_{\rm{pad}}[k] = \sum_{k=0}^{M-1}x_{\rm{pad}}[n]\mathrm{e}^{-j2\pi kn/M} \]

由于\(x_{\rm{pad}}[n]\)在\(N \leq n < M\)时为零，则上面的补零DFT可化简为：

\[X_{\rm{pad}}[k] = \sum_{k=0}^{N-1}x[n]\mathrm{e}^{-j2\pi kn/M} \]

频率间隔变为\(\dfrac{f_s}{M} < \dfrac{f_s}{N}\)，频谱点更密集。我们需要证明\(X_{\rm{pad}}[k]\)是\(𝑋(𝜔)\)的插值结果，且插值函数是\(\text{sinc}\)。

证明步骤1：补零后的频谱表达式

\[X_{\rm{pad}}[k] = \sum_{k=0}^{N-1}x[n]\mathrm{e}^{-j2\pi kn/M} \]

将其视为\(𝑋(𝜔)\)在特定频率点的值：

\[X_{\rm{pad}}[k] = X(\omega) |_{\omega=\frac{2\pi k}{𝑀}} \]

这表明\(X_{\rm{pad}}[k]\)是\(X(\omega)\)在\(𝑀\)个等间隔点（\(𝜔=\dfrac{2𝜋𝑘}{𝑀}\)）上的采样。我们需要找到\(X_{\rm{pad}}[k]\)与原始\(𝑋[𝑘]\)的关系。

证明步骤2：频域插值的定义

插值的目标是从原始DFT\(𝑋[𝑘]\)（\(𝑁\)个点）重建\(𝑋(𝜔)\)的连续形式，然后在更密集的点（\(𝑀\)个点）上采样。理想插值应满足：

🌊 ① 在原始点\(𝑘=0,1,…,𝑁−1\)上，插值结果等于\(𝑋[𝑘]\)。
🚎 ② 在新点上，插值结果是对\(𝑋(ω)\)的合理估计。

证明步骤3：从原始DFT重建

已知信号\(x[n]\)的DFT变换为\(X[k]\)，根据DFT的逆变换，原始信号\(x[n]\)可以通过\(𝑋[𝑘]\)重建：

\[x[n] = \dfrac{1}{N} \sum_{k=0}^{N-1}X[k]\mathrm{e}^{-j2\pi kn/N} \]

将其代入\(𝑋(𝜔)\)的定义：

\[X(\omega) = \sum_{n=0}^{N-1}x[n]\mathrm{e}^{-j\omega n} = \sum_{n=0}^{N-1} \left(\dfrac{1}{N} \sum_{k=0}^{N-1}X[k]\mathrm{e}^{-j2\pi kn/N}\right)\mathrm{e}^{-j\omega n} \]

交换求和顺序：

\[X(\omega) = \dfrac{1}{N} \sum_{k=0}^{N-1}X[k] \sum_{n=0}^{N-1} \mathrm{e}^{-j(\omega - 2\pi k/N) n} \]

证明步骤4：计算内层求和

内层是一个几何级数：

\[\sum_{n=0}^{N-1}e^{-j(\omega-\frac{2\pi k}{N})n}=\frac{1-e^{-j(\omega-\frac{2\pi k}{N})N}}{1-e^{-j(\omega-\frac{2\pi k}{N})}} =e^{-j(\omega-\frac{2\pi k}{N})\frac{N-1}{2}}\cdot\frac{\sin\left(\frac{(\omega-\frac{2\pi k}{N})N}{2}\right)}{\sin\left(\frac{\omega-\frac{2\pi k}{N}}{2}\right)} \]

这正是\(\mathrm{sinc}\)函数的形式（归一化后）。因此：

\[X(\omega) = \sum_{k=0}^{N-1}X[k] \dfrac{1}{N} \mathrm{e}^{-j(\omega-\frac{2\pi k}{N})\frac{N-1}{2}}\cdot\frac{\sin\left(\frac{(\omega-\frac{2\pi k}{N})N}{2}\right)}{\sin\left(\frac{\omega-\frac{2\pi k}{N}}{2}\right)} \]

这就是用\(𝑋[𝑘]\)重建\(𝑋(𝜔)\)的插值公式，插值核是周期化的\(\mathrm{sinc}\)函数（也叫Dirichlet核）。

证明步骤5：补零后的频谱与插值

将\(𝜔=\dfrac{2𝜋𝑘}{𝑀}\)代入：

\[X_{\mathrm{pad}}[k]=X(e^{j\omega})|_{\omega=\frac{2\pi k}{M}}=\sum_{l=0}^{N-1}X[l]\cdot\frac{1}{N}\cdot\frac{\sin\left(\frac{(\frac{2\pi k}{M}-\frac{2\pi l}{N})N}{2}\right)}{\sin\left(\frac{\frac{2\pi k}{M}-\frac{2\pi l}{N}}{2}\right)}\cdot e^{-j(\frac{2\pi k}{M}-\frac{2\pi l}{N})\frac{N-1}{2}} \]

整理后：

\[X_{\mathrm{pad}}[k]=\sum_{l=0}^{N-1}X[l]\cdot\mathrm{sinc}\left(\frac{\pi N(k/M-l/N)}{\pi}\right)\cdot e^{-j\phi} \]

其中，\(\text{sinc}(x) = \dfrac{\sin(\pi x)}{\pi x}\)，相位项\(\mathrm{𝑒}^{−𝑗𝜙} = \mathrm{𝑒}^{-j(\frac{2\pi k}{M}-\frac{2\pi l}{N})\frac{N-1}{2}} = \mathrm{𝑒}^{-j \pi (N- 1)(\frac{k}{M}-\frac{l}{N})}\)是由于窗口对齐产生的。

证明步骤6：为什么是\(\mathrm{sinc}\)函数？

\(\qquad \cdot\)采样定理：由香农采样定理，从离散样本重建连续信号的理想插值函数是\(\mathrm{sinc}\)。若频谱\(𝑋(𝜔)\)是带限的(带宽小于\(\dfrac{𝑓_𝑠}{2}\))，\(\mathrm{sinc}\)插值能完美重建。

\(\qquad \cdot\)补零的本质：补零后，频谱\(𝑋_{\text{pad}}[k]\)是\(𝑋(𝜔)\)在更密集点的采样，而这种插值过程正好对应\(\mathrm{sinc}\)函数的卷积效应。

理想情况下，\(\mathrm{sinc}\)插值公式为：

\[X_{\mathrm{pad}}[k]=\sum_{l=0}^{N-1}X[l]\cdot\mathrm{sinc}\left(\frac{\pi N(k/M-l/N)}{\pi}\right)\cdot e^{-j\phi} \]

补零的有限长度版本是对其的近似。

3. 物理意义：为什么sinc函数？

🚄 \(\mathrm{sinc}\)函数的来源：\(\mathrm{sinc}\)函数是理想低通滤波器（矩形窗）的傅里叶变换。在频域中，插值需要保持带限特性，\(\mathrm{sinc}\)函数是唯一满足这一条件的函数。

🚖 补零的频域效应：时域补零相当于频域与一个更宽的矩形窗卷积，其结果是用\(\mathrm{sinc}\)函数插值原始频谱。

步骤5中公式的推导整理过程：

“相位项\(𝑒^{−𝑗𝜙}\)是由于窗口对齐产生的”的解释：

参考链接：Sinc滤波器的解释 - 博客园

参考链接：傅里叶变换的波形分辨率与频率分辨率 - 微信公众号

3 实信号频谱对称，复信号频谱非对称

4 DFT/FFT 等效为窄带滤波器组

4.1 相干积累FFT实现方式：DFT/FFT 等效为窄带滤波器组

信号\(x[n]\)的离散傅里叶变换DFT为：

\[\begin{aligned} X[k] &= \sum_{n=0}^{N-1} x[n] \mathrm{e}^{-j\frac{2\pi}{N}kn} \\ &= \sum_{n=0}^{N-1} x[n] \mathrm{e}^{-j\frac{2\pi}{N}kn} \times 1 \\ &= \sum_{n=0}^{N-1} x[n] \mathrm{e}^{-j\frac{2\pi}{N}kn} \times \mathrm{e}^{j\frac{2\pi}{N}kN} \\ &= \sum_{n=0}^{N-1} x[n] \mathrm{e}^{-j\frac{2\pi}{N}k(N-n)} \end{aligned} \]

设信号\(x[n]\)与分别与冲激响应为\(h[n]=\mathrm{e}^{-j\frac{2\pi}{N}kn}(\text{ }k=0, 1, \cdots, N-1)\)的滤波器在时域卷积，则有：

\[y(m) = \sum_{n=0}^{N-1} x[n]h[m-n] \]

此卷积在\(m=N\)时刻的输出为：

\[\begin{aligned} y(N) &= \sum_{n=0}^{N-1} x[n]h[N-n] \\ &= \sum_{n=0}^{N-1} x[n]\mathrm{e}^{-j\frac{2\pi}{N}k(N-n)} \end{aligned} \]

可以观察到，两者是完全一样的，这样就说明了信号\(x[n]\)的DFT在\(k_0\)处的值\(X[k_0]\)等于信号\(x[n]\)与冲激响应为\(h[n]=\mathrm{e}^{-j\frac{2\pi}{N}k_0n}\)的滤波器在时域卷积。同时，我们又知道时域卷积等于频域相乘，而\(h[n]=\mathrm{e}^{-j\frac{2\pi}{N}k_0n}\)的数字频谱是一个\(\text{sinc}\)函数，因此，对于\(k=0, 1, 2, \cdots, N-1\)构成一个窄带滤波器组，中心数字角频率分别是：\(\dfrac{2\pi}{N} k\)，每个窄带滤波器对等于其中心频率的单音信号的增益为：\(N\)。

4.2 用窄带滤波器组原理解释相干积累增益

https://zhuanlan.zhihu.com/p/609960189

5 非相干积累和相干积累

5.1 非相干积累

6 相干积累-DFT/FFT-窄带滤波器组-DBF

https://zhuanlan.zhihu.com/p/609960189

7 深入浅出雷达脉冲压缩技术

https://blog.csdn.net/helloQQp/article/details/119965517

8 信号的相关性分析

首先，从离散傅立叶变换公式：

\[S(k) = \text{DFT}[(s(n)] = \sum_{n=0}^{N-1} s(n) \mathrm{e}^{-jn \frac{2\pi}{N} k}, \quad k = 0, 1, 2, ..., N-1 \]

可以看出，每个DFT输出 只是输入信号与余弦/正弦波之间逐项乘积的总和，这实际上是相关性的计算。

根据定义，相关性是两个信号之间相似性的度量。在日常生活中，我们通过与已知事物的相关性来识别事物。相关性在科学、体育、经济、商业、市场营销、犯罪学和心理学等各个领域中都发挥着重要作用，甚至可以用一本书来详细探讨这个主题。

对于夏洛克·福尔摩斯的所有推理，他在观察之后的下一步总是进行相关性分析。例如，他通过将一些观察与脑中的模板进行对比(相关性分析)，准确描述了詹姆斯·莫蒂默博士的狗：
他的牙齿痕迹非常明显。在我看来，狗的下颌在这些痕迹之间的空间显示得太宽了，对于一只梗犬来说太宽，而对于一只獒犬来说又不够宽。它可能是什么——是的，天哪，它是一只卷毛猎犬。

8.1 实信号的相关

8.1.1 理论分析

两个信号之间相关性的目的是测量这两个信号的相似程度。从数学角度来看，两个连续信号\(s(t)\)和\(h(t)\)的互相关函数为：

\[\text{Corr}(\tau) = R_{\text{sh}}(\tau) = \int_{-\infty}^{+\infty} s(t)h(t+\tau)\text{ d}t \]

两个离散信号\(s[n]\)和\(h[n]\)之间的相关性定义为：

\[\text{Corr}[m] = R_{\text{sh}}[m] = \sum_{m=-\infty}^{+\infty} s[n]h[n+m] \]

上式表示的是相关运算，是两数字序列对应项相乘再相加的运算。式中\(m\)表示位移量，\(m>0\)表示\(h(n)\)序列左移，\(m<0\)表示\(h(n)\)序列右移，不同的\(m\)得到不同的\(\text{Corr}[m]\)值。如\(\text{Corr}[1], \text{Corr}[0], \text{Corr}[-1]\)。\(\text{Corr}[m]>0\)，表示有相同成分存在，\(\text{Corr}[m]<0\)表示有相反成分存在，\(\text{Corr}[m]=0\)表示两序列正交或者相互独立。

简单来说：\(\text{Corr}[m]\)表示在延迟\(m\)时，信号\(s[n]\)和\(h[n]\)的相似程度。其幅值受信号幅值影响，因此一般用其归一化形式：

\[\text{Corr}[m] = R_{sh}[m] = \dfrac{\sum_{m=-\infty}^{+\infty} s[n]h[n+m]}{\displaystyle\sqrt{\sum_{m=-\infty}^{+\infty} s[n]^2 \displaystyle\sum_{m=-\infty}^{+\infty}h[n]^2}} \]

归一化后的互相关函数\(\text{Corr}[m] = R_{\text{sh}}[m]\)的取值范围在\([−1,1]\)之间，值越接近 1 表示两个信号在对应延迟下越相似。

直观上，逐点乘积的求和计算了每个时间移位处两个信号的相似程度。例如，如果\(s[n]\)与\(s[n]\)进行相关性计算，那么在时间移位为0时，两个信号完全相同，相关性输出将是样本平方和，即信号\(s[n]\)的能量。

8.1.2 举例计算

设有两个信号：\(s[n] = [2, -1, 1]，h[n] = [-1, 1, 2]\)，两个离散信号计算相关性的图示流程如下图所示：

8.1.3 代码及实验

在信号处理中，互相关法是一种重要的分析手段。互相关函数主要用于衡量两个信号在时间上的相似性。它通过滑动一个信号相对于另一个信号，并计算重叠部分的相似程度，来识别信号之间的时间对齐关系。
如下两个离散信号\(x[n]\)和\(y[n]\)，其互相关函数图如下：

• 从时域图可以看出，两个信号长得比较像，但怎么表达两个函数长得像呢？ 我们直观可以感觉到，信号\(y[n]\)要比\(x[n]\)延迟一些，幅值小一些，但二者形状上是比较像的。
• 通过作互相关，我们可以计算出二者在不同时间差异下的相似程度，则相似程度最大处的横坐标，即代表了二者信号的相位差。
• 我们首先对信号进行取样，取样长度\(N=1000\)，互相关函数幅值最大处\(\tau=961\)，因此\(961 − ( 1000 − 1 ) = − 38\)，即信号\(x[n]\)要比\(y[n]\)超前38个采样点的时间。

% MATLAB代码演示互相关的意义
close all;
clear all;
% 参数设置
Fs = 1000;           % 采样频率 (Hz)
t = 0:1/Fs:1-1/Fs;   % 时间向量 (1秒)
f = 5;               % 信号频率 (Hz)
delay = 40;          % y(n)相对于x(n)的延迟样本数
similarity = 0.8;    % 相似程度系数 (0到1之间)
% 生成信号x(n)
x = sin(2*pi*f*t);
% 生成信号y(n)，通过相似程度系数调整，并添加延迟
y = similarity * sin(2*pi*f*(t - delay/Fs));
% 绘制x(n)和y(n)的时域图
figure('Name','互相关函数演示','Position',[400,400,1200,800]);
subplot(2,1,1);
plot(t, x, 'b', 'LineWidth', 1.5);
hold on;
plot(t, y, 'r--', 'LineWidth', 1.5);
title('信号x(n)和y(n)的时域图');
xlabel('时间 (s)');
ylabel('幅度');
legend('x(n)', 'y(n)');
grid on;
% 计算互相关函数，使用 'biased' 选项
N = length(x);
R_xy = xcorr(x, y, 'biased');  % 使用有偏互相关
% 找到互相关函数的最大值及其对应的位置
[max_corr, max_idx] = max(R_xy);
lag_max = max_idx - (N - 1);
% 绘制互相关函数
subplot(2,1,2);
stem(0:length(R_xy)-1, R_xy, 'k', 'LineWidth', 1.5);
title(['互相关函数 R_{xy}(\tau)  (最大值 = ', num2str(max_corr), ')']);
xlabel('延迟 \tau (样本)');
ylabel('互相关幅度');
grid on;
% 标记最大值位置
hold on;
plot(lag_max, max_corr, 'ro', 'MarkerFaceColor', 'r');
text(lag_max, max_corr, ['(\tau = ', num2str(lag_max), ')'], 'Color', 'r', 'VerticalAlignment', 'bottom');
hold off;
% 显示延迟信息
fprintf('信号y(n)相对于x(n)的延迟为 %d 个样本。\n', lag_max);

8.2 复数信号的相关性

复数相关性中一个信号需要取共轭,两个连续复信号\(s(t)\)和\(h(t)\)的互相关函数为：

\[\text{Corr}(\tau) = R_{\text{sh}}(\tau) = \int_{-\infty}^{+\infty} s(t)h^*(t+\tau)\text{ d}t \]

两个离散复信号\(s[n]\)和\(h[n]\)之间的相关性定义为：

\[\text{Corr}[m] = R_{\text{sh}}[m] = \sum_{m=-\infty}^{+\infty} s[n]h^*[n+m] \]

参考链接1：信号处理：互相关函数 - CSDN
参考链接2：:通信专业中的相关性到底在讲什么？

9 功率谱和功率谱密度和相关函数

9.1 功率谱、功率谱密度

任意实信号\(x(t)\)，及其扩展函数：绝对值\(|x(t)|\)和平方\(x^2(t)\)。如果是个电压或电流信号，我们可以得到能量Energy和功率Power的相关定义。如任意时间段的总能量，瞬时功率，平均功率等。

假设电压信号\(x(t)\)作用到负载为\(1\Omega\)的电阻上，则瞬时功率\(P(t)\)、在\([t_1, t_2]\)时间段上消耗的能量\(E\)和平均功率\(P_{\text{av}}\)为：

\[\begin{aligned} & P(t) = x^2(t) \\ & E = \int_{t_1}^{t_2} x^2(t) \text{ d}t = \sum_{n=0}^{N-1} |x[n]|^2\\ & P_{\text{av}} = \dfrac{E}{t_2 - t_1} = \dfrac{1}{t_2 - t_1}\int_{t_1}^{t_2} x^2(t) \text{ d}t \end{aligned} \]

\(x(t)\)存在傅里叶变换的充分条件为，\(x(t)\)满足绝对积分（absolute integratable）且有有限个间断点，转换为公式：

\[\int_{-\infty}^{\infty} |x(t)| \text{ d}t < \infty \]

在存在傅里叶变换的条件下，根据帕萨瓦尔（Parseval）定理，我们知道同样有着能量守恒的关系:

\[E = \int_{-\infty}^{\infty} |x(t)|^2 \text{ d}t = \dfrac{1}{2\pi} \int_{-\infty}^{\infty} |X(j\omega)|^2 \text{ d}\omega = \int_{-\infty}^{\infty} |X(j 2\pi f)|^2 \text{ d}f \]

其中，\(X(j\omega)\)是\(x(t)\)的傅里叶变换，同时离散形式为：

\[E = \sum_{n = 0}^{N-1} |x[n]|^2 = \dfrac{1}{N}\sum_{k=0}^{N-1} |X[k]|^2 \]

由上式我们可了解到，信号在频域的能量求解方法与时域类似，也为模的平方后求和。但是为使时域与频域的总能量相等，需将频域的能量除以信号的长度\(N\)，也可称为周期\(T\)。此处将\(N\)称为\(T\)的原因是在离散傅里叶变换中，会将输入信号当成一个无限循环、周期为T的的周期信号。根据功率的定义，我们在离散帕塞瓦尔的公式上除以\(N\)，即可获得信号在\(T\)内的功率：

\[P = \dfrac{1}{N}\sum_{n = 0}^{N-1} |x[n]|^2 = \dfrac{1}{N^2}\sum_{k=0}^{N-1} |X[k]|^2 \]

功率谱PS与功率谱密度PSD最主要的区别就在频率分量与单位频率这两个概念上，为解释这两个概念大家可看以下图像：

假设使用24Hz的采样频率对某信号采样0.5s，一共采样了12个点，频域分辨率$\text{d}f$的计算公式如下： $$ \text{d}f = \dfrac{f_s}{N} $$ 式中，$f_s$为采样频率，$N$为信号长度。由上可得频域分辨率为$\text{d}f=2Hz$，$[2,4,6,8,10,12]$即为各频率分量。该图像的横轴为频率，纵轴为值的大小，其作用为展示功率谱PS与功率谱密度PSD值之间的关联性，并无实际物理意义。

以各频率分量为中心，频域分辨率\(\text{d}f\)为宽度，将频率划分为\([1,3],[3,5],[5,7],[7,9],[9,11],[11,13]\)六个区间。以频率分量为横轴，各区间内的功率大小为纵轴绘制的图像即为功率谱图，如图中蓝线所示。所有频率分量的功率谱之和即为信号的总功率：

\[P_{\text{all}} = \text{PS}_2 + \text{PS}_4 + \cdots + \text{PS}_{12} \]

若使用各频率分量处的功率值除以区间宽度\(\text{d}f\)，即可得功率在该区间内的密度，该值即为功率谱密度，如图中红点所示。由此我们可得以下公式：

\[P_{\text{all}} = (\text{PSD}_2 + \text{PSD}_4 + \cdots + \text{PSD}_{12}) * \text{d}f \]

图中黄线为真实的信号功率谱密度，可理解为区间宽度df趋近于无穷小时的功率谱密度。

通过以上介绍，我们可得出频率成分k处的功率谱PS与功率谱密度的计算公式，如下所示：

\[\text{PS}(k) = \dfrac{1}{N^2} |X[k]|^2 \qquad \quad \text{PSD}(k) = \dfrac{1}{N^2 \text{d}f} |X[k]|^2 \]

9.2 matlab计算功率谱的四种方法

• 功率谱计算代码

计算功率谱的四种方法，前三种方法采用归一化频率，第四种方法采用原始频率。

clear;clc;close all;
% 生成一个示例时间序列（例如正弦波加上一些随机噪声）
t = 0:0.01:100-0.1;  % 时间向量
% x = sin(t) + 0.5 * randn(size(t));  % 正弦波加噪声
x = sin(20*pi*t);

%% 使用 pwelch 计算功率谱
% 输入信号 x，窗口大小为 256，重叠为 128，采样频率为 1 Hz
[pxx, f] = pwelch(x, 256, 128, [], 1);

% 绘制功率谱
figure;
plot(f, 10*log10(pxx));  % 使用对数尺度绘制功率谱（dB）
title('功率谱 (使用 Welch 方法)');
xlabel('频率 (Hz)');
ylabel('功率 (dB)');

%% 使用 periodogram 计算功率谱
[pxx, f] = periodogram(x, [], [], 1);

% 绘制功率谱
figure;
plot(f, 10*log10(pxx));  % 使用对数尺度绘制功率谱（dB）
title('功率谱 (使用 Periodogram 方法)');
xlabel('频率 (Hz)');
ylabel('功率 (dB)');

%% 直接使用fft
% 计算信号的 FFT
X = fft(x);

% 计算功率谱
pxx = abs(X).^2 / length(x);  % 计算功率谱
f = (0:length(x)-1) / length(x);  % 频率轴

% 绘制功率谱
figure;
plot(f, 10*log10(pxx));  % 使用对数尺度绘制功率谱（dB）
title('功率谱 (使用 FFT 方法)');
xlabel('频率 (Hz)');
ylabel('功率 (dB)');

%% 使用自相关函数的方法
% 计算自相关函数
[R, lags] = xcorr(x, 'biased');
% 计算功率谱
N = length(x);  % 信号长度
Fs = 100;  % 采样频率

% 获取自相关函数的一半（因为自相关函数是对称的）
R_half = R(N:end);

% 计算自相关函数的傅里叶变换，得到功率谱
S = abs(fft(R_half));
f = (0:length(S)-1) * (Fs / length(S));  % 频率轴

% 可选：将功率谱转换为dB
S_dB = 10 * log10(S);
figure;
plot(f, S_dB);  % 显示dB功率谱
xlabel('Frequency (Hz)');
ylabel('Power (dB)');
title('Power Spectrum using Autocorrelation');

• 功率谱密度计算代码
  • Matlab实例：频谱、功率谱和功率谱密度计算 - 微信公众号
  • 离散信号功率谱密度、功率谱、频谱、能量谱的matlab简单计算 - 知乎
  • MATLAB中功率谱密度计算pwelch函数使用详解 - CSDN

9.3 自相关函数和信号的功率谱密度的关系

在物理实际中，我们碰到的很多噪声信号可近似为宽平稳随机过程（wide-sense stationary process）。有信号自相关函数的定义:

\[R_x(t_1, t_2) = \mathbb{E}[x(t_1)x(t_2)] \]

其中，\(\mathbb{E}\)为求期望。如果自相关函数\(R_x(t_1, t_2)\)不仅依赖于时间差\(τ=t_2​−t_1​\)，还与具体的时间\(t_1\)或\(t_2\)相关，则过程是非平稳的。常见的非平稳过程包括：维纳过程（Wiener Process，布朗运动）‌；如果\(x(t)\)是宽平稳随机过程，那么\(R_x(t_1, t_2)\)仅和时间差\(\tau = t_2 - t_1\)有关，可简写为：\(R_x(\tau)\)。

平稳过程‌：分析简单，可通过傅里叶变换研究功率谱密度（自相关函数的傅里叶变换）。
‌非平稳过程‌：需使用时变分析工具（如短时傅里叶变换、小波变换），或对信号进行平稳化预处理（如去趋势、差分）。

维纳-辛钦定理（wiener khinchin theorem）揭示了自相关函数和信号的功率谱密度的关系。信号\(x(t)\)的自相关函数和功率谱密度函数互为傅里叶变换对。并且当\(τ=0\)时，功率谱密度\(S_x(f)\)的频域积分为\(R_x(0)\)，可以看做\(x(t)\)的平均功率。当然还有更特殊的case，如果信号\(x(t)\)均值为零，那么\(R_x(0)\)还可以看做\(x(t)\)的方差（variance），这是一个比较特殊的结论。

\[\left\{ \begin{array}{l} {S_x}(\omega ) = \displaystyle\int_{ - \infty }^\infty {{R_x}(\tau ){{\rm{e}}^{ - {\rm{j}}\omega \tau }}{\rm{d}}\tau } \\ {R_x}(\tau ) = \dfrac{1}{{2{\rm{\pi }}}}\displaystyle\int_{ - \infty }^\infty {{S_x}(\omega ){{\rm{e}}^{{\rm{j}}\omega \tau }}{\rm{d}}\omega } \end{array} \right. \quad \text{OR} \quad \left\{ \begin{array}{l} {S_x}(f) = \displaystyle\int_{ - \infty }^\infty {{R_x}(\tau ){{\rm{e}}^{ - {\rm{j2\pi }}f\tau }}{\rm{d}}\tau } \\ {R_x}(\tau ) = \displaystyle\int_{ - \infty }^\infty {{S_x}(f){{\rm{e}}^{{\rm{j2\pi }}f\tau }}{\rm{d}}f} \end{array} \right. \]

计算平均功率，令\(\tau = 0\)：

\[R_x(0) = \int_{-\infty }^\infty S_x(f) \rm{d}f \]

对于功率谱密度，我们最熟悉的莫过于噪声的PSD。例如热噪声，闪烁(1/f)噪声，散粒（shot）噪声。不同的噪声类型的PSD不一样。在电路中我们经常考虑在1单位电阻上的功率谱密度，所以在表示噪声的功率谱密度是经常使用V^2/Hz的单位。其实比较符合PSD定义的单位还是W/Hz。也有转化为对数坐标的dBm/Hz等。

当然扩展一下，在电学以外，比如振动力学中研究随机振动，会有位移PSD，速度PSD，加速度PSD，单位分别表示为m^2/Hz,(m/s)2/Hz，(m/s²⁾2/Hz或者g^2/Hz。看来能够研究的随机变量或信号都可以有PSD的概念。想想也是，其实就是把一些随机过程转化到频域通过PSD研究其确定性性质。

实信号的PSD分布通常关于频率0左右偶对称。这样的正负频率范围的功率谱密度称之为“双边”（Two-Side）PSD。当然我们也可以仅用正频率范围表示，这时候称之为“单边”（One-Side）PSD。如下图所示。

如果“单边”PSD关于某频率$f_0$左右对称，如果把该频率平移至频率0Hz处，这种功率谱我们称之为“单边带”（Single sideband）PSD。对于SSB PSD，我们可以进一步简化为仅用正频率表示。这时候称之为“双边带”（Double Sideband）PSD。

注意功率谱幅度的变化，将“双边”转化为“单边”时，功率谱密度需要加倍，也就是增加3dB。而从“单边”转化为“单边带”时，仅仅是载波频率的平移，幅度功率谱幅度并不改变。而从“单边带”转化为“双边带”时，功率谱密度也需要增加3dB。

参考链接1：功率谱、功率谱密度介绍及matlab代码 - splk的文章 - 知乎
参考链接2：matlab计算功率谱的四种方法 - CSDN
参考链接3：功率谱: 功率谱密度(PSD)、功率谱 - CSDN

10 雷达信号处理的一个小技巧

在做雷达信号仿真的时候，经常会碰到将发射信号延时一段时间的仿真，比如接收信号就是发射信号的时延。以前的处理办法是计算时延，然后想办法把信号平移过去。

最近在仿真机载雷达杂波建模的时候，看别人的代码，发现一种新的办法，这种办法本质是利用：任意信号\(f(t)\)与冲激函数卷积，等效于将\(f(t)\)平移至冲激函数所在时刻。代码示例如下：

%% 计算杂波点对应回波
Echo = zeros(1, Nprf*Np);
for i = 1:Np
    H = zeros(1, Nprf*Np);
    for j = 1:N_Clutter_block
        % 计算第i个慢时间采样点时候第j个杂波单元离雷达的距离
        Rij = norm(clutter_pos(j, :) - air_coord_instant(i, :));
        % 计算第i个慢时间采样点时候第j个杂波单元的回波的延迟采样点索引
        td_num = round(2*Rij/c*Fs)+(i-1)*Nprf;
        % 计算第j个杂波块的回波功率
        Pr_clutter_ij = Pt*AntGain(j)^2*lambda^2*RCS(j)/((4*pi)^3*Rij^4);
        % 计算回波幅度
        A_clutter_ij = sqrt(Pr_clutter_ij);
        % 将当前杂波点的回波贡献加到时间索引 td_num 处
        H(td_num) = H(td_num) + A_clutter_ij*exp(-1j*2*pi*f0*2*Rij/c);
    end
    % 前面只是得到了单个点的回波，这里将单点目标响应转换为完整的回波波形 # TODO 这种方法值得记录
    Echo_temp = ifft(fft(H).*fft(Source));
    Echo = Echo + Echo_temp;
end

for i = 1:Np
    % 对每个脉冲进行脉冲压缩，y(i, :) 为匹配滤波输出，即快时间脉冲压缩，去除一部分点？？？
    y(i, :) = ifft((fft(Echo(1+(i-1)*Nprf:i*Nprf)).*conj(fft(Source(1:Nprf)))));
end