频谱泄露、栅栏效应、补零实验

理想情况


假设一个信号由两个频率不一样的余弦波组成 c o s ( 2 π f t ) = c o s ( 2 π f 1 t ) + c o s ( 2 π f 2 t ) cos(2\pi ft)=cos(2\pi f_1t)+cos(2\pi f_2t) cos(2πft)=cos(2πf1t)+cos(2πf2t)组成,其中 f 1 = 50 H z , f 2 = 65 H z f_1=50Hz,f_2=65Hz f1=50Hz,f2=65Hz我们以采样率 f s = 500 H z f_s=500Hz fs=500Hz对该信号进行采集,采样时间长度 T = 0.4 s T=0.4s T=0.4s,对该信号做FFT后得到的信号如下图所示,此时在频谱上可以清楚地找到两个信号频率对应的位置。
在这里插入图片描述
代码:

%% 理想情况
% 当采样时间正好为周期的整数倍时,没有频谱泄露;
% 因为序列时实值信号,所以FFT结果关于N/2对称,fftshift之后关于0对称
f0 = 50; f1 = 65; fs = 500; T = 0.4;
S = @(t) cos(2*pi*f0.*t) + cos(2*pi*f1.*t);
t = 0:1/fs:T-1/fs;
y = S(t);

nsamps = numel(y);
y_spec = fft(y);
y_spec = fftshift(y_spec)./nsamps;
f = ((0:nsamps-1)-nsamps/2)/T; 

figure(1);
subplot(2,1,1);plot(y,'linewidth',1.5); xlabel('time (s)'); ylabel('y');title('理想情况');
subplot(2,1,2);plot(f, abs(y_spec),'linewidth',1.5); xlabel('frequency (Hz)'); ylabel('Amplitude');
axis([0,inf,0,inf]);

频谱泄露


频谱泄露是对无限长信号进行截取操作时引入的旁瓣干扰,其在频域上的表现为除了主要频率之外还有不希望存在的频率成分。当对非周期信号进行截取操作时,无论怎么截取都会引起频谱泄露当对周期信号进行截取操作时,如果截取时间长度为原信号周期的整数倍,则不会产生频谱泄露反之则会。(因为DFT默认对截断的时域信号进行了周期延拓,在整数倍周期处截断再进行延拓可以完整地恢复出原始信号,反之如果在非周期处截断信号并延拓得到就不是原先的信号了)。

这里依然以上述信号为例。假设此时采样时间变成了0.3s,那么在该时间窗口内, c o s ( 2 π f 1 t ) cos(2\pi f_1t) cos(2πf1t)可以被截断为整数倍信号周期,然而 c o s ( 2 π f 2 t ) cos(2\pi f_2t) cos(2πf2t)则不是整数倍的信号周期。此时对采样信号做FFT结果如下。可以看到频谱上 f = 50 H z f=50Hz f=50Hz处幅度大于 f = 65 H z f=65Hz f=65Hz的幅度, c o s ( 2 π f 2 t ) cos(2\pi f_2t) cos(2πf2t)的能量被泄露到了其他频点上。
在这里插入图片描述
代码:

%% 频谱泄露
% 当采样时间不是周期的整数倍时,发生了频谱泄露;
f0 = 50; f1 = 67; fs = 500; T = 0.3;
S = @(t) cos(2*pi*f0.*t) + cos(2*pi*f1.*t);
t = 0:1/fs:T-1/fs;
y = S(t);

nsamps = numel(y);
y_spec = fft(y);
y_spec = fftshift(y_spec)./nsamps;
f = ((0:nsamps-1)-nsamps/2)/T; 

figure(2);
subplot(2,1,1);plot(y,'linewidth',1.5); xlabel('time (s)'); ylabel('y');title('频谱泄露');
subplot(2,1,2);plot(f, abs(y_spec),'linewidth',1.5); xlabel('frequency (Hz)'); ylabel('Amplitude');
axis([0,inf,0,inf]);

频谱泄露的解决方法1:改变采样时长


既然频谱泄露是由于采样时长不是信号整数倍周期导致的,那么只要将采样时长变成各主频信号周期的最大公倍数就可以了。比如我们将采样时长缩短为0.2s,此时做FFT得到下图,虽然采样时长变短导致频率分辨率降低,但是频谱泄露的问题解决了。
在这里插入图片描述
代码:

%% 抑制频谱泄露——方案一:采样到两个信号成分的整数倍周期
% 改变采样时长,使采样时长为各主频信号周期的整数倍
f0 = 50; f1 = 65; fs = 500; T = 0.2;
S = @(t) cos(2*pi*f0.*t) + cos(2*pi*f1.*t);
t = 0:1/fs:T-1/fs;
y = S(t);

nsamps = numel(y);
y_spec = fft(y);
y_spec = fftshift(y_spec)./nsamps;
f = ((0:nsamps-1)-nsamps/2)/T; 

figure(3);
subplot(2,1,1);plot(y,'linewidth',1.5); xlabel('time (s)'); ylabel('y');title('采样到整数倍周期');
subplot(2,1,2);plot(f, abs(y_spec),'linewidth',1.5); xlabel('frequency (Hz)'); ylabel('Amplitude');
axis([0,inf,0,inf]);

频谱泄露的解决方法2:加窗


对时域信号采样实际上默认对原信号增加了矩形窗,但是矩形窗的旁瓣能量较高,因此在采样时长不是信号整数倍周期时频谱泄露的效应比较明显。为了降低频谱泄露对主频率的影响,可以通过加入各种窗口的方法抑制旁瓣能量,常用的窗口包括汉宁窗、哈明窗、布莱克曼窗等,这里以哈明窗为例展示加窗作用。可以看到相比原先存在频谱泄露的频谱,加窗以后频谱能量更加集中了。(比如看75Hz处的幅度,加窗之后明显比没加窗时幅度低)
在这里插入图片描述
代码:

%% 抑制频谱泄露——方案二:加窗抑制频谱泄露
f0 = 50; f1 = 65; fs = 500; T = 0.3;
S = @(t) cos(2*pi*f0.*t) + cos(2*pi*f1.*t);
t = 0:1/fs:T-1/fs;
y = S(t);

w = hamming(fs * T);
y_hamm = y.*w.';

nsamps = numel(y_hamm);
y_spec = fft(y_hamm);
y_spec = fftshift(y_spec)./nsamps;
f = ((0:nsamps-1)-nsamps/2)/T; 

figure(4);
subplot(2,1,1);plot(y_hamm,'linewidth',1.5); xlabel('time (s)'); ylabel('y');title('hamming窗');
subplot(2,1,2);plot(f, abs(y_spec),'linewidth',1.5); xlabel('frequency (Hz)'); ylabel('Amplitude');
axis([0,inf,0,inf]);

更加直观地理解加窗的作用,可以认为加窗后,强行将阶段信号两端的信号幅度收敛到了0附近,增加了延拓后信号在衔接处的连贯性。 比如原先没加窗的信号是这样的:
在这里插入图片描述
在这里插入图片描述
加入窗函数之后变成:
在这里插入图片描述
强行使得信号连接上,从而降低了阶段处非整数倍信号周期的影响。

栅栏效应——补零


我们注意到加窗后得到的频谱依然无法得到 65 H z 65Hz 65Hz的频率分量。这是因为信号的采样时长为0.3s,因此物理分辨率是 f s / N = 1 / T = 3.333 H z f_s/N=1/T=3.333Hz fs/N=1/T=3.333Hz。此时65Hz不是整数倍的分辨率,因此无法在频谱上体现出来。这就是栅栏效应,解决该问题的常用方法为时域补零。(注意:时域补零不会增加信号的物理分辨率,但是会增加信号的计算分辨率,信号的物理分辨率只和采样时长有关)我个人觉得这里可以理解为DFT是对DTFT频域上的采样,补零相当于增加了采样点,如果DTFT本身就无法将两个频率成分区别开的话,那么频域上插再多的值(时域上补再多的零)也没有用。这里我们知道两个信号之间的频差为15Hz,大于3.33Hz所以理论上可以在频谱上将两个频点区分开,但是由于采样点数的限制,65Hz落在了两个采样的频点之间,为了将65Hz的频率成分提取出来,这里可以进行补零操作。结果如下:

在这里插入图片描述
代码:

%% 出现了栅栏效应——补零
f0 = 50; f1 = 65; fs = 500; T = 0.3;
S = @(t) cos(2*pi*f0.*t) + cos(2*pi*f1.*t);
t = 0:1/fs:T-1/fs;
y = S(t);

w = hamming(fs * T);
y_hamm = y.*w.';
y_hamm_add_0 = [y_hamm, zeros(1,fs*T/3)];

nsamps = numel(y_hamm_add_0);
y_spec = fft(y_hamm_add_0);
y_spec = fftshift(y_spec)./nsamps;
f = ((0:nsamps-1)-nsamps/2)/T; 

figure(5);
subplot(2,1,1);plot(y_hamm_add_0,'linewidth',1.5); xlabel('time (s)'); ylabel('y');title('补零');
subplot(2,1,2);plot(f, abs(y_spec),'linewidth',1.5); xlabel('frequency (Hz)'); ylabel('Amplitude');
axis([0,inf,0,inf]);

通过在时域信号上补零,65Hz的频率成分被分辨出来了。但是这里依然存在频谱泄露。说明补零只能缓解栅栏效应,但是不能解决频谱泄露的问题。

栅栏效应和频谱泄露的关系


(这是我的猜测,如果有不对的地方,还望大佬纠正。)如果不考虑补0操作,那么栅栏效应一定会导致频谱泄露,反之也成立 证明:

设信号对应的角频率为 ω 0 \omega_0 ω0,采样点数为 N N N。那么当信号频率不是频率分辨率的整数倍时,会产生栅栏效应,此时可得
ω 0 ≠ k 2 π N \omega_0 \neq k\frac{2\pi}{N} ω0=kN2π
ω 0 = 2 π f 0 / f s \omega_0=2\pi f_0/f_s ω0=2πf0/fs可得
2 π f 0 / f s ≠ k 2 π N 2\pi f_0/f_s \neq k\frac{2\pi}{N} 2πf0/fs=kN2π
f s f_s fs为采样频率, f 0 f_0 f0为原信号频率

f 0 ≠ k f s N f_0 \neq k\frac{f_s}{N} f0=kNfs

T T T为采样时长, T 0 T_0 T0为原始信号周期
1 T 0 ≠ k 1 T \frac{1}{T_0} \neq k\frac{1}{T} T01=kT1
T ≠ k T 0 T \neq kT_0 T=kT0

所以可得采样时长不是信号周期的整数倍,所以同样会产生频谱泄露,反之也成立。其中栅栏效应可通过补零的操作缓解,而频谱泄露可通过改变采样时长或者加窗缓解。

总结

解决频谱泄露和栅栏效应最有效的方法是改变采样时长,使得采样时长为整数倍信号周期,但是这在实际过程中常常难以实现,所以需要用加窗和补零的操作分别缓解频谱泄露和栅栏效应。

  • 22
    点赞
  • 170
    收藏
    觉得还不错? 一键收藏
  • 4
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值