www.gusucode.com > GPS仿真Matlab编程源码程序 > GPS仿真Matlab编程源码程序/LS_QR_16.m
clc
clear all
close all
fs=20e6; % sampling freq
fc=0; %
N=1200;
NN=700;
load s1.mat % 训练序列
load s2.mat % interfere signal,AM
load s4.mat
% load s6.mat
sss=[s4(1:NN) s1(1:N) s4(1:1e4)]; % 发射数据 s1训练序列，s4有用信息
s3=hilbert(s2); % change to complex
%s3、s6干扰信号
sir=-10;
pii=10^(-sir/20);
s=[sss;pii*s3(1:N+NN+1e4)]; %s3干扰信号
% a1=60; %有用信号方向
% a2=20; % 干扰信号方向 20
% fai1=40;%有用信号俯仰角
% fai2=10;% 干扰信号俯仰角 10
a1=120; %有用信号方向 120
a2=90; % 干扰信号方向 90
fai1=20;%有用信号俯仰角 20
fai2=60;% 干扰信号俯仰角 60
NTrial=1 %试验次数
for nTrial = 1:NTrial
M=4;
a11=exp(-j*pi*[0:M-1]'*(cos(a1*pi/180)*cos(fai1*pi/180)));
a22=exp(-j*pi*[0:M-1]'*(sin(fai1*pi/180)));
a=kron(a11,a22); %有用信号导向矢量
b11=exp(-j*pi*[0:M-1]'*(cos(a2*pi/180)*cos(fai2*pi/180)));
b22=exp(-j*pi*[0:M-1]'*(sin(fai2*pi/180)));
b=kron(b11,b22); %干扰信号导向矢量
% c=exp(-j*pi*[0:M-1]'*sin((a1+5)*pi/180));
A=[a,b];
snr=10; % snr: in dB
% randn('state',000)
x=A*s+10^(-snr/20)*randn(4*M,N+NN+1e4);
x=floor(2^10*x);
K=400; %一个周期的采样点
K1=128;
K2=250; % 参与LS运算的快拍数
gold_rate=6.4*1e6;
w0=eye(4*M,1);
t1=0; % 采样开始时刻, 20 samples shift
T1=5e-3; % 采样持续时间
m=0;
%--------xiangguan LS ------------
for n=1:6
y=w0'*x(:,1+(n-1)*K:K+(n-1)*K);
[yy,ptinit1]=xiangguan(y,fs,gold_rate,1,K,8,K1);
% ptinit1=ptinit1+3;
aa(n)=max(abs(yy));
if aa(n)>4
delay1(n)=ptinit1/fs;
sr=signalgen(1,t1,fs,T1,delay1(n),fc,gold_rate,K1);
break
end
end % end of xiangguan
r=sr(1:K2);
ADData=x(:,1+(n-1)*K:K2+(n-1)*K);
SprMseqLocal=r;
NLS=K2;
NULA=4*M;
RxFloor = ADData(:,1:NLS)*ADData(:,1:NLS)';
pFloor = ADData(:,1:NLS)*SprMseqLocal(1:NLS)';
%%%%%%%%%%%%%%%%%%均方域矩阵直接求逆%%%%%%%%%%
p=x(:,1+(n-1)*K:K2+(n-1)*K)*r';
R_1=floor((x(:,1+(n-1)*K:K2+(n-1)*K))*(x(:,1+(n-1)*K:K2+(n-1)*K))');
w_shift=0;
w=((2^w_shift*inv(R_1))*p);
% tic
w_cholesky=zxg_16(R_1,p,NLS);
% toc
%% 按行逐列地清零，变成三角阵，与快拍进来的次序等同（行Nrow>列Nclo）
ADData1 = ADData.';
Rx =ADData1(1:NLS, 1:NULA);
sHous = floor(Rx'*SprMseqLocal(1:NLS)'/2^4); % 需要取共轭
%%%------------------------------定点---------------------------------
[Nrow, Nclo] = size(Rx);
Q = eye(Nrow);
T = zeros(2);
B = Rx(1, :);
%%%%%%%%%%%%%%%%%%%%%%%%%----- Householder变换-----%%%%%%%%%%%%%%%%%%%%%%%%%
H = eye(Nrow);
NA = 16;
B = Rx;
% B = floor(2^NA * B);
B = floor(B);
for k=1:Nclo
e1 = [1; zeros(Nrow-k,1)];
alpha1 = norm(B(k:Nrow, k))*B(k, k)/norm(B(k:Nrow, k)'*e1); % 使得alpha1*B(:,1)'*e1为实数
alpha1 = floor(alpha1);
b1(k) = log2(max(max(abs(alpha1))));
u1_tmp = B(k:Nrow, k) - alpha1*e1;
u1_tmp = floor(u1_tmp);
b2(k) = log2(max(max(abs(u1_tmp))));
u1 = u1_tmp;
b3(k) = log2(max(max(abs(u1))));
H1 = zeros(Nrow);
H1(1:k-1, 1:k-1) = 2^NA * eye(k-1);
% H1(k:Nrow, k:Nrow) = 2^NA * (eye(Nrow-k+1) - 2*u1*u1'/(norm(u1_tmp))^2);
H1(k:Nrow, k:Nrow) = 2^NA * eye(Nrow-k+1) - 2^NA *2*u1*u1'/(norm(u1_tmp))^2;
H1(k:Nrow, k:Nrow) = floor(H1(k:Nrow, k:Nrow));
H = floor(H1*H/2^NA);
B = floor(H1*B/2^NA);
end
R=B;
Q=H;
maxv=max(max(Q*R-Rx));
%%%---------------------------共同部分-----------------------------------------
%%% 前向回代求y
C = R';
C = C(1:Nclo, 1:Nclo); %NA为C扩大的倍数
NB=7;
yHous = zeros(Nclo, 1);
yHous(1) = 2^NB * sHous(1)/C(1, 1);
yHous(1) = floor(2^NB*yHous(1));
b5(1) = log2(max(max(abs(yHous(1)))));
for k=2:Nclo
yHous(k) = floor(2^NB*2^NB*sHous(k)/C(k,k)-C(k,1:k-1)*yHous(1:k-1)/C(k,k));
% yHous(k) = 2^NB*sHous(k)/C(k,k)-C(k,1:k-1)*(yHous(1:k-1)/2^NB)/C(k,k);
% yHous(k) = floor(2^NB * yHous(k));
b5(k) = log2(max(max(abs(yHous(k)))));
end
%%后向回代求w
WHous = zeros(Nclo,1);
NC = 18;
WHous(Nclo) = yHous(Nclo)/R(Nclo, Nclo);
WHous(Nclo) = floor(2^NC * WHous(Nclo));
for k = Nclo-1:-1:1
WHous(k) = floor(2^NC * (yHous(k)/R(k, k))-R(k, k+1:Nclo)*(WHous(k+1:Nclo)/R(k, k)));
b6(k) = log2(max(max(abs(WHous(k)))));
end
WHous = conj(WHous);
w44=reshape(WHous,4,4);
Q=300;
theta=linspace(0,pi,Q);
fai=linspace(0,pi/2,Q);
A=exp(-j*pi.*[0:M-1]'*(kron(cos(theta),cos(fai))));
B=exp(-j*pi.*[0:M-1]'*(kron(ones(1,Q),sin(fai))));
for i=1:4
Grow(i,:)=conj(w44(i,:))*A;
Gcol(i,:)=w44(:,i)'*B;
end
GROW=sum(Grow);
GCOL=sum(Gcol);
G1=GROW.*GCOL;
G1=reshape(G1,Q,Q);
figure(1), mesh(theta*180/pi,fai*180/pi,20*log10(abs(G1)./max(max(abs(G1)))))
ylabel('仰角 (度)');
xlabel('方位角 (度)');
zlabel('增益 (dB)');
hold on
plot3([a1,a1],[fai1,fai1],[-50,0],'k','linewidth',1.5)
plot3([a2,a2],[fai2,fai2],[-50,0],'linewidth',1.5)
w_QR=WHous/norm(WHous);
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w44_1=reshape(w_cholesky,4,4);
for i=1:4
Grow_1(i,:)=conj(w44_1(i,:))*A;
Gcol_1(i,:)=w44_1(:,i)'*B;
end
GROW_1=sum(Grow_1);
GCOL_1=sum(Gcol_1);
G1_1=GROW_1.*GCOL_1;
G1_1=reshape(G1_1,Q,Q);
figure(2), mesh(theta*180/pi,fai*180/pi,20*log10(abs(G1_1)./max(max(abs(G1_1)))))
ylabel('仰角 (度)');
xlabel('方位角 (度)');
zlabel('增益 (dB)');
hold on
plot3([a1,a1],[fai1,fai1],[-50,0],'k','linewidth',1.5)
plot3([a2,a2],[fai2,fai2],[-50,0],'linewidth',1.5)
w_ls=w_cholesky/norm(w_cholesky);
SINR_LS_QR(nTrial) = 10*log10((norm(w_QR'*a))^2/(((norm(w_QR'*b))^2*10^((-sir)/10)+10^(-snr/10))));
SINR_LS_cholesky(nTrial) =10*log10((norm(w_ls'*a))^2/(((norm(w_ls'*b))^2*10^((-sir)/10)+10^(-snr/10))));
end