功率譜
功率譜是功率譜密度函數的簡稱，它定義為單位頻帶內的信號功率。它表示了信號功率隨着頻率的變化情況，即信號功率在頻域的分布狀況。功率譜表示了信號功率隨着頻率的變化關系 。
常用於功率信號（區別於能量信號）的表述與分析，其曲線（即功率譜曲線）一般橫坐標為頻率，縱坐標為功率。周期性連續信號x(t)的頻譜可表示為離散的非周期序列XnXnXn,它的幅度頻譜的平方│Xn│2│Xn│2│Xn│2所排成的序列，就被稱之為該周期信號的“功率譜”。
Matlab 使用
fft做出來是頻譜，psd做出來是功率譜；功率譜丟失了頻譜的相位信息；頻譜不同的信號其功率譜是可能相同的；功率譜是幅度取模后平方，結果是個實數。matlab中自功率譜密度直接用psd函數就可以求，按照matlab的說法，psd能實現Welch法估計，即相當於用改進的平均周期圖法來求取隨機信號的功率譜密度估計。psd求出的結果應該更光滑吧。
psd簡介
PSD(power spectrum analysis)功率譜分析，PSD在給定頻帶上的積分計算信號在該頻帶上的平均功率。與均值-平方譜相反，這個光譜中的峰值並沒有反映出給定頻率的能量。
1 直接法
直接法又稱周期圖法，它是把隨機序列x(n)x(n)x(n)的NNN個觀測數據視為一能量有限的序列，直接計算x(n)x(n)x(n)的離散傅立葉變換，得X(k)X(k)X(k)，然后再取其幅值的平方，並除以NNN，作為序列x(n)x(n)x(n)真實功率譜的估計。
periodogram函數是用來計算功率譜密度的
[Pxx,f]=periodogram(x,window,nfft,Fs);
X：所求功率譜密度的信號；
window：所使用的窗口，默認是問boxcar，其長度必須答與x的長度一致；
nfft：采樣點數；
clear;
Fs=1000; %采樣頻率
n=0:1/Fs:1;
%產生含有噪聲的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
window=boxcar(length(xn)); %矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
subplot(1,2,1);
plot(xn);
subplot(1,2,2);
plot(f,10*log10(Pxx));
2、間接法
間接法先由序列x(n)x(n)x(n)估計出自相關函數R(n)R(n)R(n)，然后對R(n)R(n)R(n)進行傅立葉變換，便得到x(n)x(n)x(n)的功率譜估計。
clear;
Fs=1000; %采樣頻率
n=0:1/Fs:1;
%產生含有噪聲的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
cxn=xcorr(xn,'unbiased'); %計算序列的自相關函數
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
subplot(1,2,1);
plot(xn);
subplot(1,2,2);
plot(k,plot_Pxx);
3、Bartlett平均周期圖的方法是將NNN點的有限長序列x(n)x(n)x(n)分段求周期圖再平均。
clear;
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(length(n)); %矩形窗
noverlap=0; %數據無重疊
p=0.9; %置信概率
[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
subplot(1,2,1);
plot(k,plot_Pxx);
subplot(1,2,2);
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);
4、Welch法對Bartlett法進行了兩方面的修正，一是選擇適當的窗函數w(n)w(n)w(n)，並再周期圖計算前直接加進去，加窗的優點是無論什么樣的窗函數均可使譜估計非負。二是在分段時，可使各段之間有重疊，這樣會使方差減小。
clear;
Fs=1000;
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %數據無重疊
range='half'; %頻率間隔為[0 Fs/2]，只計算一半的頻率
[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f1]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f2]=pwelch(xn,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
subplot(1,3,1);
plot(f,plot_Pxx);
subplot(1,3,2);
plot(f1,plot_Pxx1);
subplot(1,3,3);
plot(f2,plot_Pxx2);
————————————————
以上內容幾乎為復制而來，寫於此處用於記錄。
版權聲明：本文為CSDN博主「MissXy_」的原創文章，遵循 CC 4.0 BY-SA 版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/MissXy_/java/article/details/94972476