分析时域窗长度和FFT计算点数对频率分辨率和栅栏效应的影响

目录

 

频率分辨率

栅栏效应

频谱泄漏

实验结果：

窗长度改变：

改变fft计算点数

分析：

代码：

---

频率分辨率

频率分辨率是指将两个相邻谱峰分开的能力。在实际应用中是指分辨两个不同频率信号的最小间隔。果采样频率为fs，采样时间间隔为t，采样点数为N，采样时间为t（完成一组样本的采集所需要的时间），则频率分辨率Δf=fs/N，研究数字频谱最有效方法通常是离散傅里叶变换。

栅栏效应

栅栏效应是因为DFT计算的频谱被限制在基频的整数倍而不可能将频谱视为一个连续函数而产生的。就一定意义而言，栅栏效应表现为用DFT计算整个频谱时，就好像通过一个“栅栏”来观看一个图景一样，只能在离散点的地方看到真实图景。

频谱泄漏

所谓频谱泄漏 ,就是信号频谱中各谱线之间相互影响 ,使测量结果偏离实际值 ,同时在谱线两侧其他频率点上出现一些幅值较小的假谱。

产生频谱泄露的原因

由于计算机只能处理有限长的数据，所以需要对采集的信号进行截断，相当于对原始信号做了加窗处理。对信号加窗就是对信号在时域上乘以一个窗函数，时域的乘积对应频域的卷积，而窗函数的频域包括主瓣和旁瓣，旁瓣造成了信号频谱的泄漏。频域泄漏不可避免，只能减小。

如何抑制

可以取更长的数据点，与原始数据越接近越好，但缺点就是运算量加大；

可以选择窗谱的旁瓣能量较小的窗函数。

典型的窗函数中，矩形窗的频率分辨率最高，旁瓣泄露最大。

实验结果：

窗长度改变：

N=180;%信号长度

fs=9;%采样频率

l=10;%窗长度为l

M1=300;%加窗信号fft计算点数

M=M1;

 

L=40

L=80

L=140

L=180

改变fft计算点数

%加窗信号fft计算点数

N=180;%信号长度

fs=9;%采样频率

l=180;%窗长度为l

M1=50;%加窗信号fft计算点数

M=M1; 

140

180

185

190

200

280

400

800

分析：

频率分辨率有两种说法：机械频率分辨率和物理频率分辨率。机械频率分辨率由采样频率和采样点所决定的，为fs/N，N为采样点数，fs为采样频率。物理频率分辨率指的是频谱上两个谱峰之间的间隔。比如说：当信号为两个频率相近的正弦信号的叠加的时候，频谱信号会有两个相近的谱峰，那么频谱分辨率大小就决定了两个谱峰是否可以被分辨出。

数字信号处理中的频谱分辨率为1/T，T为信号长度，而与采样点数无关。

当矩形窗的长度从0增大时，DTFT的频谱分辨率增大，DFT频谱的栅栏效应会减弱。

由于矩形窗函数的频谱函数有旁瓣值，所以加窗后的信号的频谱函数也会有频谱泄露，但随着窗的长度不断变大，频谱泄露的情况不断改善，但是不可能完全消除。

当fft计算点数小于采样点数时，加窗函数的DFT的频谱图都有失真，当计算点数大于采样点数且不断增大时，DFT频谱函数与DTFT函数相差不大，频率分辨率会逐渐稳定为fs/N，此时，也会存在栅栏效应，若想减弱栅栏效应，可以通过增加fft计算点数来使DFT频谱更加光滑，此时采样频率f(s)会随之成正比上升，又由于频率分辨率F=f(s)/N，频率分辨率不改变，也就是说，补零不改变频率分辨率。

但是若像提高频率分辨率，是不能只增大计算点数M的，因为增大M点数，相当于在信号序列后面补0，即相当于在DFT频域上进行插值，虽然最小间隔变小了，但只是把DFT频谱变得更加光滑了；所以，若想提高频谱分辨率，需要增大采样信号的长度。

代码：

%分析时域窗长度和FFT计算点数对频率分辨率和栅栏效应的影响。对含有两个间隔较小的频率分量的余弦序列，

%时域加窗后计算FFT，画出DTFT和DFT频谱图，对比分析不同的时域窗长度和FFT计算点数。

clear all;

% t=10;%信号长度

N=180;%信号长度

fs=9;%采样频率

l=180;%窗长度为l

M1=800;%加窗信号fft计算点数

M=M1;

% M=100;%原信号fft计算点数

n=(0:N-1)./fs;%离散时间序列

x=cos(0.5.*pi.*n)+cos(0.7.*pi.*n);

% x=10.*cos(1.3.*pi.*n);

%原始信号

figure,subplot(2,2,1),plot(n,x);title('原始信号时域函数')

%原信号DTFT变换

w=linspace(-2.*pi,2.*pi,1000);%抽样点为1000，两个周期内的dtft频域函数

martrix = w'.*n;

X=exp(-j.*(martrix))*x';%dtft变换

subplot(222),plot(w,X),title('原信号DTFT变换');

%原信号DFT变换,M为fft计算点数

% M=500;

% k=0:M;

% Xn=x*(exp(-j*2.*pi/M)).^(n'*k);

% subplot(223),plot(2.*pi.*k./M,Xn),title('原信号DFT变换');

Xn=fft(x,M);

subplot(223),stem(linspace(0,2.*pi,length(Xn)),fftshift(Xn),'MarkerSize',3),title('原信号DFT变换');

%窗函数

ht=ones(1,l);

tem=zeros(1,length(x)-l);

ht=[ht tem];

H=fft(ht,100);%窗的DFT变换

% figure,subplot(221),plot(0:length(ht)-1,ht,'r'),title('时域窗')

figure,subplot(221),stem(0:length(H)-1,fftshift(H),'MarkerSize',3),title('窗DFT频响函数')

%加窗

x1=x.*ht;

subplot(222),plot((0:length(x1)-1)./fs,x1),title('加窗后时域图')

%加窗信号DTFT变换

w1=linspace(-2.*pi,2.*pi,1000);%抽样点

martrix1 = w1'.*n;

X1=exp(-j.*(martrix1))*x1';

subplot(223),plot(w1,X1),title('加窗信号DTFT变换');

%加窗信号DFT变换

Xn1=fft(x1,M1);

subplot(224),stem(linspace(0,2.*pi,length(Xn1)),fftshift(Xn1),'MarkerSize',3),title('加窗信号DFT变换');