clear all;clc;
N=256;
SNR1=30;SNR2=30;SNR3=27;%设定信噪比；
A1=getAk(SNR1);A2=getAk(SNR2);A3=getAk(SNR3);%求得信号的幅度；
f1=0.1;f2=0.25;f3=0.27;%设定频率值
s1=zeros(1,N);s2=zeros(1,N);s3=zeros(1,N);%初始化sk
s1=getsk(A1,f1,N);%求sk
s2=getsk(A2,f2,N);
s3=getsk(A3,f3,N);
v=randn(1,N);% 构建高斯白噪声；
u=s1+s2+s3+v;%构建了时域上的信号
N2=2*N;
u2=zeros(1,N2);
u2(1:N)=u;%对信号序列进行补零
UW=fft(u2);%求DFT变换
SW1=zeros(1,N2);
for w=1:N2
SW1(w)=((abs(UW(w)))^2)/N;
end
r0=zeros(1,N2);
r0=ifft(SW1);
r=zeros(1,N2);
r(1:N)=r0(N+1:N2);
r(N+1:N2)=r0(1:N);
p=16;%设定阶数
R=zeros(p,p);
for i=1:p
for j=1:p
R(i,j)=r(N+1+i-j);%得到了p阶的矩阵R 注意：这里面r(0)对应这我们的r(N+1),因为需要的r是从-N到N-1的，但是matlab中是不能存在r(-N)的，就将r(-N)对应这r(1),需要的r(0)对应着r(N+1).
end
end
rp=zeros(p,1);
for i=1:p
rp(i)=r(N+1+i);
end
ep=zeros(p,1);
ep=-(R^-1)*rp;
tempr=zeros(1,p);%为了求δ2而设定的临时行向量
for i=1:p
tempr(i)=r(N-i+1);
end
b2=r(N+1)+tempr*ep;%δ2
M=1024;
a=zeros(1,M);
a(1)=1;
for i=2:(p+1)
a(i)=ep(i-1);
end
SW=fft(a);
temp=zeros(1,M/2);
temp=SW(1:M/2);
SW(1:M/2)=SW((M/2+1):M);
SW((M/2+1):M)=temp;%将0-2*pi的序列折换成-pi到pi的序列。
SAR=zeros(1,M);
for w=1:M
SAR(w)=log10(abs(b2/((abs(SW(w)))^2)));
end
d=1/(M-1);
x=zeros(1,M);
for i=1:M
x(i)=-0.5+(i-1)*d;
end
plot(x,SAR);
评论0