线性预测(liner_prediction.m)

【核心代码】

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 8、线性预测算法
%%%%%%@@@@@@@@@
%%%%%%空间谱估计中的线性预测算法，包括前向预测算法，后向预测算法，双向预测算法，多阶线性预测算法◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎
clc;
clear;
close all;
tic
M=16;%阵元数目
N=3;%信源数目
B=8;%形成的波束个数
snap=1000;%快拍数目
C=3e8;
f0=10e6;
lamda=C/f0;
d=0.5*lamda;
% k=d/lamda;
theta0=20;
theta1=30;
theta2=-10;
fs=1000;
ts=1/fs;
t=(0:snap-1)*ts;
a=[0:M-1]';%阵列矢量
u0=5;
u1=10;
u2=20;


s0=exp(j*2*pi*0.5*u0*t.^2);
s1=exp(j*2*pi*0.5*u0*t.^2);
s2=exp(j*2*pi*0.5*u0*t.^2);

%阵列流行矢量

a_theta0=exp(j*2*pi*d/lamda*a*sin(theta0/180*pi));
a_theta1=exp(j*2*pi*d/lamda*a*sin(theta1/180*pi));
a_theta2=exp(j*2*pi*d/lamda*a*sin(theta2/180*pi));

A=[a_theta0 a_theta1 a_theta2];
S=[s0;s1;s2];
X0=A*S;
SNR=15;
randn('state',0)
real_noise=randn(size(X0));
randn('state',3);
imag_noise=randn(size(X0));
noise0=(real_noise j*imag_noise)/2^0.5;
noise=10^(-SNR/20)*noise0;
X=X0 noise;%%%%%%%%%%------------=========完整的基带信号
Rx=X*X'/length(t);
%%%@@@@@@@@@@@@@@@@下面为一阶线性预测算法
%%%%%%%%%%前向预测算法
Xf=flipud(X(1:M-1,:));
Xft=transpose(Xf);%转置
Xfh=Xf';%共轭转置
Rf=Xf*Xfh/snap;
rf=Xf*X(M,:)'/snap;
Wflp=conj((inv(Rf)*rf));%复共轭
%%%%%%%%%%后向预测算法
Xb=X(2:M,:);
Xbt=transpose(Xb);%转置
Xbh=Xb';%共轭转置
Rb=Xb*Xbh/snap;
rb=Xb*X(1,:)'/snap;
Wblp=inv(Rb)*rb;
%%%%%%%%%%双向预测算法
Xfbt=[Xft;conj(Xbt)];Xfb=transpose(Xfbt);Xfbh=Xfb';
Rfb=Xfb*Xfbh/snap;rfb=Xfb*[X(M,:)';transpose(X(1,:))]/snap;
Wfblp=conj(inv(Rfb)*rfb);

%%%@@@@@@@@@@@@@@@@下面为多阶线性预测算法
p=10;P=M-p;%线性预测的阶数
J=fliplr(eye(P));%P阶交换矩阵
RF0=zeros(P,P);%前向预测的协方差阵
Rf0=zeros(P,1);
af0=zeros(P,1);
RB=zeros(P,P);%后向预测的协方差阵
Rb=zeros(P,1);


theta=-90:0.01:90;
for ii=1:length(theta)
    a_theta=exp(j*2*pi*d/lamda*a*sin(theta(ii)/180*pi));
    Pflp(ii)=10*log10(1/(abs(a_theta'*[1;-Wflp]))^2);
    Pblp(ii)=10*log10(1/(abs(a_theta'*[1;-Wblp]))^2);   
    Pfblp(ii)=10*log10(1/(abs(a_theta'*[1;-Wfblp]))^2); 
    for jj=1:p
        RF=J*Rx(jj:P jj-1,jj:P jj-1)*J;
        Rf=J*Rx(jj:P jj-1,P jj);
        WFLP=conj((inv(RF)*Rf));
        PFLP(jj,ii)=10*log10(1/(abs(a_theta(1:P 1)'*[1;-WFLP]))^2);
        RB=Rx(jj 1:P jj,jj 1:P jj);
        Rb=Rx(jj 1:P jj,jj);
        WBLP=(inv(RB)*Rb);
        PBLP(jj,ii)=10*log10(1/(abs(a_theta(1:P 1)'*[1;-WBLP]))^2); 
        RFB=RF conj(RB);
        Rfb=Rf conj(Rb);
        WFBLP=conj((inv(RFB)*Rfb));   
        PFBLP(jj,ii)=10*log10(1/(abs(a_theta(1:P 1)'*[1;-WFBLP]))^2);
    end          
end

figure(1);
plot(theta,Pflp);grid on;xlabel('角度');ylabel('峰值');title('一阶前向预测');
figure(2);
plot(theta,Pblp);grid on;xlabel('角度');ylabel('峰值');title('一阶后向预测');
figure(3);
plot(theta,Pfblp);grid on;xlabel('角度');ylabel('峰值');title('一阶双向预测');
figure(4);
plot(theta,sum(PFLP,1)/p);grid on;xlabel('角度');ylabel('峰值');title('多阶前向预测');
figure(5);
plot(theta,sum(PBLP,1)/p);grid on;xlabel('角度');ylabel('峰值');title('多阶后向预测');
figure(6);
plot(theta,sum(PFBLP,1)/p);grid on;xlabel('角度');ylabel('峰值');title('多阶双向预测');
toc