clear all;close all;clc;
% 噪声方差1,均值0,实高斯白噪声
sigma=1;
% 三个信号的信噪比
snr=[30 30 27];
% 归一化频率
f=[0.1 0.25 0.27];
% 各个信号的随机相位
for ii=1:3
phi(ii)=unifrnd(0,2*pi);
end
% 根据信噪比计算三个信号的幅度
a=sqrt(10.^(snr/10)*2*sigma^2);
nfft=1024;
N=256;
n=0:N-1;
% AR模型的阶数 p,参数p可能带来的问题:p过大,可能出现假峰,p过小,可能不能准确估计峰值
p=16;
% 生成各信号s(i)
s=zeros(3,N);
for ii=1:3
s(ii,:)=a(ii)*cos(2*pi*f(ii)*n+phi(ii));
end
% 生成噪声w
w=randn(1,N);
un=sum(s)+w;
% 计算观测数据的自相关函数 cun
cun=xcorr(un,'unbiased');
% 取自相关函数的cun(0)~ cun(p)
cun_ar=cun(N:N+p);
% 计算AR模型的系数ak(1)~ak(p)
Rp=toeplitz(cun_ar(1:p));
rp=cun_ar(2:p+1)';
ak=-inv(Rp)*rp;
% 估计AR模型的输入白噪声方差
var=abs(cun_ar*[1,ak']');
w=linspace(-pi,pi,nfft);
e_w=zeros(1,nfft);
for m=1:p
e_w=e_w+ak(m)*exp(-j*w*m);
end
S=var./(abs(1+e_w)).^2;
plot_S=10*log10(S/max(S));
w_norm=w/2/pi;
plot(w_norm,plot_S);
xlabel('w/2pi');
ylabel('归一化功率谱/dB');