选择信号的采样频率和信号长度的技巧

在实际信号分析中经常会遇到要分辨出频率间隔为△f的两个分量，在这种情形中如何选择采样频率和信号的长度呢？
设有一个信号x（t）由三个正弦信号组成，其频率分别是f1=1Hz，f2=2.5Hz，f3=3Hz,即
x(t)=sin(2π*f1*t)＋sin(2π*f2*t)＋sin(2π*f3*t)
下面介绍如何选择采样频率fs和信号长度N。
因为信号的最高频率fc为3Hz，故按采样定理fs≥2fc=6，选择fs=10Hz。由频域分辨分析可知，若要区分1kHz和2.5kHz的频率分量，则按式N≥[fs/△fmin]，最小采样长度N1必须满足
N1≥[fs/△fmin]=10/2.5-1=6.6
若要区分2.5kHz和3kHz的频率分量，则最小采样长度N2必须为
N2≥[fs/△fmin]=10/3-2.5=20
因此，为了能区分各频率的峰值，信号最小长度为20。

案例1：某信号由3个正弦信号组成，频率分别为1Hz、2.5Hz、3Hz，采样频率为10Hz。分别以数据长度N=20、40、128来分析该信号。程序如下：

clear all; clc; close all;

M=256; fs=10;                                   % 设置数据长度M和采样频率fs
f1=1; f2=2.5; f3=3;                             % 设置3个正弦信号的频率
t=(0:M-1)/fs;                                   % 设置时间序列
x=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t); % 计算出信号波形

X1=fft(x,20);                                   % FFT变换
X2=fft(x,40);
X3=fft(x,128);
freq1=(0:10)*fs/20;                              % 计算3个信号在频域的频率刻度
freq2=(0:20)*fs/40;
freq3=(0:64)*fs/128;
% 作图
plot(freq1,abs(X1(1:11)),'r--');hold on;
plot(freq2,abs(X2(1:21)),'k-');
plot(freq3,abs(X3(1:65)),'b-','linewidth',2);
legend('N=20','N=40','N=128');
title('不同N值的DFT变换');
xlabel('频率/Hz'); ylabel('幅值');
set(gcf,'color','w');

运行结果如下：

从图中可以看出，当N=20点时，虽然2.5Hz和3Hz这两个峰值大致能分开，但还是不太明显，可以认为是两个峰值，也可能被误认为有一个峰值在这两点之间。

当N=40时这两个峰值就十分明显了，因为N增加一倍后在这两点之间增加了一个谷值，从而突出了峰值。

而当N=128时峰值更明显了，但由于栅栏现象和矩形窗泄漏存在，3个正弦信号虽然输入时幅值相同，但从频域上反映出的幅值各不相同。

在以上例子中3个正弦信号的幅值是相同的，但大多数情形中相邻正弦信号幅值不相同，有时差值很大。
案例2：已知信号x(t)中50Hz频率分量的幅值为311，46Hz频率分量的幅值为1.55，采样频率fs=8kHz。要求46Hz信号的幅度分析精度不低于5%，试问：
①选择何种类型窗函数较合适？
②采样长度N应为多少？
③分析信号的实际频谱。
由于50Hz频率分量幅度远大于46Hz频率分量，需要防止50Hz频率分量对46Hz频率分量的泄漏，同时考虑46Hz频率分量分析精度的要求，允许50Hz频率分量的最大泄漏为
20lg|(1.55/311)*5%|=-80dB
以上介绍的窗函数中没有一个窗函数的第一旁瓣衰减能达到一80dB；但可以选择旁瓣衰减大、高频衰减速度快的窗函数，以满足实际衰减要求。根据下表可知，

选择布莱克曼窗比较合适。布莱克曼窗第3个旁瓣衰减为（-58-18×3）dB=-112dB。
此时，采样长度的选择不但要考虑窗函数的主瓣宽度，还要考虑旁瓣位置：
N≥fs/△fmin*(K+M)=[8000/(50-46)]*（3+3）=12000
式中：K为窗函数的主瓣宽度与矩形窗的主瓣宽度之比，M为旁瓣位置。按上面表中的信息矩形窗的主瓣宽度为4π/N，而布莱克曼窗的主瓣宽度为12π/N，所以K=3；又取布莱克曼窗第3个旁瓣，M=3。
程序如下：

clear all; clc; close all;

f1=50; a1=311.46;                        % 设置第1个分量的频率与幅值
f2=46; a2=1.57;                          % 设置第2个分量的频率与幅值
N=12000;                                 % 设置数据长度N
fs=8000;                                 % 设置采样频率fs
t=(0:N-1)/fs;                            % 设置时间刻度
x=a1*cos(2*pi*f1*t)+a2*cos(2*pi*f2*t);   % 设置信号
freq=(0:N/2)*fs/N;                       % 设置频率刻度
wind=blackman(N)';                       % 给出布莱克曼窗函数
X=fft(x.*wind);                          % FFT
plot(freq,20*log10(abs(X(1:N/2+1))),'k');% 作图
grid; xlim([0 100])
xlabel('频率/Hz'); ylabel('幅值/dB');
title('信号谱图');
set(gcf,'color','w');

运行结果如下：

 参考文献：MATLAB数字信号处理85个实用案例精讲——入门到进阶；宋知用（编著）