MATLAB功率谱函数psd与pwelch的使用
目录
一、遇到的问题
在网上看到这样一个问题，然后搬过来记一下。
大家好，我用matlab的psd函数和pwelch函数求功率谱密度，但是得到的结果为什么下差大约30db，这个是怎么回事？
摘自：https://www.ilovematlab.cn/thread-339674-1-1.html
clear;
Fs=1000;
n=0:1/Fs:1;
y=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
Nfft=1000;
window=hanning(100);
noverlap=20;
dflag='none';
[Pxx,f]=psd(y,Nfft,Fs,window,noverlap);
figure
plot(f,10*log10(Pxx));
xlabel('frequency/Hz');ylabel('PSD/dB');
title('PSD-Welch method')
hold on;
[Pxx1,f1]=pwelch(y,window,noverlap,Nfft,Fs);
plot(f1,10*log10(Pxx1),'r');
针对这个问题的答案是：
| 程序中用PSD的用法是较早版本中的函数用法，在现在的新版本中已不这样用了，是结合了spectrum函数一起使用。在较早版本中PSD函数得到的是功率谱，不是功率谱密度，所以与pwelch的结果（计算功率谱密度）会有30db的差别。 |
怎么解决这个问题呢？
| [Pxx1,f1]=pwelch(y,window,noverlap,Nfft,Fs); %加入下面一段 Pxx1=Pxx1*Fs/2.0 Pxx1(1)=Pxx1(1)*2; Pxx1(length(Pxx1))=Pxx1(length(Pxx1))*2; 结果完全一样 |
二、功率谱函数psd与pwelch的用法与区别
psd与pwelch均为求功率谱的函数，但是psd求得的功率谱不可直接使用还需要进行一定的变换。
现给出二者的等价关系：
[px,fx]=psd(x,nfft,fs,Windows,noverlap);
px=px/fs*2;%---自谱
px([1,end])=px([1,end])/2; %---直流分量
与之等价的pwelch函数表达为：
[px,fx]=pwelch(x,Windows,noverlap,nfft,fs);
其中，
x --------- 时程信号；
nfft--------傅里叶变换（fft）点数；
Windows--窗函数，通常设为hanning(nfft)
nooverlap-重叠的点数，通常为nfft/2；
px---------求得的信号x的功率谱密度；
fx----------对应px的频率序列。
需要说明的是，若中间有个变量不知如何设定（如noverlap），却仍然想设定该变量后面的值，
可采用中括号[]代替，[]则表示该变量并未设定具体值，而是采用函数默认值来取。
下面以一具体算例给大家做演示：
clear
x=randn(10000,1);%生成随机信号
fs=10; %---------采样频率
nfft=1024;%-------fft点数
%---psd求谱
[px,fx]=psd(x,nfft,fs,hanning(nfft),nfft/2);
px=px/fs*2;
px([1,end])=px([1,end])/2;
%---pwelch求谱
[px2,fx2]=pwelch(x,hanning(nfft),nfft/2,nfft,fs);
%---绘图
figure
semilogy(fx,px);hold on
semilogy(fx2,px2,'r')
legend('psd','pwelch')
title('psd与pwelch方法求谱结果比较')运行以上程序得到下图，由图可见psd曲线与pwelch曲线完全重合。
