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

矩阵乘法

xTΛx=(x1,,xn)(λ1λn)(x1xn)=i=1nλixi2

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

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

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

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

2.1 Wigner分布的实现

2.1.1 理论介绍

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


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

Wx(t,Ω)=x(t+τ2)x(tτ2)ejΩτ dτ

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

  • 离散化

要离散化,就要对信号进行采样(把时间离散化),对于信号x(t),我们令采样间隔为Ts,即t=nTs,则有:

Wx(t,Ω)=x(nTs+τ2)x(nTsτ2)ejΩτ dτ

进一步,要用计算机计算,连续的积分变量也要离散化,因此我们令τ2=kTs,即τ=2kTs,故而上式对τ的积分变成对k的求和:

Wx(t,Ω)=2Tsk=x(nTs+kTs)x(nTskTs)ej2kΩTs

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

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

Wx(1,0.5)=x(1+τ2)x(1τ2)ej0.5τ dτ

这时候我们只要把τ离散化,比如,令τ2=kTs,即τ=2kTs,则有:

Wx(1,0.5)=2Tsk=x(1+kTs)x(1kTs)ej2k×0.5Ts

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

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

Wx(t,Ω)=2Tsk=x(nTs+kTs)x(nTskTs)ej2kΩTs

PS:还有一点要补充地是,为什么t离散化为t=nTsτ离散化为τ2=kTs呢,为什么不是τ2=kTτ呢?

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

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

Wx(n,ω)=2k=x(n+k)x(nk)ej2kω

注意1:这里n=012是离散化了,但是ω=ΩTs还是没有离散化,因为虽然Ts离散化了,但是Ω是连续的,一个连续的值Ω乘以一个离散固定值Ts仍然是连续值。

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

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

  • 有限性

现在还剩下两个问题需要解决,① 一是频率ω的离散化;② 二是确定式中k的求和范围。

现在令:

rx(n,k)=x(n+k)x(nk)

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

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

这里直接通过具体计算来展示,比如我们令N=6,则有:

下标索引n rx(n,k)的值 k的取值范围
n=0 rx(0,k)=[x0x0] k=0
n=1 rx(1,k)=[x0x2,x1x1,x2x0] k=1,0,1
n=2 rx(2,k)=[x0x4,x1x3,x2x2,x3x1,x4x0] k=2,1,0,1,2
n=3 rx(3,k)=[x1x5,x2x4,x3x3,x4x2,x5x1] k=2,1,0,1,2
n=4 rx(4,k)=[x3x5,x4x4,x5x3] k=1,0,1
n=5 rx(5,k)=[x5x5] k=0

根据上表可以看出,由于下标采用的是n=0,1,,N1,因此限制了不同n取值下k的取值范围,总之:一条准则,k的取值范围不能使信号x?的下标出现超出n=0,1,,N1的范围。

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

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

下标索引n rx(n,k)的值后填充0 k的取值范围
n=0 rx(0,k)=[x0x0,0,0,0,0,0] k=0
n=1 rx(1,k)=[x0x2,x1x1,x2x0,0,0,0] k=1,0,1
n=2 rx(2,k)=[x0x4,x1x3,x2x2,x3x1,x4x0,0] k=2,1,0,1,2
n=3 rx(3,k)=[x1x5,x2x4,x3x3,x4x2,x5x1,0] k=2,1,0,1,2
n=4 rx(4,k)=[x3x5,x4x4,x5x3,0,0,0] k=1,0,1
n=5 rx(5,k)=[x5x5,0,0,0,0,0] k=0

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

下标索引n rx(n,k)的值 k的取值范围
n=0 rx(0,k)=[0,0,x0x0,0,0,0] k=0
n=1 rx(1,k)=[0,x0x2,x1x1,x2x0,0,0] k=1,0,1
n=2 rx(2,k)=[x0x4,x1x3,x2x2,x3x1,x4x0,0] k=2,1,0,1,2
n=3 rx(3,k)=[x1x5,x2x4,x3x3,x4x2,x5x1,0] k=2,1,0,1,2
n=4 rx(4,k)=[0,x3x5,x4x4,x5x3,0,0] k=1,0,1
n=5 rx(5,k)=[0,0,x5x5,0,0,0] k=0

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

Wx(n,ω)=2k=(N+1)/2(N+1)/2rx(n,k)ej2kω

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

Wx(n,ω+π)=2krx(n,k)ej2k(ω+π)=Wx(n,ω)

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

ω=0,1Nπ,2Nπ,,N1Nπ

所以,最终公式的完全离散化表达为:

Wx(n,l)=2k=(N+1)/2(N+1)/2rx(n,k)ej2kπNl,l=0,1,,N1

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)');
  • 然后,调用一个函数来计算rx(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)对应的Rx矩阵中间添加了一个全0列矩阵,不这样添加的话,对于偶数信号的Rx矩阵就不是方阵了。
② 第2:上面的图是保留函数最后3行的结果,自己也可以尝试去掉最后3行运行查看结果,差异就是最后一列在前在后的区别。
③ 第3:通过与我自己用visio绘制的图形可知,该函数的输出结果是采用方式2,并且经过fftshift之后的结果,fftshift就相当于把k=[(N+1)/2,(N+1)/2]变成k=[1,2,...,N2,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函数:

X(l)=DFT[(x(k)]=k=0N1x(k)ej2kπNl,l=0,1,2,...,N1

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

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

2.2 离散傅里叶DFT的实现

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

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

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

F(Ω)=f(t)ejΩt dt

首先,对积分变量t进行离散化:t=nTs,则dt=(n+1)TsnTs=Ts,故有:

F(Ω)=n=f(nTs)ejΩnTsTs

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

F(ω)=Tsn=f(n)ejnω

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

F(ω)=Tsn=0N1f(n)ejnω

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

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

0,2πM,22πM,,(M1)2πM

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

F(k)=Tsn=0N1f(n)ejnk2π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.00MHz1.05MHz两个频率成分合成的正弦波实信号。

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

  • 两种分辨率的辨析

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

因此,波形频率分辨率可定义为:

ΔRw=1T

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

同样,FFT分辨率可以定义为:

ΔRf=fsNfft

其中,fs为采样频率 (the sampling frequency),Nfft为FFT的点数。

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

例1:信号时长为1秒,采样率1000Hz,FFT点数1000

波形频率分辨率:ΔRw=1T=1/1=1Hz(理论上可分辨间隔1Hz的两个频率)

FFT分辨率:ΔRf=fsNfft=1000/1000=1Hz(频谱中每1Hz显示一个频点)

例2:信号时长为0.5秒,采样率1000 Hz,FFT点数2000(补零)

波形频率分辨率:ΔRw=1T=1/0.5=2Hz(实际只能分辨间隔2Hz的频率)

FFT分辨率:ΔRf=fsNfft=1000/2000=0.5Hz(频谱每0.5Hz显示一个频点,但实际无法分辨0.5Hz间隔的频率!)

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

现象:通过补零增加FFT点数(如例2),频谱显示的频点间隔更小(FFT分辨率提高)。

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

验证:若两个频率间隔小于1T,即使补零到Nfft,频谱主瓣仍会重叠,无法分辨两个信号。

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

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

提升FFT分辨率:仅需增加Nfft或降低fs ,但这只是优化显示效果,并没有很特别的意义。

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

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

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

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

对于一个连续时间信号x(t),其傅里叶变换为:X(f)=x(t)ej2πft dt,在实际应用中,我们通常处理的是有限长度的信号(即截断信号),假设信号持续时间为T。此时,信号可以表示为:

xT(t)=x(t)w(t)w(t)={10tT0otherwise

对这个有限长度信号进行傅里叶变换,得到的频谱为:XT(f)=0Tx(t)ej2πft dt,也就相当于加了窗函数w(t)

窗函数的频域影响:矩形窗w(t)的傅里叶变换是sinc函数:W(f)=Tsinc(fT)=Tsin(πfT)πfsinc函数的主瓣宽度决定了能够分辨的最小频率间隔。

sinc函数的主瓣宽度(从零点到零点的距离)是2T,而主瓣的半宽(中心到第一个零点的距离)是1T。根据卷积定理,时域乘积对应频域卷积:

XT(f)=X(f)W(f)

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

假设信号x(t)是两个频率接近的正弦波之和:x(t)=cos(2πf1t)+cos(2πf2t),它的频谱 X(f) 是两个冲激函数:𝛿(𝑓𝑓1)+𝛿(𝑓𝑓2)。经过有限时长T的截断后,频谱变成:

XT(f)=X(f)W(f)=[𝛿(𝑓𝑓1)+𝛿(𝑓𝑓2)]Tsinc(fT)

结果是两个sinc函数,分别以f1f2为中心,宽度约为1T。① 如果|f1f2|<1T :两个sinc函数重叠严重,无法区分两个峰值。② 如果|f1f2|>1T:两个sinc函数的主瓣分开,可以分辨出两个频率。

因此,频率分辨率定义为1T,它是能够分辨两个频率的最小间隔。

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

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

1. 问题的核心

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

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

(1)原始信号及其频谱

假设有一个离散信号𝑥[𝑛],长度为𝑁,采样频率为𝑓s。其离散时间傅里叶变换(DTFT)为:

X(ω)=n=0N1x[n]ejωn

X(ω)是连续的频域函数,ω是归一化频率(范围[𝜋,𝜋])。离散傅里叶变换(DFT)是对X(ω)𝑁个等间隔点采样:

X[k]=X(ω)|ω=2πkn=n=0N1x[n]ej2πkn/N,k=0,1,2,,N1

频率间隔为Δ𝑓=𝑓𝑠/𝑁

(2)补零后的信号

补零后,信号变为𝑥pad[𝑛],长度为𝑀(M>N)

xpad[n]={x[n],0n<N0,Nn<M

其DFT为:

Xpad[k]=k=0M1xpad[n]ej2πkn/M

由于xpad[n]Nn<M时为零,则上面的补零DFT可化简为:

Xpad[k]=k=0N1x[n]ej2πkn/M

频率间隔变为fsM<fsN,频谱点更密集。我们需要证明Xpad[k]𝑋(𝜔)的插值结果,且插值函数是sinc

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

Xpad[k]=k=0N1x[n]ej2πkn/M

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

Xpad[k]=X(ω)|ω=2πk𝑀

这表明Xpad[k]X(ω)𝑀个等间隔点(𝜔=2𝜋𝑘𝑀)上的采样。我们需要找到Xpad[k]与原始𝑋[𝑘]的关系。

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

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

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

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

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

x[n]=1Nk=0N1X[k]ej2πkn/N

将其代入𝑋(𝜔)的定义:

X(ω)=n=0N1x[n]ejωn=n=0N1(1Nk=0N1X[k]ej2πkn/N)ejωn

交换求和顺序:

X(ω)=1Nk=0N1X[k]n=0N1ej(ω2πk/N)n

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

内层是一个几何级数:

n=0N1ej(ω2πkN)n=1ej(ω2πkN)N1ej(ω2πkN)=ej(ω2πkN)N12sin((ω2πkN)N2)sin(ω2πkN2)

这正是sinc函数的形式(归一化后)。因此:

X(ω)=k=0N1X[k]1Nej(ω2πkN)N12sin((ω2πkN)N2)sin(ω2πkN2)

这就是用𝑋[𝑘]重建𝑋(𝜔)的插值公式,插值核是周期化的sinc函数(也叫Dirichlet核)。

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

𝜔=2𝜋𝑘𝑀代入:

Xpad[k]=X(ejω)|ω=2πkM=l=0N1X[l]1Nsin((2πkM2πlN)N2)sin(2πkM2πlN2)ej(2πkM2πlN)N12

整理后:

Xpad[k]=l=0N1X[l]sinc(πN(k/Ml/N)π)ejϕ

其中,sinc(x)=sin(πx)πx,相位项𝑒𝑗𝜙=𝑒j(2πkM2πlN)N12=𝑒jπ(N1)(kMlN)是由于窗口对齐产生的。

证明步骤6:为什么是sinc函数?

采样定理:由香农采样定理,从离散样本重建连续信号的理想插值函数是sinc。若频谱𝑋(𝜔)是带限的(带宽小于𝑓𝑠2),sinc插值能完美重建

补零的本质:补零后,频谱𝑋pad[k]𝑋(𝜔)在更密集点的采样,而这种插值过程正好对应sinc函数的卷积效应

理想情况下,sinc插值公式为:

Xpad[k]=l=0N1X[l]sinc(πN(k/Ml/N)π)ejϕ

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

3. 物理意义:为什么sinc函数?

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

🚖 补零的频域效应:时域补零相当于频域与一个更宽的矩形窗卷积,其结果是用sinc函数插值原始频谱。

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

“相位项𝑒𝑗𝜙是由于窗口对齐产生的”的解释:

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

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

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

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

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

信号x[n]的离散傅里叶变换DFT为:

X[k]=n=0N1x[n]ej2πNkn=n=0N1x[n]ej2πNkn×1=n=0N1x[n]ej2πNkn×ej2πNkN=n=0N1x[n]ej2πNk(Nn)

设信号x[n]与分别与冲激响应为h[n]=ej2πNkn( k=0,1,,N1)的滤波器在时域卷积,则有:

y(m)=n=0N1x[n]h[mn]

此卷积在m=N时刻的输出为:

y(N)=n=0N1x[n]h[Nn]=n=0N1x[n]ej2πNk(Nn)

可以观察到,两者是完全一样的,这样就说明了信号x[n]的DFT在k0处的值X[k0]等于信号x[n]与冲激响应为h[n]=ej2πNk0n的滤波器在时域卷积。同时,我们又知道时域卷积等于频域相乘,而h[n]=ej2πNk0n的数字频谱是一个sinc函数,因此,对于k=0,1,2,,N1构成一个窄带滤波器组,中心数字角频率分别是:2πNk,每个窄带滤波器对等于其中心频率的单音信号的增益为: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

posted @   博客侦探  阅读(36)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 地球OL攻略 —— 某应届生求职总结
· 周边上新:园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源!
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】
点击右上角即可分享
微信分享提示