0
点赞
收藏
分享

微信扫一扫

雷达信号处理中的离散傅里叶变换(DFT)以及加窗


本文编辑:@调皮连续波,保持关注调皮哥,获得更多雷达学习资料和建议!


大家好,我是调皮哥,今天继续给大家分享干货,助力大家轻松、快乐、有方向地学习雷达。

离散傅里叶变换 (DFT) 计算离散时间信号的频谱,一般会采用快速傅里叶变换(FFT)的高效算法来计算=但单纯地使用FFT通常会产生频谱泄露,从而导致频谱的计算不准确。

如果将窗函数与DFT结合使用,则可以大大减少频谱泄漏。在本文中,我们将了解产生频谱泄漏的原因,然后介绍窗口函数并展示采用窗函数法如何改善频谱,让各位读者更清晰的理解。

本文的内容可以用于任何采用数字信号处理的地方,例如雷达、声呐、通信等领域。

1.频谱泄露产生原因

1.1 频谱分辨率不足

对于离散时间序列 ,DFT 定义为:

()其中,   是离散时间序列   是  和  的样本数,也叫做采样点数,  是时间索引 : , 

我们看到,根据定义,DFT 适用于N个样本的有限长度信号,对于 的采样时间,离散时间变量由下式给出: 

 是复数,因为傅里叶变换会引入欧拉公式中的虚数  ,傅氏变换后得到的复数,实部就代表该频率下的余弦信号分量,虚部就代表该频率下的正弦信号分量。当  为实函数的时候,其频谱函数的实部为偶函数,虚部为奇函数。当 



如果 是实信号,根据三角公式: 

那么,  是正负频谱对称的双边谱,w是角频率。



DFT 能够计算持续时间为N个样本或更短的任何信号的频谱,由于我们感兴趣的许多信号没有明确定义的边界,因此我们最终采用长度为N的一段信号来计算 DFT。

例如,考虑8和9Hz的正弦波,采样频率为 , 个样本点,为了方便令L=N。下面给出了8Hz正弦波信号模型以及计算DFT 的 Matlab 代码,其中我们使用 Matlab 中的fft()函数来计算公式(1) 的 DFT。该函数计算离散频谱索引值,即 :

fs= 128 ;                      % Hz 采样频率
Ts= 1 /fs;                    % s 采样时间
L= 64 ;                        % 样本总数
f0= 8 ;                        % Hz 正弦频率
n= 0 :L -1 ;                        % 时间索引
x= sin ( 2 * pi *f0*n*Ts)* 2 /L;    % 离散时间正弦
X = fft(x,L);                       % 离散频谱通过 DFT 
Xmag=abs(X);                
k= 0 :L -1 ;                       % 频率索引
f= k*fs/L;                        % Hz 离散频率


在上述程序中,我们将正弦信号的幅度缩放了2/L,这样可以频谱分量的幅度为1。请注意,对于复向量X,Matlab 采用abs()函数计算X的每个元素的幅值大小。


离散时间正弦曲线  绘制在图 1 的顶部,我们看到在 64 个样本中正好是四个周期的正弦波。对于从 0 到  的离散频率,DFT 绘制在中间。  是纯虚数,这是因为  是奇函数。图1显示了  的大小,频谱绘制在 0 到 

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_卷积_02


图1. 正弦曲线的 DFT。f 0 = 8 Hz,f s = 128 Hz,L = 64 个样本。顶部:离散时间正弦波。中间:X(f) 从0到fs的频谱,底部:|X(f)| 对于 f = 0 到 f s /2。


如公式 3 所示,离散频率或区间的间隔为 

f0 = 9 ; % Hz

正弦曲线绘制在图 2 的顶部,现在 64 个样本中有 4.5 个周期。  的大小绘制在底部图中。明显看到,频谱分布在许多的频率单元上:这就是前面提到的频谱泄漏现象。另外,还存在一个明显的问题:在 9 Hz 处没有频率单元。

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_时域_03


图 2. 正弦曲线的 DFT。f 0 = 9 Hz,f s = 128 Hz,L = 64 个样本。顶部:离散时间正弦波。底部:|X(f)| 对于 f = 0 到 f s /2。


上述结论表示,当截取的信号长度很短时,频谱的分辨率存在不足,无法满足区分相邻单元频率的目的,当信号长度L提高时,情况会有很大的改善。当:

L= 64*2 ;                        % 样本总数

得到下图,既可以区分8Hz和9Hz。

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_卷积_04

1.2.分段处理

频谱泄露原因:实际应用中,我们对模拟信号进行采样然后进行分段处理,这里的“分段”处理,相当于对信号x(n)乘了一个R(N)的N点矩形窗,R(N)的频谱如下,造成了信号频谱能量向两边扩散。

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_时域_05

图3 显示了捕获L个正弦波样本的视图,顶部的图是正弦波的样本,中间的图是一个矩形函数,其中 L = 64 个样本等于 1.0,底部的图显示了我们逐个样本将正弦波乘以矩形函数时的结果。

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_正弦波_06


图 3. 矩形窗口函数的窗口。顶部:要捕获的正弦波。中部:矩形窗函数。底部:加窗正弦波。


我们建立了一个可以捕获 L 个正弦波样本点的矩形函数窗口,每当捕获信号的持续时间超过  时,就会发生窗口化。这对于现实世界的信号来说很常见:在数学上,加窗是正弦波与窗函数的逐个元素相乘,即窗函数对正弦波的调制。

使用 DFT 时,我们可以通过在正弦波信号末端补零(图 3 中 n = 63 以上的零值样本)更清楚地看到矩形窗口的效果。这将为我们提供一个较小的  离散频率步长,其中N是样本总数,包括补零。

下面的 Matlab 代码像以前一样生成 64 个正弦波样本,  。然后补零以获得总共 8*64 = 512 个样本。然后我们找到这个补零信号的 DFT。(注意:虽然我们在这里明确地补零,但如果 x 的长度小于 N,Matlab FFT 函数将自动对信号 x补零)。

fs= 128 ;                     % Hz 采样频率
Ts= 1 /fs;                    % s 采样时间
L= 64 ;                       % 正弦波样本数
f0= 8 ;                       % Hz 正弦频率
n= 0 :L -1 ;                    % 时间索引
u= sin ( 2 * pi *f0*n*Ts)* 2 /L;    % 离散时间正弦曲线
R= 8 ;                              % 零填充因子
N=R*L;                             %补零
x= [u zeros( 1 ,NL)] 的总样本百分比;         % 零填充信号
X = fft(x,N);                   % 离散光谱通过 DFT 
Xmag= abs (X);                  % X 
k= 0 :N -1 ;                    % 频率索引
f= k*fs/N;                      % Hz 离散频率

补零信x绘制在图 4 的顶部,频谱X的幅度绘制在底部,我们将频率轴限制为0到32Hz 

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_卷积_07


图 4. 正弦波的矩形补零窗口,f0 = 8Hz,fs = 128 Hz。顶部:具有零填充的窗口正弦波。底部:加窗正弦波


现在让:

f0 = 9 ; % 赫兹

补零信号x绘制在图 5 的顶部,频谱X的幅度绘制在底部。我们现在可以看到,图 2 的频谱泄漏分量都落在频谱的旁瓣峰值处。此外, 的旁瓣幅度随频率下降非常缓慢,矩形窗口已经把频谱弄得一团糟,存在很多其他分量的频率!

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_正弦波_08


图 5. 正弦波的矩形补零窗口,f0 = 9 Hz,fs = 128 Hz。顶部:具有补零的窗口正弦波。底部:加窗正弦波的幅度谱。


2.矩形窗的频谱

图 4 和图 5 显示了加窗正弦曲线的频谱,现在让我们找到矩形窗口本身的频谱,当然,这其实在上述文中简单说过,下面是详细的过程。

下面的 Matlab 代码生成一个 L= 64 个的窗口,然后补零以获得总共 N = 8*64 = 512 个样本,如图 6 的上图所示,找到补零后的矩形窗口,然后把  的中心频率移到0的位置上,我们在图 6 的底部图中绘制 了频率单位为  的偏移幅度响应。

fs= 1 ;                    % Hz 采样频率
L= 64 ;                    % 矩形窗长度
R= 8 ;                     % 零填充因子
N= R*L;                   具有零填充
x= [ones( 1 ,L) zeros( 1 ,NL)]/L 的总样本百分比;% 具有零填充的矩形窗口
X = fft(x,N);             % 离散光谱通过 DFT 
Xmag= abs (X);             X 
k= 0 :N -1 ;                 % 频率指数
f= k*fs/N;                % Hz 频率
% 频谱以 0 Hz 为中心
Xmag_shift = fftshift(Xmag);    % 交换 Xmag 的右半部分和左半部分
f_shift= f- fs/ 2 ;               % 将频率范围转移到 -fs/2:fs/2

图 6 底部绘制的频谱幅度的函数形式由下式给出: 

()最后,请注意主瓣带宽与 L 成反比:零点到零点的带宽为  ,因此当截取的信号长度越长,带宽越窄,频谱分辨率越高。

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_卷积_09


图 6. 矩形窗口的频谱。顶部:补零窗口。底部:频谱幅度以 0 Hz 为中心。


将图4或图5与图6进行比较,为什么图 6 的窗口频谱是对称的,而窗口正弦波的频谱却不是?时域中的乘法(调制)等价于频域中的卷积,窗口与单个频谱分量的卷积应该是对称的,连续时间信号符合这种情况。但对于离散时间信号而言“时域中的乘法相当于频域中的循环卷积,因此落在0到fs范围之外的卷积分量,在线性卷积环绕下又叠加到循环卷积下的频谱上,从而导致图 4 和 图5 存在不对称性。


注:由于离散傅里叶变换(DFT)的实质是周期序列变换到频域的描述,可以证明:两个有限长序列在时域的循环卷积,其DFT等于在频域两个序列相应的DFT的乘积。


图 B.1 的上图显示了具有补零的8 Hz正弦波的频谱 U,采样率为 128 Hz,频谱绘制在 0 到 fs 之间。中间图显示了具有补零的矩形窗口的幅度谱  ,下图显示了 U 和 X 的循环卷积的幅度 

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_卷积_10

3.使用窗函数改进频谱

频谱泄露问题没有完美的解决方案,但我们可以可以采用算法改进。鉴于矩形窗口对 DFT 有如此有害的影响,我们不禁要问:有什么办法可以改进吗?

答案是:有,但需要改变窗口的形状,因此我们需要一个比矩形窗口具有更窄的窗口。

图7 的顶部绘制了正弦波的样本。在之前的图 3 中,中间的图显示了一个具有 L = 64 个样本的窗函数,该函数具有平滑的形状,其幅度开始和结束接近于零,图7显示了我们将逐个正弦波样本乘以窗函数的结果。

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_卷积_11

图 7 所示的窗口函数是常用的 Hanning 窗口:

()例如,对于 L = 8,我们有:w = [0, 0.1464 ,0.5, 0.8536 ,1,0.8536 ,0.5 ,0.1464]  请注意,公式 5 由 Matlab 函数hann(L,'periodic') 实现。

下面让我们使用汉宁窗获取正弦波的频谱。同样,我们将使用长度为 L = 64 个样本的 9Hz 正弦波,其中 fs= 128Hz。在下面的 Matlab 代码中,我们创建了一个长度为L的 Hanning 窗口,然后逐个样本将正弦与窗口相乘(Matlab 运算符 .* 执行此操作)。然后我们对加窗正弦波 u 进行补零以获得总长度为 512 个样本的信号 x(x 计算中的因子 2 将频谱幅度缩放为最大值 1.0)。

图 8 顶部显示了加窗的正弦波 x。请记住,补零可用于提高频谱的频率分辨率,但不是计算 DFT 的要求。最后,我们取 x 的 DFT。得到的频谱幅度显示在图 8 的底部,我们将频率轴限制为 0 到 32 Hz (fs/4)。

fs= 128 ;                      % Hz 采样频率
Ts= 1 /fs;% s 采样时间
f0 = 9 ;                       % Hz 正弦频率
L= 64 ;                        % 窗口长度
R= 8 ;                         % 零填充因子
N= R*L;                       % 具有零填充的总样本
n= 0 :L -1 ;                     零填充之前的 % 时间索引
win= .5 *( 1 -cos ( 2 * pi *n/L));    % 汉宁窗长度 L 
u= sin ( 2 * pi *f0*n*Ts)* 2 /L;     % 离散时间正弦
u_win=win.*u;                % 加窗正弦
x= 2 *[u_win zeros( 1 ,NL)];    % 加窗的零填充正弦
X = fft(x,N);                 % 离散光谱通过 DFT 
Xmag= abs (X);                 X 
k= 0 :N -1 ;                     % 频率指数
f= k*fs/N;                    % 赫兹频率

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_正弦波_12


图 8. 具有补零的汉宁加窗正弦波。f0 = 9 Hz,fs = 128 Hz。顶部:具有补零的窗口正弦波。底部:加窗正弦波的幅度谱。


将图8的频谱与图 5 的频谱进行比较,我们看到旁瓣要低得多。然而,由于窗口函数的调制,我们仍然有接近 9 Hz 的频谱展开(它甚至比矩形窗口宽一些),因为任何窗口函数都存在这样的情况。

我们可以使用以下方法在 dB 标度上绘制两个频谱:

XdB= 20 * log10 (Xmag);

dB 频谱如图 9 所示。同样,如果您只对捕获的正弦波进行 DFT,则可以得到标记为“rectangular”的频谱:它是由矩形窗函数调制的正弦波频谱,标有“hanning”的频谱是在进行 DFT 之前用汉宁窗调制正弦波的结果。

如前所述,使用汉宁窗极大地改善了旁瓣,但主瓣却展宽了,-3 dB 带宽几乎为 3 Hz。请注意,主瓣带宽与 L 成反比。因此,以增加 L 为代价,我们可以减少该带宽,比如在图 10 绘制了 f0= 9.25 Hz 和 L = 256 的窗口正弦波频谱,汉宁窗看起来更接近一个理想的正弦频谱。

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_卷积_13


图 9. 使用 Hanning 和矩形窗的正弦波频谱。f0 = 9 Hz,fs = 128 Hz,L= 64,N = 8*L = 512。


我们得到一个经验可以表述如下:假设获取的信号是未知的,在使用汉宁窗进行了补零DFT后,得到了如图 10 所示的频谱,该图并不是信号本身的一个特征,其旁瓣是由窗函数引起的。图 10 的频谱使用了 DFT 长度 N = 8*L 的补零。

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_时域_14


图 10. 使用 Hanning 和矩形窗口的正弦波频谱。f0 = 9.25 Hz,fs = 128 Hz,L= 256,N = 8*L = 2048。


补零对于获得有效频谱并不是必需的。例如,图 11 显示了没有补零的频谱,其中L = N = 256,离散频率分辨单元fs /L = 128/256 = 0.5 Hz(因此,我们选择的 f0 = 9.25 Hz落在离散频率分辨单元的中间)。两个频谱的峰值都低于补零情况的峰值,发生这种情况是因为 f0落在补零情况下的频率单元上。

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_卷积_15


图 11. 使用 Hanning 和矩形窗的正弦波频谱。f0 = 9.25 Hz,fs = 128 Hz,L = N = 256(无补零)。


比较窗函数本身的频谱,图 12 绘制了矩形函数和汉宁函数的 dB 频谱,以 0 Hz 为中心。两个函数的带宽都与 L 成反比,窗函数的两个最重要的属性是带宽(即噪声带宽)和旁瓣电平

矩形窗口具有非常高的旁瓣,相对于频率下降得相当缓慢(第一个旁瓣在主瓣 -13 dB 处),汉宁窗的旁瓣较低(第一旁瓣为 -31.5 dB),但带宽较宽。通常,具有较低旁瓣电平的窗函数具有较宽的带宽,当试图在不同频率的大信号存在的情况下显示小信号时,旁瓣电平特别重要

雷达信号处理中的离散傅里叶变换(DFT)以及加窗_时域_16


图 12. Hanning 和矩形窗的光谱


有许多窗口函数的边带比 Hanning 窗口 低,Kaiser 窗具有可调节的旁瓣电平。

总而言之,以时域的乘积等于频域的卷积这个基本性质作为基础,一个矩形窗函数对其进行DFT变换过后,其频域响应  可以由  函数来近似,所以  就是那个所谓旁瓣效应的函数。


那么为什么会有泄漏?通俗地讲,其实就是当我们在一段无限长的序列中取一段有限长度的序列,其实我们人为已经对这个无限长序列加了一个矩形窗,也就是这个无限长序列的频谱已经和一个矩形窗函数的频谱做了卷积了,自然输出的这个有限长序列输出的频谱就变形了。而如果所取的这个有限长序列的N个点,恰好序列内有k个周期,那么由上面的公式可知,就没有泄漏发生。日常应用中,我们对一段信号进行分帧处理的时候,频谱泄漏会影响分析,所以才会用到窗函数,而经过窗函数处理的时域信号,其初始点和结束点的时域振幅都接近0。


4.雷达信号处理中的加窗处理

如上所述,为了改进频谱估计结果,所以需要加窗,因为 DFT 中使用的周期性假设要求估计的频率是频谱分辨率的整数倍,这也就是为什么上述8Hz能够估计出来,而9Hz会出现泄露。但是一般频率是不可能都是保证频谱分辨率的整数倍,结果是频谱不连续,进而导致能量扩散到多个频率区间,表现为旁瓣(泄露),这个问题可以通过使用窗函数来改善。

窗函数从零附近或零开始,然后在采样数据序列的中心增加到最大值,然后再次减小 。从卷积原理来看,加窗具有平滑信号频率响应的作用,与没有加窗的频谱相比旁瓣的幅度会降低。然而,加窗也会导致主瓣变宽,因此我们需要权衡利弊,并不是所有技术都只有好处没有坏处。

时域加窗: 

 频域卷积: 

 4.1 雷达信号处理加窗

在脉冲压缩、相干积累、数字波束形成等大多数雷达信号处理中,输出均为 sinc函数的形式,其峰值主副瓣比只有 -13.4dB。在工程实际中,高副瓣将影响对目标的检测,容易产生虚警或受到干扰等。因此,在雷达信号处理中,经常采用加窗技术抑制副瓣,提高峰值主副瓣比,进而提高信噪比SNR。

在雷达信号处理中,经常会执行DFT运算,比如FMCW雷达距离估计、速度估计,以及3D-FFT角度估计、MTD等算法中,而执行DFT就会发生不可预料的频谱泄露。



(1)脉冲压缩、相干积累、数字波束形成

(2)FMCW雷达距离估计(1D-FFT)

(3)FMCW雷达速度估计(2D-FFT)

(4)3D-FFT角度估计(3D-FFT)



4.2 典型的窗函数

常用的窗有:矩形窗、汉明窗(Hamming)、汉宁窗(Hanning)、布莱克曼(Blackman)等,我们可以通过Matlab很容易的对序列进行加窗,也可以简单的看到窗函数的频率响应。

(1)矩形窗

优点是主瓣比较集中,缺点是旁瓣较高,并有副旁瓣,导致加窗过程中带进了高频干扰和频谱泄漏。

(2)汉宁窗

又称升余弦窗,汉宁窗使主瓣加宽并降低,旁瓣则显著减小,从减小泄漏观点出发,汉宁窗优于矩形窗。但汉宁窗使主瓣加宽,相当于频谱带宽加宽,导致频率分辨力下降。汉宁窗函数的最大旁瓣值衰减-31dB,但是主瓣宽度比矩形窗函数的主瓣宽度增加了1倍。

(3)汉明窗

余弦窗的一种,又称改进的升余弦窗,与汉宁窗的加权系数不同。汉明窗加权的系数能使旁瓣达到更小。汉明窗的第一旁瓣衰减为-42dB,其旁瓣衰减速度比汉宁窗衰减速度慢。汉明窗与汉宁窗都是很有用的窗函数。

除了上述这些之外,还有一些别的窗函数,下面的文章中有谈到:https://zhuanlan.zhihu.com/p/551373297

关于窗函数还有很多内容需要讲解,今天暂时就到这里,后面有机会再继续。

参考资料



【1】https://www.dsprelated.com/showarticle/1433.php




举报

相关推荐

0 条评论