pwelchhanshu.m

对不同的功率谱密度分析方法的研究

Fs=200;
NFFT=280;
n=0:1/Fs:1;
load weicha1.txt
x=weicha1;
%  x=sin(2*pi*100*n)+4*sin(2*pi*200*n)+sin(2*pi*400*n)+randn(size(n));
% Cx=xcorr(x,'unbiased');%计算序列的自相关函数
% Cxk=fft(Cx,NFFT);%计算FFT
% pxx=abs(Cxk);%求解PSD
% t=0:round(NFFT/2-1);
% k=t*Fs/NFFT;
% p=10*log(pxx(t+1));
% plot(k,p)
% xlabel('频率(Hz)');
% ylabel('相对功率谱密度(dB/Hz)');
% h=fir1(20,0.3);
% load 100a.txt
% x=X100a;
% x=filter(h,1,r);
window1=boxcar(100);
window2=hamming(100);
window3=blackman(100);
noverlap=20;
[Pxx1,f1]=pwelch(x,window1,noverlap,NFFT,Fs);
[Pxx2,f2]=pwelch(x,window2,noverlap,NFFT,Fs);
[Pxx3,f3]=pwelch(x,window3,noverlap,NFFT,Fs);

% [P,f]=pyulear(x,20,NFFT,Fs);
% % [P2,f]=pmcov(x,20,[],Fs);
% Pxx=10*log10(P);
Pxx1=10*log10(Pxx1);
Pxx2=10*log10(Pxx2);
Pxx3=10*log10(Pxx3);
subplot(3,1,1)
% figure
plot(f1,Pxx1)
title('矩形窗');
subplot(3,1,2)
plot(f2,Pxx2)
title('海明窗');
subplot(3,1,3)
plot(f3,Pxx3)
title('布莱克曼窗');
xlabel('频率(Hz))');
ylabel('相对功率谱密度(dB/Hz)');

% legend('协方差法','改进的协方差法');