以下的matlab程序分別使用周期圖法、相關函數法以及AR譜方法計算信號的功率譜。
% power spectrum estimated
clear all;
clc;
close all;
Fs=1000; % 採樣頻率
nfft = 1024; % fft計算點數
%產生含有噪聲的序列
n=0:1/Fs:1;
xn=cos(2*pi*100*n)+3*cos(2*pi*200*n)+(2*randn(size(n)));
subplot(2,2,1);plot(xn);title('加噪信號');grid on
% 周期圖法
window=boxcar(length(xn)); %矩形窗
[psd1,f]=periodogram(xn,window,nfft,Fs); %直接法
psd1 = psd1 / max(psd1);
subplot(2,2,2);plot(f,10*log10(psd1+0.000001));title('周期圖法');grid on
% 自相關法
cxn=xcorr(xn,'unbiased'); %計算序列的自相關函數
CXk=fft(cxn,nfft);
psd2=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
psd2 = psd2 / max(psd2);
psd2=10*log10(psd2(index+1)+0.000001);
subplot(2,2,3);plot(k,psd2);title('自相關法');grid on
% AR譜
psd3 = pyulear(xn, Fs, nfft);
psd3=psd3/max(psd3);
psd3 = psd3 / max(psd3);
index=0:round(nfft/2-1);
psd3=10*log10(psd3(index+1)+0.000001);
subplot(2,2,4);plot(k, psd3);title('AR譜預計');grid on;
如今就此來說是關於功率譜的幾點理解:
功率譜的數據都是相對值,他無法給出信號的實際絕對幅值,一般只要看峰值之間的比值即可了,也能夠對數據歸一化
功率譜中的峰值代表的是信號中的周期成分,隱含的周期信號能量要比隨機信號大(這個一般能理解的)。
比方,上面xn是
xn=cos(2*pi*100*n)+3*cos(2*pi*200*n)+(2*randn(size(n)));
其周期包含200Hz及100Hz,200Hz的幅值比100Hz大。再看功率譜圖中。清晰的體現了信號中含有200Hz和100Hz的信號。
直接法又稱周期圖法,它是把隨機序列x(n)的N個觀測數據視為一能量有限的序列。直接計算x(n)的離散傅立葉變換。得X(k),然后再取其幅值的平方,並除以N,作為序列x(n)真實功率譜的預計
間接法先由序列x(n)預計出自相關函數R(n)。然后對R(n)進行傅立葉變換。便得到x(n)的功率譜預計
直接採用平均周期圖法求功率譜時,功率普形狀呈鋸齒形,譜峰點的准確位置不大好定。於是能夠採用其它的方法對譜進行平滑操作,平滑化,不過為了使圖形光滑,並不會使得波的本質受到歪曲和畸變。
反過來說,因為不純的東西去掉了。本質的東西必定會更加顯示出來!平滑化的程度能夠依據所分析的信號,選擇合理的窗函數和帶寬!
平均周期圖法和其它方法求出的結果,參數條件取得一樣的話。能夠得到全然同樣的結果。