本书与天线理论方面的教科书配套使用,适合散返高等院校通信、无线电等相关专业的教学应用,也可供相关领域的工程技术人员参考。书中的所有天线模型都可在原书配套光盘中找到,便于读者的应用和演练。
你这里有两个程序,第二个程序与第一个实质上是一样的,区别就是信号与导向矢量的写法有点不同,这里我就不注释了。还有,我下面附了一段我自己的写的程序,里面有SIM算法。G-S正交化算法等。是基于圆阵形式的,你的算法是基于线阵的,他们程序上的区别在于导向矢量的不同。我的算法是某项目中的,保证好使。建议学习波束形成技术,注意把程序分块,例如分成,求导向矢量;最优权值;形成波束等等。程序如下:
4单元均匀线阵自适应波束形成图
clear
clc
format long
v=1
M=4
N=1000%%%%%%%快拍数
f0=21*10^3%%%%%%%%%%%信号与干扰的频率
f1=11*10^3
f2=15*10^3
omiga0=2*pi*f0%%%%%%%信号与干扰的角频率
omiga1=2*pi*f1
omiga2=2*pi*f2
sita0=08%信号方向
sita1=04%干扰方向1
sita2=21%干扰方向2
for t=1:N %%%%%%%%%%%%信号
adt(t)=sin(omiga0*t/(N*f0))
a1t(t)=sin(omiga1*t/(N*f1))
a2t(t)=sin(omiga2*t/(N*f2))
end
for i=1:M%%%%%%%%%%%%信号的导向矢量:线阵的银旦形式
ad(i,1)=exp(j*(i-1)*pi*sin(sita0))
a1(i,1)=exp(j*(i-1)*pi*sin(sita1))
a2(i,1)=exp(j*(i-1)*pi*sin(sita2))
end
R=zeros(M,M)
for t=1:N
x=adt(t)*ad+a1t(t)*a1+a2t(t)*a2%阵列对信号的完整响应
R=R+x*x'%信号的协方差矩阵
end
R=R/N%%%%%%%%%协方差矩阵,所有快拍数的平均
miu=1/(ad'*inv(R)*ad)%%%%%%这个貌似是LMS算法的公式,具体我记不太清,这里是求最优权值,根据这个公式求出,然后加权
w=miu*inv(R)*ad
%%%%%%形成波束%%%%%%%%%%%%%%%%%%%
for sita=0:pi/100:pi
for i=1:M
x_(i,1)=exp(j*(i-1)*pi*sin(sita))
end
y(1,v)=w'*x_%%%%%%%对信号进行加权,消除干扰
v=v+1
end
y_max=max(y(:))%%%%%%%%%%%%%%%归一化
y_1=y/y_max
y_db=20*log(y_1)
sita=0:pi/100:pi
plot(sita,y)
Xlabel(‘sitaa’)
Ylabel(‘天线增益db’)
4单元均匀线阵自适应波束形成
目标
clear
clc
format long
v=1
M=4阵元数
N=100
f0=21*10^3
omiga0=2*pi*f0
sita0=06%信号方向
for t=1:N
adt(t)=sin(omiga0*t/(N*f0))
end
for i=1:M
ad(i,1)=exp(j*(i-1)*pi*sin(sita0))
end
R=zeros(4,4)
r=zeros(4,1)
for t=1:N
x=adt(t)*ad
R=R+x*x'
end
R=R/N
miu=1/(ad'*inv(R)*ad)
w=miu*inv(R)*ad
for sita=0:pi/100:pi/2
for i=1:M
a(i,1)=exp(j*(i-1)*pi*sin(sita))
end
y(1,v)=w'斗贺*a
v=v+1
end
sita=0:pi/100:pi/2
plot(sita,y)
xlabel('sita')
ylabel('天线增益’)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%我的程序%%%%%%%%%%%%%%%
function jieshousignal
%期望信空搏派号数:1个
%干扰信号数:4个
%信噪比已知
%干燥比已知
%方位角已知
clc
clear all
close all
%//参数设置===========================================
1=0
2=0
3=0
% for rrr=1:16000
signal_num=1 %signal number
noise_num=5 %interference number
R0=06 %圆的半径
SP=2000 %Sample number
N=8 %阵元数
snr=-10%Signal-to-Noise
sir1=10 %Signal-to-Interference one
sir2=10 %Signal-to-Interference two
sir3=10 %Signal-to-Interf
sir4=10
sir5=10
%//================noise Power-to-signal Power====================
factor_noise_1=10^(-sir1/10)
factor_noise_2=10^(-sir2/10)
factor_noise_3=10^(-sir3/10)
factor_noise_4=10^(-sir4/10)
factor_noise_5=10^(-sir5/10)
factor_noise_targe=10^(-snr/10)
% //======================== ===============
d1=85*pi/180%%干扰1的方位角
d2=100*pi/180%干扰2的方位角
d3=147*pi/180%干扰3的方位角
d4=200*pi/180%干扰4的方位角
d5=250*pi/180%干扰5的方位角
d6=150*pi/180%目标的方位角
e1=15*pi/180%%干扰1的俯仰角
e2=25*pi/180%干扰2的俯仰角
e3=85*pi/180%干扰3的俯仰角
e4=50*pi/180%干扰4的俯仰角
e5=70*pi/180%干扰5的俯仰角
e6=85*pi/180%目标的俯仰角
% //====================目标信号==========================
t=1:1:SP
fc=2e7
Ts=1/(3e10)
S0=5*cos(2*pi*fc*t*Ts)%目标信号
for kk=1:N
phi_n(kk)=2*pi*(kk-1)/N
end
%//==================== *** 纵矢量==========================================
A=[conj(exp(j*2*pi*R0*cos(d6-phi_n)*sin(e6)))conj(exp(j*2*pi*R0*cos(d1-phi_n)*sin(e1)))conj(exp(j*2*pi*R0*cos(d2-phi_n)*sin(e2)))conj(exp(j*2*pi*R0*cos(d3-phi_n)*sin(e3)))conj(exp(j*2*pi*R0*cos(d4-phi_n)*sin(e4)))conj(exp(j*2*pi*R0*cos(d5-phi_n)*sin(e5)))]'
A1=[conj(exp(j*2*pi*R0*cos(d1-phi_n)*sin(e1)))conj(exp(j*2*pi*R0*cos(d2-phi_n)*sin(e2)))conj(exp(j*2*pi*R0*cos(d3-phi_n)*sin(e3)))conj(exp(j*2*pi*R0*cos(d4-phi_n)*sin(e4)))conj(exp(j*2*pi*R0*cos(d5-phi_n)*sin(e5)))]'
% //==========================================================Power of the interference
% // depending on the signal power and SIR
Ps1=0
Ps2=0
Ps3=0
Ps4=0
Ps5=0
S1=zeros(1,SP)
S2=zeros(1,SP)
S3=zeros(1,SP)
S4=zeros(1,SP)
S5=zeros(1,SP)
Ps0=S0*S0'/SP % signal power
Ps1=Ps0*factor_noise_1
Ps2=Ps0*factor_noise_2
Ps3=Ps0*factor_noise_3
Ps4=Ps0*factor_noise_4
Ps5=Ps0*factor_noise_5
% //==========================干扰信号的随机包络=========================
S1=normrnd(0,sqrt(Ps1/2),1,SP)+j*normrnd(0,sqrt(Ps1/2),1,SP)
S2=normrnd(0,sqrt(Ps2/2),1,SP)+j*normrnd(0,sqrt(Ps2/2),1,SP)
S3=normrnd(0,sqrt(Ps3/2),1,SP)+j*normrnd(0,sqrt(Ps3/2),1,SP)
S4=normrnd(0,sqrt(Ps4/2),1,SP)+j*normrnd(0,sqrt(Ps4/2),1,SP)
S5=normrnd(0,sqrt(Ps5/2),1,SP)+j*normrnd(0,sqrt(Ps5/2),1,SP)
%//
S=[S0S1S2S3S4S5]
SS1=[S1S2S3S4S5]
X=A*S%信号加干扰
XX2=A1*SS1%接收到的干扰
Pw_noise=sqrt(Ps0*factor_noise_targe/2)
a1=randn(N,SP)
a2=randn(N,SP)
a1=a1/norm(a1)
a2=a2/norm(a2)
W=Pw_noise*(a1+sqrt(-1)*a2)
X=X+W
% //--------------------------SMI算法----------------------------------------
Rd=X*S0'/SP
R=X*X'/(SP*1)
Wc_SMI=pinv(R)*Rd/(Rd'*pinv(R)*Rd)%权向量
Wc_SMI=Wc_SMI/norm(Wc_SMI)
Y_SMI=Wc_SMI'*X %SMI算法恢复出来的信号
%//-------------------------------------GS算法------------------
m=1
for i=1:400:2000
X2(:,m)=XX2(:,i)
m=m+1
end
a=zeros(1,8)
phi_n=zeros(1,8)
phi=0:pi/180:2*pi
theta=0:pi/180:pi/2
for kk=1:8
a(kk)=1
phi_n(kk)=2*pi*(kk-1)/8
end
x1=zeros(1,8)
x2=zeros(1,8)
x3=zeros(1,8)
x4=zeros(1,8)
x5=zeros(1,8)
x1=X2(:,1)'
x2=X2(:,2)'
x3=X2(:,3)'
x4=X2(:,4)'
x5=X2(:,5)'
Z1=x1
Z1_inner_product=Z1*conj(Z1)
Z1_mode=sqrt(sum(Z1_inner_product))
Y1=Z1/Z1_mode
Inner_product=sum(x2*conj(Y1))
Z2=x2-Inner_product*Y1
Z2_inner_product=sum(Z2*conj(Z2))
Z2_mode=sqrt(Z2_inner_product)
Y2=Z2/Z2_mode
Inner_product1=sum(x3*conj(Y1))
Inner_product2=sum(x3*conj(Y2))
Z3=x3-Inner_product1*Y1-Inner_product2*Y2
Z3_inner_product=sum(Z3*conj(Z3))
Z3_mode=sqrt(Z3_inner_product)
Y3=Z3/Z3_mode
Inner_product1_0=sum(x4*conj(Y1))
Inner_product2_0=sum(x4*conj(Y2))
Inner_product3_0=sum(x4*conj(Y3))
Z4=x4-Inner_product1_0*Y1-Inner_product2_0*Y2-Inner_product3_0*Y3
Z4_inner_product=sum(Z4*conj(Z4))
Z4_mode=sqrt(Z4_inner_product)
Y4=Z4/Z4_mode
Inner_product1_1=sum(x5*conj(Y1))
Inner_product2_1=sum(x5*conj(Y2))
Inner_product3_1=sum(x5*conj(Y3))
Inner_product4_1=sum(x5*conj(Y4))
Z5=x5-Inner_product1_1*Y1-Inner_product2_1*Y2-Inner_product3_1*Y3-Inner_product4_1*Y4
Z5_inner_product=sum(Z5*conj(Z5))
Z5_mode=sqrt(Z5_inner_product)
Y5=Z5/Z5_mode
%Y1
%Y2
%Y3
%Y4
%Y5
w0=zeros(1,8)
w=zeros(1,8)
for mm=1:8
w0(mm)=exp(-j*2*pi*R0*cos(d6-phi_n(mm))*sin(e6))
end
dd1=sum(w0*conj(Y1))*Y1
dd2=sum(w0*conj(Y2))*Y2
dd3=sum(w0*conj(Y3))*Y3
dd4=sum(w0*conj(Y4))*Y4
dd5=sum(w0*conj(Y5))*Y5
w=w0-dd1-dd2-dd3-dd4-dd5
Wc_GS=w
Wc_GS=Wc_GS/(norm(Wc_GS))
Y_GS=Wc_GS*X %GS算法恢复出来的图像
%//----------------------------------MMSE算法-----------------------
Rd=X*S0'/SP
R=X*X'/(SP*1)
Wc_MMSE=pinv(R)*Rd
Wc_MMSE=Wc_MMSE/norm(Wc_MMSE)
Y_MMSE=Wc_MMSE'*X %MMSE算法恢复出来的信号
S0=S0/norm(S0)
Y_GS=Y_GS/norm(Y_GS)
Y_SMI=Y_SMI/norm(Y_SMI)
Y_MMSE=Y_MMSE/norm(Y_MMSE)
% figure(1)
% plot(real(S0))
% title('原始信号')
% xlabel('采样快拍数')
% ylabel('信号幅度')
% figure(2)
% plot(real(Y_SMI))
% title('运用SMI算法处理出的信号')
% xlabel('采样快拍数')
% ylabel('信号幅度')
% figure(3)
% plot(real(Y_GS))
% title('运用G-S算法处理出的信号')
% xlabel('采样快拍数')
% ylabel('信号幅度')
% figure(4)
% plot(real(Y_MMSE))
% for i=1:SP
% ss(i)=abs(S0(i)-Y_SMI(i))^2
% end
% q_1=mean(ss)
% for i=1:SP
% ss1(i)=abs(S0(i)-Y_GS(i))^2
% end
% q_2=mean(ss1)
% for i=1:SP
% ss2(i)=abs(S0(i)-Y_MMSE(i))^2
% end
% q_3=mean(ss2)
%
% 1=1+q_1
% 2=2+q_2
% 3=3+q_3
% end
% 1/16000
% 2/16000
% 3/16000
phi=0:pi/180:2*pi
theta=0:pi/180:pi/2
%
% % //------------------------ 形成波束-----------------------------------------
F_mmse=zeros(91,361)
F_smi=zeros(91,361)
F_gs=zeros(91,361)
for mm=1:91
for nn=1:361
p1=sin(theta(mm))
p2=cos(phi(nn))
p3=sin(phi(nn))
q1=sin(e6)
q2=cos(d6)
q3=sin(d6)
for hh=1:8
w1=cos(phi_n(hh))
w2=sin(phi_n(hh))
zz1=q2*w1+q3*w2
zz2=p2*w1+p3*w2
zz=zz2*p1-zz1*q1
F_mmse(mm,nn)= F_mmse(mm,nn)+conj(Wc_MMSE(hh))*(exp(j*2*pi*R0*(zz2*p1)))
F_smi(mm,nn)=F_smi(mm,nn)+conj(Wc_SMI(hh))*(exp(j*2*pi*R0*(zz2*p1)))
F_gs(mm,nn)=F_gs(mm,nn)+conj((Wc_GS(hh))')*(exp(j*2*pi*R0*(zz2*p1)))
end
end
end
F_MMSE=abs(F_mmse)
F_SMI=abs(F_smi)
F_GS=abs(F_gs)
figure(5)
mesh(20*log10(F_MMSE))
figure(6)
mesh(20*log10(F_SMI))
title('SMI算法波束形成图')
xlabel('方位角')
ylabel('俯仰角')
zlabel('幅度/dB')
figure(7)
mesh(20*log10(F_GS))
title('G-S算法波束形成图')
xlabel('方位角')
ylabel('俯仰角')
zlabel('幅度/dB')
SNR=0A=10^(SNR/10)
C=zeros(9,9)
sum=zeros(9,9)
for Nt=1:9
for Nr=1:9
if (Nt>=Nr)
Im=eye(Nr)
else Im=eye(Nt)
end
%sum=0
H=randn(Nr,Nt)
if(Nr<Nt)
Q=H*H'
else Q= H'*H
end
sum(Nt,Nr)= log2(det(Im+A*Q/Nt))
C(Nt,Nr)=sum(Nt,Nr)
end
end
C
按上面语此袜逗句森卖试试好晌
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)