各位是如何用matlab求PSD曲线的,急用

seabird08

各位是如何用matlab求PSD曲线的,急用,教教小弟吧

ch_j1985

能否介绍一下PSD曲线？

seabird08

就是自功率谱密度函数曲线

yannan1003

你可以查看一下matlab中自带的帮助文件，help psd

wwbeyondww

用fft可以吧？

vincentsuen

% 功率谱分析
clear;
clc;
fs=50;
ns=512;
t=1:ns;
t=t/fs;
z=sin(4*pi*t);
[Pxx,w]=periodogram(z,hamming(512),ns,fs,'onesided');
figure(2);
plot(w,Pxx);
title('功率谱函数');
xlabel('频率（Hz）');
ylabel('S(w)(m2/s3)');
grid;

[ 本帖最后由 vincentsuen 于 2008-11-28 21:49 编辑 ]

ChaChing

看到楼上的, 想到顺便考考以下的结果
clc; fs=50;
for ns=[256,512,1024,2048]
  t=1:ns; t=t/fs; z=sin(4*pi*t);
  [Pxx,w]=periodogram(z,hamming(ns),ns,fs,'onesided');
  figure; plot(w,Pxx); grid; title('spectrum'); xlabel('Hz'); ylabel('S(w)');
end

想想为何峰值都不同!?

vincentsuen

原帖由 ChaChing 于 2008-11-28 22:51 发表

                               
登录/注册后可看大图

看到楼上的, 想到顺便考考以下的结果
clc; fs=50;
for ns=[256,512,1024,2048]
  t=1:ns; t=t/fs; z=sin(4*pi*t);
  [Pxx,w]=periodogram(z,hamming(ns),ns,fs,'onesided');
  figure; plot(w,Pxx); grid; titl ...

峰值不同的原因，在于我们用有限长的样本来估计无限长的随机过程，当然会引起误差，会把数据包括范围之外的能量泄露掉，等于加了一个窗（0到ns之外的数据都被乘了0，之内的数据的都被乘了1）。从sin(4*pi*t)这个例子中可以看到，信号应该是只有2一个频率的，就是说只有在横坐标等于2的时候才有值，但是为什么画出来的psd图却是在2的周围一个很小带宽内都有值呢？这就是谱的泄露。所以当信号样本取得越长，2周围的这个带宽越小。同时，由于psd图对频率积分应该等于一个确定的值：均方值，就是说psd图所围的面积应该相等。所以带宽越小，峰值就越大。不信你可以做个试验，不管你的ns取多少，Pxx下所围的面积都是差不多的。解说完毕！

[ 本帖最后由 vincentsuen 于 2008-11-28 23:16 编辑 ]

seabird08

高手啊,学习下,继续讨论.
对于地震波记录来说,其功率谱如果用其中自带函数,好像不是太合适

chenjc18

下面是一些不同方法计算同一信号的matlab 程序！
直接法：
直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
Matlab代码示例：
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)的功率谱估计。
Matlab代码示例：
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)分段求周期图再平均。
Matlab代码示例：
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)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。
Matlab代码示例：
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);

chenjc18

以上是基于傅里叶变换（fft）的psd估计方法
根据上述方法计算出对应频率或者周期下的PSD值 就能画图了 你说的做PSD曲线（通常是log-log 曲线图 具体的下点论文看看 ）

wangyourwang

905lili

回复 8 # vincentsuen 的帖子

问一下，我要求均方值，有没有必要用10*log10(Pxx)把纵坐标化成分贝的单位，我用周期图法和自相关函数法求的得均方值为什么相差很大，求解释

ChaChing

回复 13 # 905lili 的帖子

个人以为求均方值与如何绘图显示应该无关吧! 不是吗!?

905lili

回复 14 # ChaChing 的帖子

是无关，只是看有些贴子化成分贝了，