【MATLAB】MATLAB实现功率谱密度分析psd及详细解说

嵌入式职场

已于 2022-07-03 17:44:30 修改

阅读量4.8k

点赞数 2

文章标签： matlab 开发语言线性代数

于 2022-02-17 21:07:24 首次发布

本文链接：https://blog.csdn.net/fanjufei123456/article/details/122991947

版权

这篇博客探讨了不同功率谱密度估计方法，包括直接法、间接法和改进的直接法（Bartlett法与Welch法）。通过MATLAB代码展示了如何使用这些方法计算含有噪声信号的功率谱，并对比了它们的幅度差异和适用场景。文章强调了在实际信号处理中选择合适方法的重要性。

摘要由CSDN通过智能技术生成

功率谱密度幅值的具体含义？？

求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！
我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？
功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的matlab 程序！欢迎大家给点建议！
直接法：
直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以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)); 

window=boxcar(length(xn)); %矩形窗 
nfft=1024; 
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法 
plot(f,10*log10(Pxx));

间接法：
间接法先由序列x(n)估计出自相关函数R(n)，然后对R(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)); 
plot(k,plot_Pxx);

改进的直接法：
对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。

1. Bartlett法

Bartlett平均周期图的方法是将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)); 
figure(1) 
plot(k,plot_Pxx); 

pause; 

figure(2) 
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

2. Welch法

Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数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,f]=pwelch(xn,window1,noverlap,nfft,Fs,range); 
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range); 

plot_Pxx=10*log10(Pxx); 
plot_Pxx1=10*log10(Pxx1); 
plot_Pxx2=10*log10(Pxx2); 

figure(1) 
plot(f,plot_Pxx); 

pause; 

figure(2) 
plot(f,plot_Pxx1); 

pause; 

figure(3) 
plot(f,plot_Pxx2);