为什么周期图法和自相关函数法求的均方值相差很大

905lili

周期图法求的均方值sigma及图像

1. >> clear;
   
2. Fs=1000; %采样频率
   
3. n=0:1/Fs:1;
   
4. xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
   
5. window=boxcar(length(xn));
   
6. nfft=1024;
   
7. [Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
   
8. plot(f,Pxx);
   
9. xlabel('频率/Hz');
   
10. ylabel('功率谱密度');
    
11. sigma=sqrt(sum(Pxx)*Fs/nfft)

复制代码

1. 
   
2. sigma =
   
3. 
   
4.     2.4096

复制代码

自相关函数法求的sigma及图像

1. clear;
   
2. Fs=1000; %采样频率
   
3. n=0:1/Fs:1;
   
4. 
   
5. %产生含有噪声的序列
   
6. xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
   
7. 
   
8. nfft=1024;
   
9. cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
   
10. CXk=fft(cxn,nfft);
    
11. Pxx=abs(CXk);
    
12. index=0:round(nfft/2-1);
    
13. k=index*Fs/nfft;
    
14. plot_Pxx=Pxx(index+1);
    
15. plot(k,plot_Pxx);
    
16. xlabel('频率/Hz');
    
17. ylabel('功率谱密度');
    
18. sigma=sqrt(sum(Pxx)*Fs/nfft)

复制代码

1. sigma =
   
2. 
   
3.   181.2802

复制代码

。
为什么两种 方法求得的 sigma 相差很大，各自的纵坐标为什么也不相同，还有就是有没有必要再用10*log10(Pxx)，把纵坐标化成分贝的单位。求高手指点

905lili

哦，所求sigma应该为标准差

ChaChing

Ref:
5.[原创]关于FFT的频谱对应关系 http://forum.vibunion.com/thread-23614-1-1.html
  Matlab编程实现FFT实践及频谱分析 http://forum.vibunion.com/thread-51332-1-1.html
  FFT方法和FFT的应用注意点 http://forum.vibunion.com/thread-77716-1-1.html
  [原创]关于FFT的结果为什么要除以N http://forum.vibunion.com/thread-23665-1-1.html
  频谱与功率谱有什么区别? http://forum.vibunion.com/thread-8098-1-1.html
  信号经过FFT后，该怎样计算幅值和相位 http://forum.vibunion.com/thread-53683-1-1.html
  matlab中关于FFT的使用(理解频率分辨率、补零问题) http://forum.vibunion.com/thread-98535-1-1.html
  如何用matlab求PSD曲线 http://forum.vibunion.com/thread-72741-1-1.html
From http://forum.vibunion.com/home-space-uid-63979-do-blog-id-18250.html

905lili

回复 3 # ChaChing 的帖子

所给的资料 里面有我想要的答案?先谢谢了

ChaChing

回复 4 # 905lili 的帖子

个人水平/时间有限, 仅瞄了下, 并没有仔细试试看!:@L
感觉LZ看过这些帖应该就知道(不是有特别涂色!?):@)

905lili

回复 5 # ChaChing 的帖子

周期图法(periodogram）是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得x(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。难道你是说上面的自相关函数法没有除N？我理解 好像就是这样的

ChaChing

回复 6 # 905lili 的帖子

个人也并非这方面专家, 但LZ的理解方向对了
不过请LZ再仔细研读下并试试看, 尤其红色那帖许多信号处理专家的讨论
等LZ分享成果:@)

905lili

回复 7 # ChaChing 的帖子

谢谢了，我试过了，就是要除N，但结果还不完全相同，但问题就出在那了

ChaChing

回复 8 # 905lili 的帖子

请LZ再仔细研读下
试下*2/N

905lili

回复 9 # ChaChing 的帖子

也试过了，结果还是相差有点大，

1. clear;
   
2. Fs=1000; %采样频率
   
3. n=0:1/Fs:1;
   
4. 
   
5. %产生含有噪声的序列
   
6. xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
   
7. 
   
8. nfft=1024;
   
9. cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
   
10. CXk=fft(cxn,nfft);
    
11. Pxx=abs(CXk)*2/nfft;
    
12. index=0:round(nfft/2-1);
    
13. k=index*Fs/nfft;
    
14. plot_Pxx=Pxx(index+1);
    
15. plot(k,plot_Pxx);
    
16. xlabel('频率/Hz');
    
17. ylabel('功率谱密度');
    
18. sigma=sqrt(sum(Pxx)*Fs/nfft)

复制代码

1. 
   
2. sigma =
   
3. 
   
4.     8.0444

复制代码

我想两种 方法的误差是会有，但结果不会相差这么大，不知道怎么解释

905lili

有谁用功率谱密度做过疲劳寿命分析呢，

ab77977

回复 11 # 905lili 的帖子

MATLAB中xcorr求解的相关函数存在一定的错误，其原因在什么地方

905lili

回复 12 # ab77977 的帖子

不知道在什么地方:@(

zengwj0423

分享的资料太有价值了