matlab——功率谱密度计算 2020年四月

功率谱

功率谱是功率谱密度函数的简称，它定义为单位频带内的信号功率。它表示了信号功率随着频率的变化情况，即信号功率在频域的分布状况。功率谱表示了信号功率随着频率的变化关系 。

常用于功率信号（区别于能量信号）的表述与分析，其曲线（即功率谱曲线）一般横坐标为频率，纵坐标为功率。周期性连续信号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