wenti.m

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

Fs=1000;
NFFT=64;
n=0:1/Fs:1;
load weicha1.txt
x=weicha1;
% x=sin(2*pi*100*n)+sin(2*pi*300*n)+sin(2*pi*450*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);
w=hanning(NFFT)';
Pxx=(abs(fft(w.*x(1:NFFT))).^2+... 
    abs(fft(w.*x(NFFT*1/2+1:NFFT*3/2))).^2+... 
    abs(fft(w.*x(NFFT*2/2+1:NFFT*4/2))).^2+... 
    abs(fft(w.*x(NFFT*3/2+1:NFFT*5/2))).^2+... 
    abs(fft(w.*x(NFFT*4/2+1:NFFT*6/2))).^2+... 
    abs(fft(w.*x(NFFT*5/2+1:NFFT*7/2))).^2)/(norm(w)^2*6);
f=(0:(NFFT-1))./NFFT*Fs;%hongzuobiao

% [P,f]=pburg(x,18,NFFT,Fs);
% % [P2,f]=pmcov(x,20,[],Fs);
PXX=10*log10(Pxx);
% Pxx2=10*log10(P2);
figure
plot(f,PXX)

xlabel('频率(Hz))');
ylabel('相对功率谱密度(dB/Hz)');
grid;

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