最近做設(shè)計(jì)，要對(duì)采集到的數(shù)據(jù)進(jìn)行頻譜分析。剛開始將所有采樣數(shù)據(jù)全部送入Matlab進(jìn)行分析，效果總是很不理想。翻閱課本、搜查網(wǎng)頁，總結(jié)出使用Matlab對(duì)采樣數(shù)據(jù)進(jìn)行頻譜分析的理論和方法，特整理于此，和大家分享。。。

1、采樣數(shù)據(jù)導(dǎo)入Matlab

采樣數(shù)據(jù)的導(dǎo)入至少有三種方法。
第一就是手動(dòng)將數(shù)據(jù)整理成Matlab支持的格式，這種方法僅適用于數(shù)據(jù)量比較小的采樣。
第二種方法是使用Matlab的可視化交互操作，具體操作步驟為：File --> Import Data，然后在彈出的對(duì)話框中找到保存采樣數(shù)據(jù)的文件，根據(jù)提示一步一步即可將數(shù)據(jù)導(dǎo)入。這種方法適合于數(shù)據(jù)量較大，但又不是太大的數(shù)據(jù)。據(jù)本人經(jīng)驗(yàn)，當(dāng)數(shù)據(jù)大于15萬對(duì)之后，讀入速度就會(huì)顯著變慢，出現(xiàn)假死而失敗。
第三種方法，使用文件讀入命令。數(shù)據(jù)文件讀入命令有textread、fscanf、load等，如果采樣數(shù)據(jù)保存在txt文件中，則推薦使用 textread命令。如 [a,b]=textread('data.txt','%f%*f%f'); 這條命令將data.txt中保存的數(shù)據(jù)三個(gè)三個(gè)分組，將每組的第一個(gè)數(shù)據(jù)送給列向量a，第三個(gè)數(shù)送給列向量b，第二個(gè)數(shù)據(jù)丟棄。命令類似于C語言，詳細(xì)可查看其幫助文件。文件讀入命令錄入采樣數(shù)據(jù)可以處理任意大小的數(shù)據(jù)量，且錄入速度相當(dāng)快，一百多萬的數(shù)據(jù)不到20秒即可錄入。強(qiáng)烈推薦！

2、對(duì)采樣數(shù)據(jù)進(jìn)行頻譜分析

頻譜分析自然要使用快速傅里葉變換FFT了，對(duì)應(yīng)的命令即 fft ，簡(jiǎn)單使用方法為：Y=fft(b,N)，其中b即是采樣數(shù)據(jù)，N為fft數(shù)據(jù)采樣個(gè)數(shù)。一般不指定N，即簡(jiǎn)化為Y=fft(b)。Y即為FFT變換后得到的結(jié)果，與b的元素?cái)?shù)相等，為復(fù)數(shù)。以頻率為橫坐標(biāo)，Y數(shù)組每個(gè)元素的幅值為縱坐標(biāo)，畫圖即得數(shù)據(jù)b的幅頻特性；以頻率為橫坐標(biāo)，Y數(shù)組每個(gè)元素的角度為縱坐標(biāo)，畫圖即得數(shù)據(jù)b的相頻特性。典型頻譜分析M程序舉例如下：

clc
fs=100;
t=[0:1/fs:100];
N=length(t)-1;%減1使N為偶數(shù)
%頻率分辨率F=1/t=fs/N
p=1.3*sin(0.48*2*pi*t)+2.1*sin(0.52*2*pi*t)+1.1*sin(0.53*2*pi*t)...
+0.5*sin(1.8*2*pi*t)+0.9*sin(2.2*2*pi*t);
%上面模擬對(duì)信號(hào)進(jìn)行采樣，得到采樣數(shù)據(jù)p，下面對(duì)p進(jìn)行頻譜分析

figure(1)
plot(t,p);
grid on
title('信號(hào) p(t)');
xlabel('t')
ylabel('p')

Y=fft(p);
magY=abs(Y(1:1:N/2))*2/N;
f=(0:N/2-1)'*fs/N;
figure(2)
%plot(f,magY);
h=stem(f,magY,'fill','--');
set(h,'MarkerEdgeColor','red','Marker','*')
grid on
title('頻譜圖 （理想值：[0.48Hz,1.3]、[0.52Hz,2.1]、[0.53Hz,1.1]、[1.8Hz,0.5]、[2.2Hz,0.9]） ');
xlabel('f (Hz)')
ylabel('幅值')

對(duì)于現(xiàn)實(shí)中的情況，采樣頻率fs一般都是由采樣儀器決定的，即fs為一個(gè)給定的常數(shù)；另一方面，為了獲得一定精度的頻譜，對(duì)頻率分辨率F有一個(gè)人為的規(guī)定，一般要求F<0.01，即采樣時(shí)間ts>100秒；由采樣時(shí)間ts和采樣頻率fs即可決定采樣數(shù)據(jù)量，即采樣總點(diǎn)數(shù)N=fs*ts。這就從理論上對(duì)采樣時(shí)間ts和采樣總點(diǎn)數(shù)N提出了要求，以保證頻譜分析的精準(zhǔn)度。

3、數(shù)據(jù)長(zhǎng)度的選擇

頻率分辨率F，顧名思義就是頻譜中能夠區(qū)分出的最小頻率刻度。如F=0.01，則頻譜圖中橫坐標(biāo)頻率的最小刻度為0.01，即0.02Hz和 0.03Hz是沒有準(zhǔn)確數(shù)據(jù)的，但Matlab在畫圖時(shí)對(duì)其進(jìn)行了插值，故而plot作圖時(shí)看到的頻譜是連續(xù)的。但用stem來作圖就可以看出頻率是離散的，stem對(duì)了解F的含義非常有幫助。
由此，我們可以進(jìn)一步思考。如果信號(hào)所包含的頻率分量不是F的整數(shù)倍，那么這個(gè)頻率分量就不會(huì)得到正確的反映。如信號(hào)包含1.13Hz頻率分量，而 F=1/ts=fs/N=0.02，則1.13/0.02=56.5，不等于整數(shù)，即在頻譜圖中找不到準(zhǔn)確的刻度，而只能在第56和57個(gè)頻率刻度上分開顯示其幅值，這自然就不準(zhǔn)確了。
因此，請(qǐng)大家在頻譜分析時(shí)一定要使F能夠被頻率精度整除。如要求頻率精確度為0.01，則F最大為0.01，也可取值為0.02、0.05、0.001等數(shù)據(jù)，使0.01/F=整數(shù)。而F僅僅由采樣時(shí)間ts（也稱數(shù)據(jù)長(zhǎng)度）決定，因此一定要選擇好ts，且要首先確定ts的值。
作為驗(yàn)證，對(duì)上面的程序做一個(gè)修改：將t=[0:1/fs:100];改為t=[0:1/fs:83];即ts由100改為83，則F=1/ts由0.01變?yōu)?.012。二者分別作出頻譜圖對(duì)比如下：

上圖1 頻譜圖：ts=100s，F(xiàn)=1/ts=0.01

上圖2 頻譜圖：ts=83s，F(xiàn)=1/ts=0.012

對(duì)比上面兩個(gè)圖即可發(fā)現(xiàn)，圖2中由于f/F不是整數(shù)，在橫坐標(biāo)中找不到對(duì)應(yīng)的刻度，從而使得各個(gè)頻率的幅值泄漏到了其他頻率。

總結(jié)上面的結(jié)論，在保證采樣定理所要求的二倍頻的前提下，并不是采樣頻率fs或采樣點(diǎn)數(shù)N越大越好，而是要控制好數(shù)據(jù)長(zhǎng)度ts，使頻率分辨率F滿足頻率精度。