即以此功德,庄严佛净土。上报四重恩,下济三途苦。惟愿见闻者,悉发菩提心。在世富贵全,往生极乐国。
缘起:现实中周期信号有的很广泛,比如雷达和声纳的载波都是周期信号。对周期信号的数字采样是数字信号处理的首要环节,采样中要考虑采样频率Fs>2*B就可保证信号在频域不发生混叠;在频谱分析中要NFFT的点数是整周期采样点数N的整数倍时即不发生频谱泄露。这两个概念对于初学者来讲如果不按释迦摩尼佛讲佛经的方法:譬喻+观想来形象化的理解的话是很难指导实际工作的。但是网上很少有人透彻分析其细节和因果,给初学者造成很多困惑。因此为了帮助初学者理解这两个关键性概念,真正的学懂数字信号处理,并能够很有信心的在实践中用对数字信号处理,本人特此花时间写篇博客。其目的:一方面希望把初学者引入数字信号处理的大门,可以说搞定混叠和泄露就类似于佛家中的明心见性了,才可以说真正入门了;另一方面希望大家多分享自己的知识,学习佛法,自利利他,断贪嗔痴,秀戒定慧,得大智慧,当有助于社会的有用之才!南无阿弥陀佛!
(一) 混叠和泄露的概念
混叠 从取样定理知道,在时域欠取样情况下,由于出现频谱混叠无法恢复原信号频谱,因而人们不能从时域样点准确的重建远连续信号。同理,在频率欠取样情况,由于出现信号波形混叠无法恢复原频谱所对应的信号,因而人们不能从频率样点重建原连续频谱。对周期信号而言,混叠所造成的影响结论也一样,只是这是的频谱是离散的而且具有谐波性。根据离散时间傅立叶级数是一个有限级数这一特点,从理论上阐明了离散时间周期信号的频谱是一个周期性且只具有有限数字频率分量的离散频谱。因此对那些具有无限频谱分量的连续时间周期信号如矩形和三角形等脉冲串,不然无法准确地从有限样点求得原始周期信号的频谱,而只能通过恰当地提高取样率,增加样点数,来减少混叠对频谱分析所造成的影响。
泄漏 离散时间周期信号除了因取样速率低于取样定理要求,使频谱分析出现混叠误差以外,还由于截取波形的时间长度不恰当造成泄漏误差。这情况在实际中往往由于事先不知道信号的确切周期所致。
(二)验证混叠和泄露
从(一)中往往不能真正学会怎样排除混叠和泄露的方法。这就跟整天抱着佛经看而不实践是一样的,譬如说食不能饱,画饼充饥,纸上谈兵一样。下面在MATLAB中对此两个概念分别做仿真验证,我的感觉是做仿真验证下就学会了,知识不是学来的,是实践来的。
% 时间间隔t=0.001s,那么f=1/0.001s = 1000Hz=1KHz;
% 所以:1KHz的采样率来近似模拟信号xa, 看看真空妙有,注意体会数学中的因果和真空妙有和释迦摩尼佛讲的一模一样;
% 实际长度为2s;
t = 0:0.001:2;
% 因为ω=2πf,所以y = sin(ωt) = sin(2πft);
% 周期T=2*pi/ω=1/f,故:f=5时有T=0.2s,说明xa是数字周期T=0.2s的周期信号;
xa = sin(2*pi*5*t);
plot(t,xa) % 画出“模拟信号”, 该信号的带宽B = 5Hz
hold on
%设置标题:
%www.fxian.org换为自己想要的名称即可
set(figure(1),'NumberTitle','off','Name','www.fxian.org');
fs = 15;
% 在信号总持续时间Tp=2秒内以fs为采样间隔一共采了N=Tp*Fs=2*7.5=15个点;
% 1/fs=1/15, 那么在[0,2]之间的样本点个数:n=(2-0)/(1/15)=30
ts = 0:1/fs:2;
% numel(ts)=16是因为采样点从1开始到15为15个采样点,起始0点不算采样点,或者理解为:16个端点形成15个采样间隔
% 1) 不混叠的情况
% 用 fs = 15Hz 的采样率对模拟信号xa进行采样,fs=15表示1秒内采15个样本点;
% 在[0,2]内共采集30+1个点:
xs1 = sin(2*pi*5*ts);
plot(ts,xs1,'ro-')
% 2) 混叠的情况
%用 fs = 7.5Hz 的采样率,fh=5Hz 混叠到 (-fs/2, fs/2)内的频率为 5 - fs = -2.5 Hz
fs = 7.5;
%混叠为-2.5Hz的正弦信号,注意前面的负号代表相位的反向
ts = 0:1/fs:2;
xs2 = sin(2*pi*5*ts);
%使用绿色的虚线绘制:
plot(ts,xs2,'--g')
hold off;
- 下面看看几种情况:
% ①NFFT=N不会产生频谱泄露,因为是NFFT是N的1倍是正周期采样;
% 用sptool做fft时一定要切记NFFT=N才出正确结果
N = 2*fs;
NFFT = N;
% 注意频率轴的刻度怎么标,由于没人指导,我以前在这个地方困扰了好几天。
f = [0:NFFT-1]*fs/NFFT;
plot(f, 2/N*abs(fft(xs2, NFFT)));
% ② NFFT=10*N 不会产生频谱泄露
N = 2*fs;
NFFT = 10*N;
f = [0:NFFT-1]*fs/NFFT;
% 注意matlab中DFT矩阵不是酉矩阵,故f>0处的信号幅度要除以N/2进行修正后
% 才是原始信号的幅度,注意不是除以NFFT/2;
plot(f, 2/N*abs(fft(xs2, NFFT)));
% ③ NFFT=16 会产生频谱泄露,因为 Fs/NFFT 除不尽,故可知 FFT 的点数不一定选 2 的整数次幂。
NFFT = 16;
f = [0:NFFT-1]*fs/NFFT;
plot(f, 2/N*abs(fft(xs2, NFFT)));