問題簡介:
對於頻率為fs的正弦序列,它的頻譜應該只是在fs處有離散譜。但是,在實際利用DFT求它的頻譜時,對時域做了截斷,結果使信號的頻譜不只是在fs處有離散譜,而是在以fs為中心的頻帶范圍內都有譜線出現,它們可以理解為是從fs頻率上“泄漏”出去的,這種現象稱 為頻譜“泄漏”。
不發生泄露的條件:F_c=(mF_s)/N (m為整數,F_c為信號頻率,F_s為采樣頻率)
現象分析:
以7K采樣為例:
DFT本質是以一個離散的窗函數對信號進行取樣,也就是在時域相乘的過程。時域相乘又等價於頻域的卷積,並對頻域的幅值產生影響,在這里只需要討論頻譜形態而不必討論幅值,由於本次采樣周期為500點,我們畫出一個500點矩形序列的頻譜圖
可以看出每個頻域的每個小峰寬度都為14Hz,理論上說他們的間隔應為
Fs/N=7000/500=14
與理論相符合
下一步是與信號的頻譜進行卷積,由於原信號的理論頻譜為在F_c及〖-F〗_c的沖激信號,我們只要將目前的頻譜向左向右移動F_c即可,為了方便,只研究右半塊。
時域采樣率為F_s那么在頻域中就以F_s/N=7000/500=14Hz,從零點開始等間隔取樣,而調制過的信號每個小峰的寬度也為F_s/N=7000/500=14Hz,主峰的寬度為小峰寬度的兩倍,也就是28Hz,此時在主峰中會以357*14=4998Hz及358*14=5012Hz進行采樣,並且在小峰中也會在峰體中采樣,從而顯出原信號不存在的頻率。
小峰寬度與頻域采樣寬度相等這就表明每個頻域取樣點可以正好落在每個小峰的零點上,而落在最高峰的中間,完美顯示出原信號的頻譜分布,此時采樣信號頻率應滿足:
F_c=(mF_s)/N (m為整數)
如果采樣頻率不滿足這一公式,那么頻率域采樣點就會處於各個旁瓣中,從而在頻域中表現出原信號不存在的頻率,這些旁瓣的幅值與原信號頻率距離成反比,越靠近信號頻率造成的諧波振幅越大。
改進措施
增大FFT變換的點數N。 通過增大N,一方面提高頻域分辨率,更大可能滿足F=m*Fs/N, 另一方面,壓縮采樣信號的瓣寬,降低泄露水平。
下圖為以不同采樣點數得出的頻譜:
可以明顯看出500點采樣比64點采樣效果好很多。
選用合適的窗函數。根據不同的需求來選擇不同特性的窗函數,主瓣寬但旁瓣衰減大的窗或是主瓣窄但旁瓣相對衰減小的窗。
選取7K, 64點采樣的頻域信號進行測試,選取常用的Hanning窗:
可以看出雖然加窗可以很好地抑制旁瓣,但對於主瓣的頻率卻不能很好地優化,我也試了一下其他的加窗函數,都不太理想,對於本題測試的窄帶信號,主瓣的信號尤為重要,只能采用增加取樣點的方法來凸顯主瓣信號。、
總結
- 要注意理論與實際的差異,並善於思考其中的原因
- 善於運用優質網絡資源
參考文檔:https://www.ni.com/white-paper/4844/zhs/
附c++采樣程序:
#include<stdio.h>
#include<math.h>
#include<iostream>
#include<fstream>
using namespace std;
void main()
{
double Fs_7 =7,Fs_10=10,Fs_11=11;//K
int i;
int const dn = 10;//延時濾波器長度
int const s_n =500 ; //序列長度
double xn[500];
for (i = 0; i<s_n; i++)
{
xn[i] =sin(2 * 3.1416 * 5 / Fs_7 * i);
}
ofstream SaveFile_a("xn_7.txt");
for (i = 0; i<s_n; i++)
SaveFile_a << " " << xn[i] << endl;
SaveFile_a.close();
for (i = 0; i<s_n; i++)
{
xn[i] = sin(2 * 3.1416 * 5 / Fs_10 * i);
}
ofstream SaveFile_a1("xn_10.txt");
for (i = 0; i<s_n; i++)
SaveFile_a1 << " " << xn[i] << endl;
SaveFile_a1.close();
for (i = 0; i<s_n; i++)
{
xn[i] = sin(2 * 3.1416 * 5 / Fs_11 * i);
}
ofstream SaveFile_a1("xn_11.txt");
for (i = 0; i<s_n; i++)
SaveFile_a1 << " " << xn[i] << endl;
SaveFile_a1.close();
}
Matlab程序
顯示采樣過信號的頻譜:
xn_f=fopen('xn_7.txt','r');
[xn,count]=fscanf(xn_f,'%f');
fclose(xn_f);
Fs=7000;
x=xn;
x=x.*hamming(64);%加窗
x = x(:,1);
x = x';
N = length(x);%
t = (0:N-1)/Fs;%
y = fft(x);
f = Fs/N*(0:round(N/2)-1);
stem(f,abs(y(1:round(N/2))));
xlabel('Frequency / (s)');ylabel('Amplitude');
title('7K²ÉÑùµÄƵÆ×');
grid;
顯示采樣矩形序列頻譜:
Fs=7000;
x=boxcar(500);%生成矩形序列
x = x(:,1);
x = x';
N = length(x);
y = fft(x,1000);
y=circshift(y,143,2);
y=abs(y)/10;
f1 = Fs/length(y)*(0:length(y)-1);
plot(f1,y);
xlabel('Frequency / (s)');ylabel('Amplitude');
title('ÐÅºÅµÄÆµÆ×');
grid;