浅水声信道模型的建立（4）----只考虑海面海底一次散射，一个目标研究

我思考了一下上一个程序为什么会失败，可能还是太远了，所以这次弄近一点。把5000米改成1000米吧。。。其他的还是不变

function的持续不变，在这里就不重复写了，具体参照上一篇博客

主函数：关于亮点1

%================================================================================

clc;close all;clear all;

%参考：数字信号处理：page:59p;

%=============================================

% r=0;z=10;%发射信号的位置。

% r1=1000;z1=15;%1号亮点

% r2=1010;z2=15;%2号亮点

% r3=1050;z3=15;%3号亮点

% rr=0;zz=20:30;%接收信号

c = 1500;            %声速             Unit：m/s

SNR = 60;            %信噪比           Unit：dB H = 100;       %水深                    Unit：m

%Sample_time = 0.1;   %采样时长         Unit：s      假设信号发射的时刻为零时刻

Re_coef_surf = -1;      %水面反射系数 

Re_coef_bottom = 0.8;   %水底反射系数 

Reflex_num =1;         %考虑最大的反射次数

%=========================================

f0=1600;           %信号频率         Unit：Hz

fs=16000;          %采样率          Unit：Hz      

T=5/f0;              %信号脉宽         Unit：s %0.0031

ts=0:1/fs:T-1/fs;   %一个脉冲包络的持续时间序列

Ns=length(ts);    %50

startt=0.5;endt=3.5;%混响持续时间

Sample_time=endt;

azm=pi/6;

%============================亮点1===========================================
%第一过程:发射换能器到亮点1   % r=0;z=10;%发射信号的位置。
H1 = 10;       %发射点水深           Unit：m
H2 = 15;       %接收点水深           Unit：m
D = 1000;       %接收与发射水平距离  Unit：m
%============================================================================
signal_11=Com_Sonar_data_generation_2(c,fs,f0,T,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);
signal_11_fft=fft(signal_11);
X1=abs(fft(signal_11));
X1=X1/max(X1);

%画图
figure(1);
plot(signal_11);
xlabel('时间s');ylabel('幅度');
string = ['信号发生器—亮点1，信噪比SNR=',num2str(SNR),'dB,考虑的最大反射次数',num2str(Reflex_num),'时域'];title(string);

%---------------------------------------------------------------------------------------------------------------
figure(2);
plot(1:length(signal_11_fft),signal_11_fft);
xlabel('频率omega');ylabel('幅度');
string = ['信号发生器—亮点1，信噪比SNR=',num2str(SNR),'dB,考虑的最大反射次数',num2str(Reflex_num),'频域'];title(string);

%--------------------------------------------------------------------------------------------------------------
figure(3);
N=length(signal_11);
f=1:N/2;   
plot(f,signal_11_fft(1:N/2));
xlabel('频率/Hz');ylabel('归一化幅值');title('频谱的一半');grid on;

%==========================================================

%过程2：从亮点散射

Np=(endt-startt)*fs; %散射系数的点数48000个点，发射信号持续时间内采样的点

for th=startt:1/fs:endt-1/fs   %混响时间 Np48000个点,th共循环48000次    

r1=c*(th-T)/2;   %T=0.005    

r2=c*th/2;     %c=1500    

N=round((r2^2-r1^2)*azm/2); %azm=0.5236（pi/6）方位角     %随着混响时间的增加，混响点数也在增加    

ain=randn(1,N);     %幅度    

fai=randn(1,N)*2*pi;  %相位      

for i=1:N        

p1(i)=ain(i)*exp(1i*fai(i));       

end    

P1(round((th-startt)*fs)+1)=sum(p1)/r1^2;   %48000

end

% 得到混响信号

R1=conv(signal_11,P1); % length(signal_11)+length(P1)-1;  56000+48000-1=103999 

Rr1=real(R1);

%==================================================================== 

figure(4);
t=startt:1/fs:endt+1240*T-6002/fs;      %: 共103999个点,这个是凑出来的
plot(t,Rr1/max(Rr1));
axis([startt endt+1240*T-6002/fs -1 1]);xlabel('时间/s');ylabel('归一化幅值');title('亮点1处--混响');

%---------------------------------------------------------------------------------------------------------------------------
figure(5);
R1_fft=fft(R1);
plot(1:length(R1_fft),R1_fft);
xlabel('omega/pi');ylabel('归一化幅值|e^j^omega|');title('亮点1处--混响---相关后傅里叶');

%-------------------------------------------------------------------------------------------------------------------------
figure(6);
R1_fft=fft(R1);
plot(1:length(R1_fft)/2,R1_fft(1:length(R1_fft)/2));
xlabel('omega/pi');ylabel('归一化幅值|e^j^omega|');title('亮点1处--混响---相关后傅里叶--频谱的一半');

%=================================================================================

%重新思考第三过程： %选取阵列第一个点为参考点，zz=20。zz=20:30,rr=0; %亮点1散射后到达阵列1  r1=5000;z1=15;%1号亮点   rr=0;zz=20:30;%接收信号

H1 = 15;       %发射点水深           Unit：m

D = 1000;       %接收与发射水平距离  Unit：m

for kkk=1:10;       %接收点水深           Unit：m        

H2=20+1*kkk;        

signal_13(kkk,:)=Com_Sonar_data_generation_2(c,fs,f0,T,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);

end    

for kkk=1:10    

zhenlie_1(kkk,:)=conv(signal_13(kkk,:),R1);

end   

total_1=zhenlie_1(1,:);

for kkk=2:10    

total_1=zhenlie_1(kkk,:)+total_1;

end

%==========================================================

figure(7);
plot(1:length(total_1),total_1);
xlabel('时间‘s');ylabel('幅度');title('阵列接收到的亮点1的信号');

% =========================================
c = 1500;            %声速             Unit：m/s
D=1;                %阵元的距离
M = 10;             % 阵元数
m=[0:M-1];          %阵元序号
d=0.2;           %系统实孔径0.2
T=5/f0;              %信号脉宽         Unit：s %0.0031
k=d/(c*T);
doa=atan(5/1000);                                          %doa为线阵法线方向与声源的夹角
w = exp(j*2*pi*k*m'*sin(doa*pi/180));    %每个阵元的时延tow=k*m‘*sin（）不同  
temp=zeros(M,length(total_1));
for kkk=1:10
    temp(kkk,1)=w(kkk,1);
end
ww=temp*total_1';     

theta = linspace(-90,90,360);                        % 扫描方位角 a = exp(-j*2*pi*k*m'*sin(theta*pi/180));           %theta是扫描角度% 导向适量

Y = abs(ww'*a); Y = 20*log10(Y/max(Y)); Y_max=max(Y);

figure(8);

plot(theta,Y);

title('未加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');

%=============================================================

%加权

quanzhi=[zeros(1,90) 0.5*ones(1,45) 1*ones(1,90) 0.5*ones(1,45) zeros(1,90)];

tempp=zeros(360,360);

for kkk=1:360    

tempp(1,kkk)=quanzhi(1,kkk);

end

YY = abs(Y*tempp);    

YY = 20*log10(YY/max(YY));

figure(9);

plot(theta,YY);

title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');

%============================================================================

YYY=2.^(3*Y);

YYY = 20*log10(YYY/max(YYY));

YYY_max=max(YYY);

figure(10);

plot(theta,YYY);

title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)'); 

%=================================================================

 最后运行其他两个亮点的程序，得到：

还是分不开。。。。。。。。。。。。。。。。。。。。。。。。。。

%=============================================================

%程序已结束===========================================

%感觉有错误，我重新思考下=========================================

我重新改了一下阵列处理那块，还是按5000m来处理，前面两个过程的程序没有改变，贴上第三个过程的程序：

%这次改进的是阵列与散射信号的运算过程，原来是相乘，我改成了卷积。
%bishe16的基础上进行修改

%=================================================================================

%重新思考第三过程：

%选取阵列第一个点为参考点，zz=20。zz=20:30,rr=0;

%亮点1散射后到达阵列1  r1=5000;z1=15;%1号亮点   rr=0;zz=20:30;%接收信号

H1 = 15;       %发射点水深           Unit：m

D = 5000;       %接收与发射水平距离  Unit：m

for kkk=1:10;       %接收点水深           Unit：m        

H2=20+1*kkk;        

signal_13(kkk,:)=Com_Sonar_data_generation_2(c,fs,f0,T,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);

end    

for kkk=1:10    

zhenlie_1(kkk,:)=conv(signal_13(kkk,:),R1);

end   

total_1=zhenlie_1(1,:);  

for kkk=2:10    

total_1=zhenlie_1(kkk,:)+total_1;  %1x159998

end

figure(7);

plot(1:length(total_1),total_1);

xlabel('时间‘s');ylabel('幅度');title('阵列接收到的亮点1的信号');

% =========================================

c = 1500;            %声速             Unit：m/s

D=1;                %阵元的距离 M = 10;             % 阵元数

m=[0:M-1];          %阵元序号

d=0.2;           %系统实孔径0.2

T=5/f0;              %信号脉宽         Unit：s %0.0031

k=d/(c*T);

doa=atan(5/5000);                                          %doa为线阵法线方向与声源的夹角

w = exp(j*2*pi*k*m'*sin(doa*pi/180));    %每个阵元的时延tow=k*m‘*sin（）不同  

theta = linspace(-90,90,360);                       

% 扫描方位角

a = exp(-j*2*pi*k*m'*sin(theta*pi/180));%theta是扫描角度% 导向适量  10x360矩阵 a=a';    %360x10矩阵  %360x10矩阵

temp=zeros(length(theta),length(ww));   %360x160007juzhen

for kkk=1:10    

temp(:,kkk)=a(:,kkk);

end

temp=temp';

Y = abs(ww*temp);  %1x160007 *160007x360 

Y = 20*log10(Y/max(Y));

Y_max=max(Y);

figure(8);

plot(theta,Y);

title('未加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');

%=============================================================

%加权

quanzhi=[zeros(1,90) 0.5*ones(1,45) 1*ones(1,90) 0.5*ones(1,45) zeros(1,90)];

tempp=zeros(360,360); for kkk=1:360    

tempp(1,kkk)=quanzhi(1,kkk); end YY = abs(Y*tempp);    

YY = 20*log10(YY/max(YY));

figure(9);

plot(theta,YY);

title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');

%============================================================================

YYY=2.^(3*Y); YYY = 20*log10(YYY/max(YYY)); YYY_max=max(YYY);

figure(10);

plot(theta,YYY);

title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)'); 

这肯定是不对的。。。。。。。。。。。。。。。。。。。。。。。。。

 %===================================================

 %距离5000改成1000m

%我思考了一下，为什么非要得到一个角度，来和得到的信号进行运算？我得到的信号不是自带角度了吗
%再改进一下

%=================================================================================

%重新思考第三过程： %选取阵列第一个点为参考点，zz=20。zz=20:30,rr=0; %亮点1散射后到达阵列1  r1=5000;z1=15;%1号亮点   rr=0;zz=20:30;%接收信号

H1 = 15;       %发射点水深           Unit：m

D = 1000;       %接收与发射水平距离  Unit：m

for kkk=1:10;       %接收点水深           Unit：m        

H2=20+1*kkk;        

signal_13(kkk,:)=Com_Sonar_data_generation_2(c,fs,f0,T,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);

end    

for kkk=1:10    

zhenlie_1(kkk,:)=conv(signal_13(kkk,:),R1);

end   

total_1=zhenlie_1(1,:);

for kkk=2:10    

total_1=zhenlie_1(kkk,:)+total_1;

end

figure(7);

plot(1:length(total_1),total_1);

xlabel('时间‘s');ylabel('幅度');title('阵列接收到的亮点1的信号');

% =========================================

c = 1500;            %声速             Unit：m/s

D=1;                %阵元的距离

M = 10;             % 阵元数 m=[0:M-1];          %阵元序号

d=0.2;           %系统实孔径0.2

T=5/f0;              %信号脉宽         Unit：s %0.0031

k=d/(c*T);

%zhenlie_1是10个阵元分别接受到的信号，10*159998的矩阵 %total_1是10个阵元的信号全部加在了一起 %先用zhenlie_1来处理程序看看结果

theta = linspace(-90,90,360);                        % 扫描方位角

a = exp(-j*2*pi*k*m'*sin(theta*pi/180));           %theta是扫描角度% 导向适量  10x360

zhenlie_1=zhenlie_1';    %159998x10矩阵

Y = abs(zhenlie_1*a);     %159998x360 %Y得加成1x360的矩阵

temp=Y(1,:);

for kkk=2:length(total_1);    

temp=Y(kkk,:)+temp;

end

Y=temp;

Y = 20*log10(Y/max(Y));

Y_max=max(Y);

figure(8);

plot(theta,Y);

title('未加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');

%=============================================================

%加权

quanzhi=[zeros(1,90) 0.5*ones(1,45) 1*ones(1,90) 0.5*ones(1,45) zeros(1,90)];

tempp=zeros(360,360); for kkk=1:360    

tempp(1,kkk)=quanzhi(1,kkk);

end

YY = abs(Y*tempp);    

YY = 20*log10(YY/max(YY));

figure(9);

plot(theta,YY);

title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');

%============================================================================

YYY=2.^(3*Y);

YYY = 20*log10(YYY/max(YYY));

YYY_max=max(YYY);

figure(10);

plot(theta,YYY);

title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)'); 

Y_max=0,所以这个思路可能也是错的，所以我把距离改更小，100m来试验一下是程序的原因还是距离太远的原因。

%============================================================================================

clc;close all;clear all; %参考：数字信号处理：page:59p;

%=============================================

% r=0;z=10;%发射信号的位置。 % r1=100;z1=15;%1号亮点% r2=110;z2=15;%2号亮点% r3=150;z3=15;%3号亮点 % rr=0;zz=20:30;%接收信号

c = 1500;            %声速             Unit：m/s

SNR = 60;            %信噪比           Unit：dB

H = 100;       %水深                    Unit：m

%Sample_time = 0.1;   %采样时长         Unit：s      假设信号发射的时刻为零时刻

Re_coef_surf = -1;      %水面反射系数 

Re_coef_bottom = 0.8;   %水底反射系数 

Reflex_num =1;         %考虑最大的反射次数

%=========================================

f0=1600;           %信号频率         Unit：Hz

fs=16000;          %采样率          Unit：Hz      

T=5/f0;              %信号脉宽         Unit：s %0.0031

ts=0:1/fs:T-1/fs;   %一个脉冲包络的持续时间序列 Ns=length(ts);    %50

startt=0.5;endt=3.5;%混响持续时间

Sample_time=endt;

azm=pi/6;            

%============================亮点1===========================================

%第一过程:发射换能器到亮点1   % r=0;z=10;%发射信号的位置。

H1 = 10;       %发射点水深           Unit：m

H2 = 15;       %接收点水深           Unit：m

D = 100;       %接收与发射水平距离  Unit：m

%============================================================================

signal_11=Com_Sonar_data_generation_2(c,fs,f0,T,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);

signal_11_fft=fft(signal_11);

X1=abs(fft(signal_11));

X1=X1/max(X1); %画图

figure(1);

plot(signal_11);

xlabel('时间s');ylabel('幅度'); string = ['信号发生器—亮点1，信噪比SNR=',num2str(SNR),'dB,考虑的最大反射次数',num2str(Reflex_num),'时域'];title(string);

figure(2);

plot(1:length(signal_11_fft),signal_11_fft);

xlabel('频率omega');ylabel('幅度'); string = ['信号发生器—亮点1，信噪比SNR=',num2str(SNR),'dB,考虑的最大反射次数',num2str(Reflex_num),'频域'];title(string);

figure(3);

N=length(signal_11); f=1:N/2;   

plot(f,signal_11_fft(1:N/2)); xlabel('频率/Hz');ylabel('归一化幅值');title('频谱的一半');grid on;

%==========================================================

%过程2：从亮点散射

Np=(endt-startt)*fs; %散射系数的点数48000个点，发射信号持续时间内采样的点 for th=startt:1/fs:endt-1/fs   %混响时间 Np48000个点,th共循环48000次    

r1=c*(th-T)/2;   %T=0.005    

r2=c*th/2;     %c=1500    

N=round((r2^2-r1^2)*azm/2); %azm=0.5236（pi/6）方位角     %随着混响时间的增加，混响点数也在增加    

ain=randn(1,N);     %幅度    

fai=randn(1,N)*2*pi;  %相位      

for i=1:N       

  p1(i)=ain(i)*exp(1i*fai(i));   

    end    

P1(round((th-startt)*fs)+1)=sum(p1)/r1^2;   %48000

end

% 得到混响信号

R1=conv(signal_11,P1); % length(signal_11)+length(P1)-1;  56000+48000-1=103999 

Rr1=real(R1);

%==================================================================== 

figure(4);

t=startt:1/fs:endt+1240*T-6002/fs;      %: 共103999个点,这个是凑出来的

plot(t,Rr1/max(Rr1));

axis([startt endt+1240*T-6002/fs -1 1]);xlabel('时间/s');ylabel('归一化幅值');title('亮点1处--混响');

figure(5);

R1_fft=fft(R1);

plot(1:length(R1_fft),R1_fft);

xlabel('omega/pi');ylabel('归一化幅值|e^j^omega|');title('亮点1处--混响---相关后傅里叶');

figure(6);

R1_fft=fft(R1);

plot(1:length(R1_fft)/2,R1_fft(1:length(R1_fft)/2)); xlabel('omega/pi');ylabel('归一化幅值|e^j^omega|');title('亮点1处--混响---相关后傅里叶--频谱的一半');

%=================================================================================

%重新思考第三过程：

%选取阵列第一个点为参考点，zz=20。zz=20:30,rr=0; %亮点1散射后到达阵列1  r1=5000;z1=15;%1号亮点   rr=0;zz=20:30;%接收信号

H1 = 15;       %发射点水深           Unit：m

D = 100;       %接收与发射水平距离  Unit：m

for kkk=1:10;       %接收点水深           Unit：m       

  H2=20+1*kkk;        

signal_13(kkk,:)=Com_Sonar_data_generation_2(c,fs,f0,T,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);

end    

for kkk=1:10   

  zhenlie_1(kkk,:)=conv(signal_13(kkk,:),R1);

end   

total_1=zhenlie_1(1,:);

for kkk=2:10    

total_1=zhenlie_1(kkk,:)+total_1;

end

figure(7);

plot(1:length(total_1),total_1);

xlabel('时间‘s');ylabel('幅度');title('阵列接收到的亮点1的信号');

% =========================================

c = 1500;            %声速             Unit：m/s

D=1;                %阵元的距离

M = 10;             % 阵元数

m=[0:M-1];          %阵元序号

d=0.2;           %系统实孔径0.2

T=5/f0;              %信号脉宽         Unit：s %0.0031

k=d/(c*T);

%zhenlie_1是10个阵元分别接受到的信号，10*159998的矩阵

%total_1是10个阵元的信号全部加在了一起

%先用zhenlie_1来处理程序看看结果

theta = linspace(-90,90,360);                        % 扫描方位角

a = exp(-j*2*pi*k*m'*sin(theta*pi/180));           %theta是扫描角度% 导向适量  10x360

zhenlie_1=zhenlie_1';    %159998x10矩阵

Y = abs(zhenlie_1*a);     %159998x360 %Y得加成1x360的矩阵

temp=Y(1,:);

for kkk=2:length(total_1);    

temp=Y(kkk,:)+temp;

end

Y=temp;

Y = 20*log10(Y/max(Y));

Y_max=max(Y);

figure(8);

plot(theta,Y);

title('未加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');

%=============================================================

%加权

quanzhi=[zeros(1,90) 0.5*ones(1,45) 1*ones(1,90) 0.5*ones(1,45) zeros(1,90)];

tempp=zeros(360,360);

for kkk=1:360    

tempp(1,kkk)=quanzhi(1,kkk);

end

YY = abs(Y*tempp);   

  YY = 20*log10(YY/max(YY));

figure(9);

plot(theta,YY);

title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');

%============================================================================

YYY=2.^(3*Y);

YYY = 20*log10(YYY/max(YYY));

YYY_max=max(YYY);

figure(10);

plot(theta,YYY);

title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)'); 

所以说这个思路还是错的。。。。。。。。。。。。。。。。。。。

%===================================================================

%接下来那直接把阵列信号全部加在一起，再进行处理

%m=100米；

%=================================================================================

%重新思考第三过程： %选取阵列第一个点为参考点，zz=20。zz=20:30,rr=0; %亮点1散射后到达阵列1  r1=5000;z1=15;%1号亮点   rr=0;zz=20:30;%接收信号

H1 = 15;       %发射点水深           Unit：m

D = 100;       %接收与发射水平距离  Unit：m

for kkk=1:10;       %接收点水深           Unit：m        

H2=20+1*kkk;        

signal_13(kkk,:)=Com_Sonar_data_generation_2(c,fs,f0,T,Sample_time,SNR,H,H1,H2,D,Re_coef_surf,Re_coef_bottom,Reflex_num);

end    

for kkk=1:10    

zhenlie_1(kkk,:)=conv(signal_13(kkk,:),R1);

end  

  total_1=zhenlie_1(1,:);

for kkk=2:10    

total_1=zhenlie_1(kkk,:)+total_1;

end

figure(7);

plot(1:length(total_1),total_1);

xlabel('时间‘s');ylabel('幅度');title('阵列接收到的亮点1的信号');

% =========================================

c = 1500;            %声速             Unit：m/s

D=1;                %阵元的距离 M = 10;             % 阵元数

m=[0:M-1];          %阵元序号

d=0.2;           %系统实孔径0.2

T=5/f0;              %信号脉宽         Unit：s %0.0031

k=d/(c*T);

%zhenlie_1是10个阵元分别接受到的信号，10*159998的矩阵

%total_1是10个阵元的信号全部加在了一起

%用total_1来处理

theta = linspace(-90,90,360);                        % 扫描方位角

a = exp(-j*2*pi*k*m'*sin(theta*pi/180));           %theta是扫描角度% 导向适量 10x360

total_1=total_1';      %1x159998矩阵

temp=zeros(length(total_1),length(theta));   %159998x360

for kkk=1:10  

   temp(kkk,:)=a(kkk,:);

end

Y = abs(total_1'*temp);     %1x159998 *159998x360

Y = 20*log10(Y/max(Y));

Y_max=max(Y);

figure(8);

plot(theta,Y);

title('未加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');

%=============================================================

%加权

quanzhi=[zeros(1,90) 0.5*ones(1,45) 1*ones(1,90) 0.5*ones(1,45) zeros(1,90)];

tempp=zeros(360,360);

for kkk=1:360    

tempp(1,kkk)=quanzhi(1,kkk);

end

YY = abs(Y*tempp);    

YY = 20*log10(YY/max(YY));

figure(9);

plot(theta,YY);

title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');

%============================================================================

YYY=2.^(3*Y);

YYY = 20*log10(YYY/max(YYY));

YYY_max=max(YYY);

figure(10);

plot(theta,YYY);

title('加权方向图'); axis([-90,90,-50,0]); grid on; xlabel('方位角/(circ)');ylabel( 'P/(dB)');