基于matlab深入形象理解频率分辨率，补零，栅栏效应，频谱泄漏

写在前面

不知大家是否对补零的作用是什么，栅栏效应，频谱泄漏，采样频率，采样点数，频率分辨率，它们之间的联系与区别有所困惑呢？请耐心读完，并用matlab来帮助我们更好的理解。

频率分辨率

可以理解为在使用DFT时，在频率轴上的所能得到的最小频率间隔
它实际是作FFT时谱图中的两条相邻谱线之间的频率间隔，也有称作步长。f0=fs/N=1/NTs=1/T,其中N为采样点数，fs为采样频率，Ts为采样间隔。所以NTs就是采样前模拟信号的时间长度T，所以信号长度越长，频率分辨率越好。是不是采样点数越多，频率分辨力提高了呢?其实不是的，因为一段数据拿来就确定了时间T，注意:f0=1/T，而T=NTs，增加N必然减小Ts ，因此，增加N时f0是不变的。只有增加点数的同时导致增加了数据长度T才能使分辨率越好。还有容易搞混的一点，我们在做DFT时，常常在有效数据后面补零达到对频谱做某种改善的目的，我们常常认为这是增加了N，从而使频率分辨率变好了，其实不是这样的，补零并没有增加有效数据的长度，仍然为T。也可以从信息论角度来分析，并未增加数据的信息熵，在另一域自然也不会多出信息熵，不仅能量守恒，信息熵也守恒。
补零好处
1)使数据N为2的整次幂，便于使用FFT。
2)补零后，其实是对DFT结果做了插值，克服“栅栏”效应，使谱外观平滑化;我把“栅栏”效应形象理解为，就像站在栅栏旁边透过栅栏看外面风景，肯定有被栅栏挡住比较多风景，此时就可能漏掉较大频域分量，但是补零以后，相当于你站远了，改变了栅栏密度，风景就看的越来越清楚了。
栅栏效应：
在DFT过程中，最后对信号的频谱采样，有些重要的峰值会被忽略，fn=(n-1)fs/N;dF=fs/N;N增大，dF减小，使X（k）做插值处理，外观变得平滑。

3)由于对时域数据的截短必然造成频谱泄露，因此在频谱中可能出现难以辨认的谱峰，补零在一定程度上能消除这种现象，但也可能加重。
频谱泄漏：
信号为无限长序列，运算需要截取其中一部分（截断），于是需要加窗函数，加了窗函数相当于时域相乘，于是相当于频域卷积，于是频谱中除了本来该有的主瓣之外，还会出现本不该有的旁瓣，这就是频谱泄露！为了减弱频谱泄露，可以采用加权的窗函数，加权的窗函数包括平顶窗、汉宁窗、高斯窗等等。而未加权的矩形窗泄露最为严重。
我们考虑窗函数主要是以下几点:
主瓣宽度B最小(相当于矩形窗时的4π/N，频域两个过零点间的宽度)。
最大边瓣峰值A最小(这样旁瓣泄露小，一些高频分量损失少了)。3.边瓣谱峰渐近衰减速度D最大(同样是减少旁瓣泄露)。
　对于周期信号，整周期截断是不发生频谱泄露的充分且必要条件，抑制频谱泄露应该从源头抓起，尽可能进行整周期截断。
N点的选取
1)由采样定理:fs>=2fh，
2)频率分辨率:f0=fs/N，所以一般情况给定了fh和f0时也就限制了N范围:N>=fs/f0。
频率/采样点数/谱线数的设置要点
1/fs代表采样周期，是时间域上两个相邻离散数据之间的时间差。
fs/N用在频率域，只在DFT以后的谱图中使用，**FFT在matlab中，横轴不是代表真实的频率，而是点数，要通过变换将点数对应频率fn=(n-1)f0=(n-1)fs/N=(n-1)/T;单位Hz,转换为角频率wn=2πfn.**注意区别fourier，它横轴是
真实的频率。
最高分析频率:Fm指需要分析的最高频率，也是经过抗混滤波后的信号最高频率。根据采样定理，Fm与采样频率Fs之间的关系一般为:Fs=2.56Fm;
2)采样点数N与谱线数M有如下的关系:
N=2.56M 其中谱线数M与频率分辨率ΔF及最高分析频率Fm有如下的关系:ΔF=Fm/M 即:M=Fm/ΔF 所以:N=2.56Fm/ΔF

基于matlab实例分析

一个信号，这个信号中仅包含两个正（余）弦波，一个是1.05MHz ，一个是 1MHz ，即 f(t)=cos(2π1050000t)+cos(2π1000000t)。设定采样频率为Fs=100Mhz，如果采 1000个点，那么时域信号的时长就有 10微秒
。
NO.1

clear;clc
close all %频率分辨率为1/NTs=fs/N=100e6/1000=100kHz>50kHz 所以分不出来
%% FFT
Fs = 100e6;            % 采样频率 Hz
T = 1/Fs;          % 采样周期 s
L0 = 1000;                  % 信号长度
L = 1000;                   % 数据长度
t0 = (0:L0-1)*T;            % 信号时间序列
x = cos(2*pi*1e6*t0) + cos(2*pi*1.05e6*t0); % 信号函数
t = (0:L-1)*T;              % 数据时间序列
%% Plot
figure(1)
plot(t*1e6,x,'b-','linewidth',1.5)
title ('\fontsize{10}\fontname{Times New Roman}Time domain signal')
xlabel('\fontsize{10}\fontname{Times New Roman}\it t /\rm \mus')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itt\rm)')
grid on;
axis([0 10 -2 2])
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')
%% FFT
Y = fft(x);                   % 快速傅里叶变换
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1 双边普的幅度两倍，对折，移到零点，fftshift是双边频谱移到零点偶对称。
P2 = abs(Y/L0);   %fft除以L0才是真实的fourier系数
P1 = P2(1:L/2+1); %第一个点到
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;% Rfft
figure(2)
plot(f,P1,'r-','Marker','.','markersize',10,'linewidth',1.5)
axis([0.5e6 1.5e6 0 1.5])
title('\fontsize{10}\fontname{Times New Roman}Power Spectrum')
xlabel('\fontsize{10}\fontname{Times New Roman}\it f /\rm Hz')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itf\rm)')
grid on;
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')
% figure(3)
% f = Fs*(0:(L-1))/L;% Rfft
% plot(f,P2,'r-','Marker','.','markersize',10,'linewidth',1.5)  %看下双边普
% figure(4)
% plot(f,abs(fftshift(Y)),'r-','Marker','.','markersize',10,'linewidth',1.5)  %看下fftshift双边普

频谱点稀疏，在1MHz附近根本无法将1 MHz 和1.05 MHz 的两个频率分开。
NO.2
在1基础上补0 可以看出时时域上前面与1相同补零6000个，在频谱上显得更加细致，
但补零操作增加了频域的插值点数，让频域曲线看起来更加光滑，即减缓了栅栏效应
频率分辨率不改变Fs/Nfft=100kHz>50kHz 增加了FFT频率分辨率使得频谱细化

clear;clc
close all
%% FFT
Fs = 100e6;            % 采样频率 Hz
T = 1/Fs;          % 采样周期 s
L0 = 1000;                  % 信号长度
L = 7000;                   % 数据长度
t0 = (0:L0-1)*T;            % 信号时间序列
x = cos(2*pi*1e6*t0) + cos(2*pi*1.05e6*t0); % 信号函数
t = (0:L-1)*T;              % 数据时间序列
y = zeros(1,L);
y(1:L0) = x;
%% Plot
figure(1)
plot(t*1e6,y,'b-','linewidth',1.5)
title('\fontsize{10}\fontname{Times New Roman}Time domain signal with Zero Padding')
xlabel('\fontsize{10}\fontname{Times New Roman}\it t /\rm \mus')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itt\rm)')
grid on;
axis([0 70 -2 2])
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')
%% FFT
Y = fft(y);                   % 快速傅里叶变换
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/L0);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;% Rfft
figure(2)
plot(f, P1,'r-','Marker','.','markersize',10,'linewidth',1.5)
axis([0.5e6 1.5e6 0 1.5])
title('\fontsize{10}\fontname{Times New Roman}Power Spectrum')
xlabel('\fontsize{10}\fontname{Times New Roman}\it f /\rm Hz')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itf\rm)')
grid on;
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')

频谱点密集了不少，但是在 1MHz 附近依然无法将 1.05MHz和 1MHz的两个频率成分分开。
可以看出，波形分辨率只与原始数据的时长T有关
也就是说，补零操作增加了频域的插值点数，让频域曲线看起来更加光滑，也就是增加了FFT频率分辨率.
NO.3
这里不补零而是真实的增加信号长度频率分辨率=14kHz<50kHz
采样频率对1MHz来说是整周期取样所以无频谱泄漏但1.05MHz有
但是其周边的点上却都有不小的幅值，这就是所谓的频谱泄露，因为数据点的个数影响，1mhz处有谱线存在，但在1.05处没有谱线存在

clear;clc
close all%% FFT
Fs = 100e6;            % 采样频率 Hz
T = 1/Fs;          % 采样周期 s
L0 = 7000;                  % 信号长度
L = 7000;                   % 数据长度
t0 = (0:L0-1)*T;            % 信号时间序列
x = cos(2*pi*1e6*t0) + cos(2*pi*1.05e6*t0); % 信号函数
t = (0:L-1)*T;              % 数据时间序列
%% Plot
figure(1)
plot(t*1e6,x,'b-','linewidth',1.5)
title('\fontsize{10}\fontname{Times New Roman}Time domain signal')
xlabel('\fontsize{10}\fontname{Times New Roman}\it t /\rm \mus')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itt\rm)')
grid on;
axis([0 70 -2 2])
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')
%% FFT
Y = fft(x);                   % 快速傅里叶变换
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/L0);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs/2*linspace(0,1,L/2+1);% Rfft
figure(2)
plot(f, P1,'r-','Marker','.','markersize',10,'linewidth',1.5)
axis([0.5e6 1.5e6 0 1.5])
title('\fontsize{10}\fontname{Times New Roman}Power Spectrum')
xlabel('\fontsize{10}\fontname{Times New Roman}\it f /\rm Hz')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itf\rm)')
grid on;
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')

NO.4
为了解决上个问题，可以设法使得谱线同时经过 1.05MHz 和 1MHz这两个频率点，找到他们的公约数。
如果原始数据不变，在后面再补充 1000 个零点：

clear;clc
close all  %在3基础上补零看看对频谱泄漏的影响
%增加了FFT频率分辨率使得频谱细化
%FFT分辨率就是 fs/8000=12.5khz ，是这两个频率的公约数
%所以谱线同时经过 1 和 1.05 这两个频率点。 减缓了频谱泄漏 但其余部分任有少量泄漏，会有一些旁瓣出现，这是因为补零影响了原始信号
%所以将L0=8000,试试 这样就不存在补零带来的误差了。
%补零仅是减小了频域采样的间隔。这样有利于克服由于栅栏效应带来的有些频谱泄露的问题 但也可能加重
%% FFT
Fs = 100e6;                % 采样频率 Hz
T = 1/Fs;          % 采样周期 s
L0 = input('please 7000 or 8000=');                  % 信号长度
L = 8000;                   % 数据长度  把它改成16000试试 频谱更清晰
t0 = (0:L0-1)*T;            % 信号时间序列
x = cos(2*pi*1e6*t0) + cos(2*pi*1.05e6*t0); % 信号函数
t = (0:L-1)*T;              % 数据时间序列
y = zeros(1,L);
y(1:L0) = x;
%% Plotfigure(1)
plot(t*1e6,y,'b-','linewidth',1.5)
title('\fontsize{10}\fontname{Times New Roman}Time domain signal')
xlabel('\fontsize{10}\fontname{Times New Roman}\it t /\rm \mus')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itt\rm)')
grid on;
axis([0 80 -2 2])
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')
%% FFT
Y = fft(y);                   % 快速傅里叶变换
% 计算双侧频谱 P2。然后基于 P2 和偶数信号长度 L 计算单侧频谱 P1。
P2 = abs(Y/L0);  %%%%%%%除以信号长度才是真正幅度  ？ 谱密度 乘以1/L0   看一下FT和DFT公式，用fft后乘以dt才是真实
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs/2*linspace(0,1,L/2+1);% Rfft
figure(2)
plot(f, P1,'r-','Marker','.','markersize',10,'linewidth',1.5)
axis([0.5e6 1.5e6 0 1.5])
title('\fontsize{10}\fontname{Times New Roman}Power Spectrum')
xlabel('\fontsize{10}\fontname{Times New Roman}\it f /\rm Hz')
ylabel('\fontsize{10}\fontname{Times New Roman}\it y\rm(\itf\rm)')
grid on;
set(gca,'FontSize', 10 ,'FontName', 'Times New Roman')
set(gcf,'unit','centimeters','position',[15 10 13.53 9.03],'color','white')

L0=7000，补零1000个点
上图会有一些旁瓣出现，这是因为补零影响了原始信号，如果，直接采8000个点作为原始数据，即将程序中的L0改为8000，那么有：

写在后面

看完这些概念和例子你是否对它们有了更清晰的理解了？
补零到底改不改变频率分辨率？
补零对于频域有什么影响？
如何增加频率分辨率？
如何减轻或解决频谱泄露的问题？
关于上述问题我们首先要先定义一个采样频率Fs,它要满足什么条件？
信号N的增大影响了频率分辨率吗？
欢迎在评论区留言，测一测你是否真的掌握了那些东西。