btestimate_hz.m

在现代谱估计中

clear all
clc
% 产生信号
sample=32;    % 采样点
p=10;         % AR模型的阶数
m=sqrt(-1);
f1=0.05;f2=0.30;f3=0.42;
a(1)=-0.850848;delta=0.101043;
ur(1)=sqrt(delta/2)*randn(1);
ui(1)=sqrt(delta/2)*randn(1);
u(1)=ur(1)+j*ui(1);z(1)=u(1);
x(1)=2*cos(2*pi*f1)+2*cos(2*pi*f2)+2*cos(2*pi*f3)+z(1);
for n=2:sample
    ur(n)=sqrt(delta/2)*randn;
    ui(n)=sqrt(delta/2)*randn;
    u(n)=ur(n)+j*ui(n);z(n)=-a(1)*z(n-1)+u(n);
    x(n)=2*cos(2*pi*f1*n)+2*cos(2*pi*f2*n)+2*cos(2*pi*f3*n)+z(n);
end 
%定义f范围
fmin=-0.5;fstep=0.001;fmax=0.5;
f=fmin:fstep:fmax;
nf=(fmax-fmin)/fstep+1;   %模拟的点数
t=sqrt(-1);

%计算自相关函数的值
rxx=zeros(1,32); rxx0=zeros(1,32);   %定义自相关函数
for k=0:sample-1       %k的取值范围
    k=k+1;
    for n=1:sample-k     %n的取值范围
        rxx(k)=rxx(k)+x(n+k-1)*conj(x(n));
    end
    rxx(k)=(1/sample)*rxx(k);
    rxx0(k)=conj(rxx(k));
end

%功率谱
pxx=zeros(1,nf);w=zeros(1,11);
for k=1:11
    w(k)=1-2*(k-1)/(21-1);
end
for j=1:nf
    for k=(-11):(-2)
        pxx(j)=pxx(j)+rxx0(-k)*w(-k)*exp(-t*2*pi*f(j)*(k+1));
    end
    for k=1:11
        pxx(j)=pxx(j)+rxx(k)*w(k)*exp(-t*2*pi*f(j)*(k-1));
    end
    pxx(j)=10*log10(pxx(j));
end
plot(f,pxx)
title('BT法f2=0.3')