bartlett.m

给出了用MATLAB语言实现功率谱估计的程序

%两种分段平均周期图法
function y=bartlett()
%求信号x(n)=sin(2*pi*f1*n)+2cos(2*pi*f2*n)+w(n), n=1,···,N-1的功率谱。其中f1=0.05,f2=0.12,w(n)为白噪声
%N=1024
%采用分段平均周期图法估计功率谱
clf;
N=1024;
Ns=256;
n=[0:N-1];
f1=0.05;
f2=0.12;
wn=randn(1,N)
xn=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn;
Pxx1=abs(fft(xn(1:256),Ns).^2)/Ns;
Pxx2=abs(fft(xn(257:512),Ns).^2)/Ns;
Pxx3=abs(fft(xn(513:768),Ns).^2)/Ns;
Pxx4=abs(fft(xn(769:1024),Ns).^2)/Ns;
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);
f=[0:length(Pxx)-1]/length(Pxx);
figure(1);
plot(f,Pxx);
xlabel('Frequency');
ylabel('Power Spectrum(dB)');
title('Averaged Periodogram (no overlap) N=4*256');
grid;
%==========================================
%Estimate PSD by the halfoverlapping signalsection
N=1024;
Ns=256;
n=[0:N-1];
f1=0.05;
f2=0.12;
wn=randn(1,N)
xn=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn;
Pxx1=abs(fft(xn(1:256),Ns).^2)/Ns;
Pxx2=abs(fft(xn(129:384),Ns).^2)/Ns;
Pxx3=abs(fft(xn(257:512),Ns).^2)/Ns;
Pxx4=abs(fft(xn(385:640),Ns).^2)/Ns;
Pxx5=abs(fft(xn(513:768),Ns).^2)/Ns;
Pxx6=abs(fft(xn(641:896),Ns).^2)/Ns;
Pxx7=abs(fft(xn(769:1024),Ns).^2)/Ns;
Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);
f=[0:length(Pxx)-1]/length(Pxx);
figure(2);
plot(f,Pxx);
xlabel('Frequency');
ylabel('Power Spectrum(dB)');
title('Averaged Periodogram (halfoverlap) N=1024');
grid;