标量/向量/矩阵的相关运算
矩阵乘法
知乎-矩阵求导术：上篇和下篇。本文基本上是这两篇文章内容的重新整理。
刘建平Pinard系列博客，这个博客主要用于查缺补漏
教材：《矩阵分析与应用》，作者张贤达
博客：矩阵求导总结
查询手册：The Matrix Cookbook
超级强大的在线矩阵求导工具：Matrix Calculus
2 信号处理中数学公式Matlab编程实现
2.1 Wigner分布的实现
2.1.1 理论介绍
首先，下面这部分内容来自书籍《现代信号处理教程（第二版）》，该书籍的相关信息请参考下面这两个链接：
| 现代信号处理教程（第二版）- 百度百科 |
| 现代信号处理教程（第二版）- 豆瓣 |
信号\(x(t)\)的自Wigner分布为：
信号处理算法要服务于实际工程，即最终的目的是要将它们应用于科研或工程的实际。因此，上面这种连续域无限长的理论公式必须要进行离散化和有限化才能在我们这个非连续有限的物理世界中实现，这时所遇到的问题就是信号的离散化及数据的有限长问题。
- 离散化
要离散化，就要对信号进行采样(把时间离散化)，对于信号\(x(t)\)，我们令采样间隔为\(T_s\)，即\(t = nT_s\)，则有：
进一步，要用计算机计算，连续的积分变量也要离散化，因此我们令\(\dfrac{\tau}{2} = kT_s\)，即\(\tau = 2kT_s\)，故而上式对\(\tau\)的积分变成对k的求和：
|
注意：上面这种描述我不知道是否准确，但是在写的过程中，还有如下一种思考： 信号\(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\)，则上式变为：
注意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\)的求和范围。
现在令：
并假定\(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问题，这个我可能穿插在代码部分，这里就先不展开了，这样表达式就变成了：
最后，看一下\(\omega\)的离散化，这里我呢也不深究，只说一下目前的简要理解，后期如果有其他想法再补充，首先我们根据复指数函数的性质可知：
也即，对于\(\omega\)变量来说，是按照\(\pi\)周期性变化的，因此，我们只需考察\(0~\pi\)内的数值即可，其他部分都是周期重复的，所以现在可把\(\omega=0~\pi\)区间内进行采样离散化，至于采样间隔多大，我也不知道，但是一般和离散信号离散点数一样即可吧，即把\(\omega=0~\pi\)区间划分\(N\)份：
所以，最终公式的完全离散化表达为：
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
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函数：
我们把DFT函数与离散化的WVD表达式对比，就很明显可以得到，把\(n\)视为固定的数，\(r_x(n, k)\)就可视为DFT函数中的\(x(k)\)，所以WVD这不就是变成了一个DFT了吗，最终的输出结果如下：
参考链接1：MATLAB辅助现代工程数字信号处理 第二版 作者 李益华 第6-10章_ 第8章.ppt - 原创力文档
2.2 离散傅里叶DFT的实现
2.2.1 连续傅里叶变换的离散化
有了上面的基础，下面我自己想试着写一下DFT的实现。
连续傅里叶变换的公式为：
首先，对积分变量t进行离散化：\(t = nT_s\)，则\(dt = (n+1)T_s-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(t)\)以\(T_s\)采样间隔得到离散化的采样序列，但是实际中不可能无限采样，故将采样序列的索引限制为：\(n = 0, 1, \cdots, N-1\)，或者说不是限制采样序列索引长度，而是我们实际就采样了\(N\)个点：
这就是离散时间傅里叶变换DTFT了。理论上，到此，我们就可以计算任意一个\(\omega\)处的\(F(\omega)\)值了，例如计算\(\omega = 0.401\)处的值，则只需把\(0.401\)代入公式然后交给计算机进行计算就行了，但是为了能有规律的计算不同\(\omega\)上的点并绘制出一条线便于观察，我们可以对\(\omega\)也进行离散化，也即：DTFT->DFT。
观察到，\(F(\omega)\)是周期性的，即：\(F(\omega)\)=\(F(\omega+2\pi)\)，因此，只需要在\([0, 2\pi]\)直接对其离散化就行了，离散多少个点呢，实际上可以随意，即：
上面的\(M\)就表示\(\omega\)在一个周期内离散化的点数，一般情况下(比如Matlab官方的fft函数)，都是默认\(M=N\)的，最终得到离散化后的公式：
至此，算是彻底完成了连续傅里叶变换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的点数直接相关的。
因此，波形频率分辨率可定义为：
其中，\(T\)表示实际有效的采样信号时长，不包括人为补零的时长。
同样，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. 直观解释：时间越长，分辨越清晰 想象你在听两个音调，比如钢琴上的两个键，它们的频率很接近（比如 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）是在时域信号末尾添加零，使信号长度从\(𝑁\)增加到\(𝑀\)。我们需要证明： 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为：
设信号\(x[n]\)与分别与冲激响应为\(h[n]=\mathrm{e}^{-j\frac{2\pi}{N}kn}(\text{ }k=0, 1, \cdots, N-1)\)的滤波器在时域卷积，则有：
此卷积在\(m=N\)时刻的输出为：
可以观察到，两者是完全一样的，这样就说明了信号\(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
