zuoye2.m

一个数字信号处理的功率谱估计的例子

clear all;
close all;
Fs=2000;
n=0:1/2000:1;
x=sin(2*pi*100*n)+4*sin(2*pi*500*n)+randn(size(n));
N=4096;


% 用 Welch 平均估计试验数据的功率谱；
xpsd1=pwelch(x,hamming(256),32,N,'twosided'); 
% f=(0:length(px)/2-1)*Fs/length(px);
% plot(f,px(1:length(px)/2));%for MATLAB 6.1
% pmax=max(xpsd1);
% xpsd1=xpsd1/pmax;
% 
xpsd1=10*log10(xpsd1);
fn=-0.5:1/4096:0.5-1/4096;
% subplot(4,2,1);
plot(fn,fftshift(xpsd1));grid on;
% xlabel('welch 法平均功率谱估计');
% 使用 Burg 算法得到功率谱估计；
[xpsd2,f2]=pburg(x,36,N,'twosided');
pmax=max(xpsd2);
xpsd2=xpsd2/pmax;
xpsd2=10*log10(xpsd2+0.000001);
% subplot(4,2,2);
figure
plot(f2/(2*pi),fftshift(xpsd2));grid on;
xlabel('Burg 算法得到功率谱估计');
% % 使用协方差算法得到功率谱估计；
% [xpsd3,f3]=pcov(x,36,N);
% pmax3=max(xpsd3);
% xpsd3=xpsd3/pmax3;
% xpsd3=10*log10(xpsd3+0.000001);subplot(4,2,3);
% plot(f3,fftshift(xpsd3));grid on;
% xlabel('协方差算法得到功率谱估计');
% % 使用自相关矩阵分解的特征向量算法得到功率谱估计；
% [xpsd4,f4,V,E]=peig(x,36,N);
% pmax4=max(xpsd4);
% xpsd4=xpsd4/pmax4;
% xpsd4=10*log10(xpsd4+0.000001);
% for i=1:N
%     xxpsd(i)=xpsd4(N+1-i);
% end
% subplot(4,2,4);
% plot(f4,fftshift(xxpsd));grid on;
% 
% xlabel('自相关矩阵分解的特征向量算法得到功率谱估计');
% % 使用改进的协方差算法得到功率谱估计；
% [xpsd5,f5]=pmcov(x,36,N);
% pmax5=max(xpsd5);
% xpsd5=xpsd5/pmax5;
% xpsd5=10*log10(xpsd5+0.000001);subplot(4,2,5);
% plot(f5,fftshift(xpsd5));grid on;
% xlabel('改进的协方差算法得到功率谱估计');
% % 使用MTM算法得到功率谱估计；
% [xpsd6,f6]=pmtm(x,4,N);
% pmax6=max(xpsd6);
% xpsd6=xpsd6/pmax6;
% xpsd6=10*log10(xpsd6+0.000001);subplot(4,2,6);
% plot(f6,fftshift(xpsd6));grid on;
% xlabel('MTM算法得到功率谱估计');
% % 使用自相关矩阵分解的 MUSIC 算法得到功率谱估计；
% [xpsd7,f7]=pmusic(x',36,N);
% pmax7=max(xpsd7);
% xpsd7=xpsd7/pmax7;
% xpsd7=10*log10(xpsd7+0.000001);
% for i=1:N
%     xxpsd(i)=xpsd7(N+1-i);
% end
% subplot(4,2,7);
% plot(f7,fftshift(xxpsd));grid on;
% xlabel('自相关矩阵分解的 MUSIC 算法得到功率谱估计');
% % 使用自相关法得到功率谱估计；
% [xpsd8,f8]=pyulear(x,36,N);
% pmax8=max(xpsd8);
% xpsd8=xpsd8/pmax8;
% xpsd8=10*log10(xpsd8+0.000001);
% subplot(4,2,8);
% plot(f8,fftshift(xpsd8));grid on;
% xlabel('自相关法得到功率谱估计');