burg.m

实现了用burg法进行功率谱估计

%题目给出的参数
N=256;
P=0.12;
a=1-sqrt(0.1);
h=[1,-a];
SNR=[10 50];   
F=[0.1 0.26];          %两个信噪比和归一化频率

%产生高斯白噪声w和通过h(n)之后的信号v
b=sqrt(P);
w=randn(1,20*N);
w=w/std(w);
w=w-mean(w);
w=w*b;

v=filter(1,h,w);
v=v(1000:1000+N-1);

%产生正弦序列
A=sqrt(2*P*10.^(SNR/10));   %根据信噪比确定两个正弦函数的幅度   
x1=A(1)*sin(2*pi*F(1)*[0:N-1]);
x2=A(2)*sin(2*pi*F(2)*[0:N-1]);
%产生题目所要求信号x
x=x1+x2+v;

%求出信号x的真实功率谱
M=2048;H=freqz(1,h,M,1);H=abs(H);H=H.*H;H=P*H;
f=0:1/M:0.5-1/M;
power_y1=pi*A(1)*A(1)/2;
power_y2=pi*A(2)*A(2)/2;
f1=0.1;
f2=0.26;
f1_n=round(f1*M);
f2_n=round(f2*M);

H(f1_n)=H(f1_n)+power_y1;
H(f2_n)=H(f2_n)+power_y2;

PX=10*log10(H);
figure;
plot(f,PX(1:1024),'r');
grid on;
hold on;

%用burg法估计x(n)的功率谱;
[xpsd,F]=pburg(x,9,M,1);
PX1=10*log10(xpsd);
plot(F,PX1,'b');
grid on;