近代天线理论课程设计-等副瓣天线

近代天线理论课程设计

一、课程设计题目要求：

利用Matlab编制遗传算法或粒子群算法程序，并实现对间距为半波长均匀直线阵综合，要求如下：

1、阵元数：20

2、副瓣电平：-25dB等副瓣

3、要求撰写设计报告，内容包括：算法原理，目标函数的设计，源代码，仿真结果（画出增益方向图，优化各单元的幅度，并与相同单元数的切比雪夫均匀直线阵做对比），参考文献。

二、遗传算法原理介绍：

1、遗传算法的生物学基础

遗传算法的生物学基础自然选择学说认为适者生存，生物要存活下去，就必须进行生存斗争。生存斗争包括种内斗争、种间斗争以及生物跟环境之间的斗争三个方面。在生存斗争中，具有有利变异的个体容易存活下来，并且有更多的机会将有利变异传给后代，具有不利变异的个体就容易淘汰，产生后代的机会就少得多。因此，凡是在生存斗争中获胜的个体都是对和环境适应性比较强的个体。

2、遗传算法的科学定义

遗传算法（GA）是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。其主要特点是直接对结构对象进行操作，不存在求导和函数连续性的限定；具有内在的隐并行性和更好的全局寻优能力；采用概率化的寻优方法，不需要确定的规则就能自动获取和知道优化的搜索空间，自适应地调整搜索方向。遗传算法以一种群体中的所有个体为对象，并利用随机化技术指导对一个被编码的参数空间进行高效搜索。其中，选择、交叉和变异构成了遗传算法的遗传操作；参数编码、初始群体的设定、适应度函数的设计、遗传操作设计、控制参数设定五个要素组成了遗传算法的核心内容。

3、遗传算法的执行过程

遗传算法是从代表问题可能潜在的解集的一个种群（population）开始的，而一个种群则由经过基因（gene）编码的一定数目的个体（individual）组成。每个个体实际上是染色体（chromosome）带有特征的实体。

染色体作为遗传物质的主要载体，即多个基因的集合，其内部表现（即基因型）是某种基因组合，它决定了个体的形状的外部表现，如黑头发的特征是由染色体中控制这一特征的某种基因组合决定的。因此，在一开始需要实现从表现型到基因型的映射即编码工作。由于仿照基因编码的工作很复杂，我们往往进行简化，如二进制编码。

初代种群产生之后，按照适者生存和优胜劣汰的原理，逐代（generation）演化产生出越来越好的近似解，在每一代，根据问题域中个体的适应度（fitness）大小选择（selection）个体，并借助于自然遗传学的遗传算子（genetic operators）进行组合交叉（crossover）和变异（mutation），产生出代表新的解集的种群。

这个过程将导致种群像自然进化一样的后生代种群比前代更加适应于环境，末代种群中的最优个体经过解码（decoding），可以作为问题近似最优解。

三、利用遗传算法综合天线

1、具体流程如下

第一步:随机产生初始种群。种群规模由需要优化的参量(例如:单元间距,馈电幅度和相位)确定,一般种群数目应取优化参数个数的十几倍。比如优化20元阵列的每个单元的馈电相位,则种群数目取200~400之间。同时确定每个单元相位的二进制编码长度,此长度需要在精度和计算时间中权衡。在本次的编码中，选取了一个种群中有50个个体，并且每个个体中具有20个阵单元.

第二步:计算每个个体的适应度值。对以上随机产生的二进制相位进行解码,生成十进制相位,代入公式计算阵列方向函数,与目标方向函数做差,并且根据不同方向的重要性对差分别加权,取每个个体中不同方向差的绝对值和 ,然后取倒数得到每个个体的适应度值。这里的适应度值越大就是跟目标低副瓣插值的求和越小。

第三步:选择。根据一定的算法,配合适应度值选择再生个体。比如可采用博彩活动中常用的赌轮盘法进行选择,这样适应度高的个体生存的可能性大,也就是说接近目标方向函数的阵列保留下来的概率大。

第四步:交叉。为了使下一代种群遗传上一代的基因,按照一定的交叉概率和交叉算法 ,生成新的个体。选择适应度最大的电流对应的编码值，然后分别与每一个个体进行交叉产生后代，分别生出三个后代，分别占比是1：1、2：1、1：2，然后分别计算适应度值，取最好的子代留下。Ii 是当前电流值的二进制，Ibest 是当前个体中适应度最高的对应的电流值的二进制。

第五步:变异。为了避免近亲繁殖,陷入局部最优,按照一定的变异概率和变异算法,生成新的个体。

第六步:判断是否达到最大遗传代数或者适应度值达到要求。若达到,则输出最佳值,否则返回第二步,继续这个过程。

2、目标函数的设计

对应着上述的第二步。首先寻找一种对问题潜在解进行“数字化”编码的方案（建立表现型和基因型的映射关系）。然后用随机数初始化一个种群，种群里面的个体就是这些数字化的编码。接下来，通过适当的解码过程之后（得到种群个体的位置坐标），用适应性函数对每一个基因个体作一次适应度评估（副瓣电平越小，适应度越高）。用选择函数按照某种规定择优选择。

在本次的目标函数的设计中，将各个副瓣的值与-25db进行均方误差，公式如下：

这里的peaks是所有峰值的取值，包括主瓣，所以在最后有一个减去25*25。这样计算的；就是所有副瓣与-25db的差值平方，采用均方误差的办法进行讨论。

3、遗传算法的流程图：

  

  对流程图的讲解，这里染色体初始化就是随机产生（0~1）之间的二进制的数值，编码数是72。适应度检测是计算当前的目标函数值是副瓣电平与25db之间的差值，看当前是否满足最小适应度。配对繁殖是父代产生子代的过程，选择当前最优的电流值跟别的电流进行繁殖，然后选出最优子代进行后续操作，而且最优父代也会保留。基因突变是随机产生一个数，看是否小于一个值，如果小于就是随机电流值的二进制某一个部位进行变化。设置的迭代次数达到之后或者适应度满足一定要求便会结束。

四、切比雪夫天线阵列

1、切比雪夫线阵计算原理

定义切比雪夫多项式：

Tnx=cos⁡(n*arccos⁡(x))

从关于阵列综合的研究中可以看到，随着电流幅度自线阵中心向两端展开，副瓣降低而主瓣展宽;但压低甚至消除副瓣将带来更多问题。在大多数应用中，既希望主瓣窄，又希望副瓣低。因此需要在主瓣宽度与副瓣电平两项指标之间实现平衡，即找到最佳折中:对于给定的主瓣宽度，副瓣尽可能低;对于给定的副瓣电平，主瓣尽可能窄。对于这种方向图就称为最佳方向图。而研究已经表明，假设副瓣所在的空域正是我们要抑制其辐射的空域，并且抑制辐射对所有副瓣都是同等重要的，则最佳设计就是使所有副瓣电平都相等。真正具有实际应用价值的阵列综合方法是实现或者逼近最佳方向图。

而对于此，若一线阵的波瓣图与某阶的切比雪夫多项式曲线相一致,则该天线阵就会具有最佳的方向性。因为根据这种多项式模拟的旁瓣的幅度都均相等,所以只要靠近主瓣的旁瓣电平不超过给定值即可。如果阵列天线单元上激励电流的幅度分布可以使得获得的阵因子函数符合于切比雪夫多项式,则其上的电流分布就为切比雪夫分布。

对于一个N单元的Chebyshev阵列的激励幅度的确定，是通过Dolph变换将其阵因子F(u)与Chebyshev多项式Tn-1(x) 等同，其中 Dolph变换为

x=b*cos⁡(u/2)

其中，u=kd*cos⁡(φ-α) ,b为一个大于1的参数，即b>1

为确定一个对称激励的N单元阵列的激励幅度，可将坐标的原点放在阵列的中心并将阵因子写为：

其中，u=kd*cos⁡(φ-α) ,α 是相邻单元的相位差，am 为激励幅度。用Dolph变换u=2*arccos⁡(xb) 并将此式带入上式，得

将上式子级数展开可得到

进而激励幅度am 的确定就可通过比较关系式Fu=TN-1b*cosu2=TN-1(x) 两端同幂项的系数得到，再将am 带入阵因子的式子可得到阵因子并绘制方向图。

2、切比雪夫线阵

当|x|>1时，切比雪夫多项式的值Tm(x)随x的增大而无限增大，故需要对x值加以限制:在|x|>1的情况下，当x=x0 ,时令Tm(x)=R。因为n元线阵的阵函数与m=n-1阶的切比雪夫多项式相对应，故在阵函数曲线的意义上，R即是主瓣最大值与副瓣最大值之比。所以给定副瓣电平后，即可求得R值，由R及阵元数n即可求得m=n-1阶切比雪夫多项式的x0 。切比雪夫多项式中相对于主瓣最大值R的x值是大于1的x0 值。我们要将阵函数和多项式之间的横坐标统一起来，则在阵函数有最大值的横坐标为1时，与其对应的切比雪夫多项式的横坐标为x0，因此当阵函数的坐标为任意x值时(|x|≤1），其相应切比雪夫多项式的横坐标为x0*x 。

根据分析可求得切比雪夫天线阵的电流加权值:

对于天线阵的辅射特性来说，起主要作用的是阵函数而不是天线本身的辐射特性。因此给定天线的阵元数，以及要求的副瓣电平即可求得相应的阵列电流加权值对电流进行加权。

五、结果展示并分析：

1、20单元均匀直线阵的增益方向图

结果图为：

20个单元线阵对应的电流值[0.1719 0.2435 0.1980 0.1980 0.3753 0.3260 0.3187 0.5381 0.4221 0.5451 0.4521 0.5012 0.4572 0.4803 0.4429 0.3431 0.3160 0.2679 0.2498 0.1766]

2、切比雪夫均匀直线阵的增益方向图

结果图为：

20个单元对应的电流为[0.567  0.371  0.474  0.579  0.682  0.778  0.862  0.929  0.976  1.000  1.000  0.976  0.929  0.862  0.778  0.682  0.579  0.474  0.371  0.567]

3、二者结果图比较分析

    自己写的代码过程画出的增益图不光滑，没有用切比雪夫法画的规整，但是好在实现了等副瓣的图形。

六、代码清单

1. 基本天线实现代码

clear;
clc;
close all;

%% 基本电磁变量设置
Frequency = 3e8;
Lightspeed = 3e8;
lamuda = Lightspeed/Frequency;
Wavenumber = 2*pi/lamuda;

%% 数组参数
N =20;     %阵列数
X = (1:N)*lamuda/2; %位置
I =  ones(1,N); %电流量
alpha = zeros(1,N); % 相位设为0

%% 数组因子
Ns =500;% 采样数
theta = linspace(-90,90,Ns);
E =zeros(1,Ns);

for num = 1:Ns
    E(num)=sum(I.*exp(1j*(Wavenumber*X*sind(theta(num))+alpha)))+1e-3;
    abs_E = abs(E)
end
abs_E
%% plot figure
figure()
plot(theta,db(abs_E)-max(db(abs_E)),'LineWidth',2);%标准化
xlabel('单元数N是20,半波长')
grid on

2、均匀直线阵代码

clear all;
close all; 
clc;

%%
%基本量设置
NP = 50;%种群初始个体数
f=3e9;c=3e8; lamda = c/f;d = lamda/2;NL = 20; %阵元数
wavenum = 2*pi/lamda;%波数
G = 200;L = 72; %G:最大遗传代数 L:电流值对应的二进制数
NS = 100;theta = linspace(-91,91,NS);%NS为采样数
alpha = zeros(1,NL);X = (0:NL-1)*lamda/2;
Pc = 0.8;Pm = 0.02; %Pc是交叉概率，Pm是变异概率
I = zeros(NP,NL);E = zeros(NP,NS);
error = [];Num = 1;

%%
%初始化种群个体的情况
pop = randi([0,1],NP,NL*L);  %种群个体对应的电流量的二进制编码
pop_1 = pop;
%对应种群个体的电流量
for i = 1:L
    Decode(i) = 2^(-i);%用于解码，变为相应的电流量
end
%将对应的电流二进制数变换成十进制
for j = 1:NP
    a = 1;
    for i = 1:L:NL*L  
        I_act = sum(Decode.*pop(j,i:(i+L-1)));%对应种群中的第j个个体的第n个电流量
        I(j,a) = I_act;
        a = a + 1;
    end
end
%计算相对应的主瓣和旁瓣的电平值
for j = 1:NP
    for num = 1:NS
        E(j,num)=sum(I(j,:).*exp(1j*(wavenum*X*sind(theta(num))+alpha)))+1e-3;
    end
end
E_abs = abs(E);
top = [];hang = [];top1 = [];
for j = 1:NP
    E_abs_max(1,j) = max(E_abs(j,:));
    Edb(j,:) = db(E_abs(j,:)) - db(E_abs_max(1,j));
    top1 = findpeaks(Edb(j,:));
    top = [top,100,top1];
end
top = [top,100];
location = find(top == 100);
for i = 1:NP
    infer = location(i);
    later = location(i+1);
    E1 = top(infer+1:later-1);
    error1 = sum((E1+25).*(E1+25))- 25.*25;
    error = [error,error1];
end
error_inital = error;
a = min(error)

%%
%遗传算法的循环
for g = 1:G
    g
    %基于轮盘赌进行选择适应度高的个体进行
      %先将错误率进行归一化
      max_error = max(error_inital);
      min_error = min(error_inital);
      
      %把最好的电流值对应的二进制拿出来
      Best = [];
      good = find(error_inital == min_error);
      for i = 1:NL*L
          Best(i) = pop(good(1),i);
      end
      
      for i = 1:length(error_inital)
          error_inital(i) = (error_inital(i)-min_error)/(max_error-min_error);
      end
      error = error_inital;
      %随机产生一个数，进行优秀个体的选择
      P = 0.5;
      for i = 1:NP
          if error(i) < P
              %满足条件按进行选择
              for j = 1:NL*L
                  I_select(Num,j) = pop(i,j);
              end
              Num = Num + 1;
              hang = [hang,i];
          end
      end
      %将归一化之后的误差矩阵恢复
      for i = 1:length(error_inital)
          error_inital(i) = error(i)*(max_error-min_error)+min_error;
      end 
      for i = 1:length(error)
          error(i) = error(i)*(max(error)-min(error))+min(error);
      end 
      [m,n] = size(I_select);%获得选择出的电流对应的行数，列数
    % 进行交叉操作 
      for j = 1:2:m-1 %对应于个体
          P = rand(1,1);
          if P < Pc
              for i = 1:NL*L  %对应于电流二进制
                 replace(i) = (I_select(j,i)+2*I_select(j+1,i))/3;
                 replace1(i) = (I_select(j,i)+I_select(j+1,i))/2;
                 replace2(i) = (2*I_select(j,i)+I_select(j+1,i))/3;
              end
              %计算三个后代中比较好的两个个体作为后代
              for a = 1:NL
                  I_replace(a) = sum(Decode.*replace(L*(a-1)+1:L*a));
                  I_replace1(a) = sum(Decode.*replace1(L*(a-1)+1:L*a));
                  I_replace2(a) = sum(Decode.*replace2(L*(a-1)+1:L*a));
              end
                  %计算相应的电平值
                  for num = 1:NS
                      E_replace(num)=sum(I_replace.*exp(1j*(wavenum*X*sind(theta(num))+alpha)));
                      E_replace1(num)=sum(I_replace1.*exp(1j*(wavenum*X*sind(theta(num))+alpha)));
                      E_replace2(num)=sum(I_replace2.*exp(1j*(wavenum*X*sind(theta(num))+alpha)));
                  end
                  %计算三个后代相应的误差
                  E_replace_abs_max = max(abs(E_replace));
                  Edb_replace = db(abs(E_replace)) - db(E_replace_abs_max);
                  [top_replace,location] = findpeaks(Edb_replace);
                  error_replace = sum((top_replace+25).*(top_replace+25))-25.*25;
                  E_replace1_abs_max = max(abs(E_replace1));
                  Edb_replace1 = db(abs(E_replace1)) - db(E_replace1_abs_max);
                  [top_replace1,location1] = findpeaks(Edb_replace1);
                  error_replace1 = sum((top_replace1+25).*(top_replace1+25))-25.*25;
                  E_replace2_abs_max = max(abs(E_replace2));
                  Edb_replace2 = db(abs(E_replace2)) - db(E_replace2_abs_max);
                  [top_replace2,location2] = findpeaks(Edb_replace2);
                  error_replace2 = sum((top_replace2+25).*(top_replace2+25))-25.*25;
                  descendant = [error_replace,error_replace1,error_replace2];
                  %比较后代的好坏
                  max_descendant = max(descendant);
                  min_descendant = min(descendant);
                  rr_max = find(descendant == max_descendant);
                  rr_min = find(descendant == min_descendant);
                  replace_all = [replace,replace1,replace2];
                  I_select(j,:) = replace_all(NL*L*(rr_min-1)+1:NL*L*(rr_min-1)+NL*L);
                  rr_midle = 6-rr_max-rr_min;
                  I_select(j+1,:) = replace_all(NL*L*(rr_midle-1)+1:NL*L*rr_midle);
          end     
      end
    
    % 进行变异操作
      for j = 1:m
          for i = 1:n
              varia = rand(1,1);
              if varia < Pm
                  I_select(j,i) = ~I_select(j,i);
              end
          end
      end
      
     %I_select选择的个体的交叉变异之后的电流值的二进制
     %计算交叉变异的旁瓣和主瓣电平值
     %将对应的电流二进制数变换成十进制
     for j = 1:m
         a = 1;
         for i = 1:L:NL*L  
             I_act = sum(Decode.*I_select(j,i:(i+L-1)));%对应种群中的第j个个体的第n个电流量
             I(j,a) = I_act;
             a = a + 1;
         end
     end
     %计算相对应的主瓣和旁瓣的电平值
     for j = 1:m
         for num = 1:NS
             E(j,num)=sum(I(j,:).*exp(1j*(wavenum*X*sind(theta(num))+alpha)))+1e-3;
         end
     end
     E_abs = abs(E);
     top = [];
     for j = 1:m
         E_abs_max(1,j) = max(E_abs(j,:));
         Edb(j,:) = db(E_abs(j,:)) - db(E_abs_max(1,j));
         top1 = findpeaks(Edb(j,:));
         top = [top,100,top1];
     end
     top = [top,100];
     location = find(top == 100);

     %计算后代相对应的误差值
     error = [];
     for i = 1:m
         infer = location(i);
         later = location(i+1);
         E1 = top(infer+1:later-1);
         error1 = sum((E1+25).*(E1+25))- 25.*25;
         error = [error,error1];
     end
     
     
     
     %将交叉变异之后的个体保留,并将差值更新
     for i = 1:m
         if error_inital(hang(i)) > error(i)
             error_inital(hang(i)) = error(i);
             for num = 1:NL*L
                 pop(hang(i),num) = I_select(i,num);
             end
         end
     end
     
     
     %将之前适应度最好的电流值对应的二进制数代替当前适应度最差的
     Max = max(error_inital);
     rr_Max = find(error_inital == Max);
     for i = NL*L
         pop(rr_Max(1),i) = Best(i);
     end
     min(error_inital)
end

%%
%将对应的电流二进制数变换成十进制
for j = 1:NP
    a = 1;
    for i = 1:L:NL*L  
        I_act = sum(Decode.*pop(j,i:(i+L-1)));%对应种群中的第j个个体的第n个电流量
        I(j,a) = I_act;
        a = a + 1;
    end
end
%计算相对应的主瓣和旁瓣的电平值
for j = 1:NP
    for num = 1:NS
        E(j,num)=sum(I(j,:).*exp(1j*(wavenum*X*sind(theta(num))+alpha)))+1e-3;
    end
end
E_abs = abs(E);
top = [];hang = [];top1 = [];
for j = 1:NP
    E_abs_max(1,j) = max(E_abs(j,:));
    Edb(j,:) = db(E_abs(j,:)) - db(E_abs_max(1,j));
    top1 = findpeaks(Edb(j,:));
    top = [top,100,top1];
end
min_error = min(error);
rr = find(error == min_error);
[m,n] = size(Edb);
Edb1 = [];
for i = 1:n
    Edb1= [Edb1,Edb(rr(1),i)];
end
I(rr,:)
%%
figure()
plot(theta,Edb1,'LineWidth',2);%标准化
xlim([-90 90])
xlabel('单元数N是20,半波长')
grid on
 

3、切比雪夫多项式系数代码

I = zeros(21,21);   %切比雪夫多项式系数存放矩阵

I(1,1) = 1; I(2,1) = 0; I(2,2) = 1;  %设定切比雪夫初值

%% 根据切比雪夫迭代公式进行多项式系数求解T（n+1）= 2x*T（n）-T（n-1）

for i = 3:21

    I(i,1) = -I(i-2,1);

    for j = 2:21

        I(i,j) = 2*I(i-1,j-1)-I(i-2,j);

    end

end

for i = 13:21

    i

    I(i,:)

end

4、切比雪夫多项式均匀直线阵

%% --------------------------------------------------------------------------
% 切比雪夫低副瓣阵列综合
% 设计一个间距为d,单元数为N，主副瓣电平比为RdB，扫描角度为theta0的切比雪夫阵列。
%
%--------------------------------------------------------------------------
%% 初始数据赋值
clear
clc
N = 20;                                  %单元数N（3<N<=13,N取整数）
if rem(N,2)==0                          %求和项数M(奇偶不同)
    M = N/2;
else
    M = (N-1)/2+1;
end
RdB = 25;                               % 主副瓣比(dB值)
lamuda = 10;                            % 波长
d = 0.5*lamuda;                         % 单元间距
theta0 = 90/180*pi;                     % 扫描角度，相对于阵列排布方向的夹角
A = [1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;       % chebyshev多项式Tn(x) = cos(nu)= f(x)系数矩阵A
    0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;        % 系数矩阵A每一行表示n，从n = 0开始
    -1,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;       % 列表示x的幂次方，从0次方开始
    0,-3,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
    1,0,-8,0,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
    0,5,0,-20,0,16,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
    -1,0,18,0,-48,0,32,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
    0,-7,0,56,0,-112,0,64,0,0,0,0,0,0,0,0,0,0,0,0,0;
    1,0,-32,0,160,0,-256,0,128,0,0,0,0,0,0,0,0,0,0,0,0;
    0,9,0,-120,0,432,0,-576,0,256,0,0,0,0,0,0,0,0,0,0,0;
    -1,0,50,0,-400,0,1120,0,-1280,0,512,0,0,0,0,0,0,0,0,0,0;
    0,-11,0,220,0,-1232,0,2816,0,-2816,0,1024,0,0,0,0,0,0,0,0,0;
    1,0,-72,0,840,0,-3584,0,6912,0,-6144,0,2048,0,0,0,0,0,0,0,0;
    0,13,0,-364,0,2912,0,-9984,0,16640,0,-13312,0,4096,0,0,0,0,0,0,0;
    -1,0,98,0,-1568,0,9408,0,-26880,0,39424,0,-28672,0,8192,0,0,0,0,0,0;
    0,-15,0,560,0,-6048,0,28800,0,-70400,0,92160,0,-61440,0,16384,0,0,0,0,0;
    1,0,-128,0,2688,0,-21504,0,84480,0,-180224,0,212992,0,-131072,0,32768,0,0,0,0;
    0,17,0,-816,0,11424,0,-71808,0,239360,0,-452608,0,487424,0,-278528,0,65536,0,0,0;
    -1,0,162,0,-4320,0,44352,0,-228096,0,658944,0,-1118208,0,1105920,0,-589824,0,131072,0,0;
    0,-19,0,1140,0,-20064,0,160512,0,-695552,0,1770496,0,-2723840,0,2490368,0,-1245184,0,262144,0;
    1,0,-200,0,6600,0,-84480,0,549120,0,-2050048,0,4659200,0,-6553600,0,5570560,0,-2621440,0,524288];
% 初始矩阵赋值
I = zeros(1,M);                         % 电流幅度矩阵
S = zeros(M,M);                         % 阵因子系数矩阵
S_compare = zeros(1,M);                 % 系数比对矩阵
R = 10^(RdB/20);                        % 非dB 值的主副瓣比
x0 = 1/2*( (R+sqrt(R^2-1))^(1/(N-1))+...% 变量代换值x0
            (R-sqrt(R^2-1))^(1/(N-1))  );
%% 求S、S_compare和I
% 从系数矩阵中择选出M个求和项对应的系数S(奇偶分开讨论)
for i = 1:M
    if rem(N,2)==0                      % 偶数情况
        for j = 1:M                     % 第i行表示x的i次方，
            S(i,j) = A(2*j,2*i);        % 第j列表示第j个求和项系数（未除x0)
        end
        S_compare(i) = A(N,2*i);        % 比对矩阵，即下标为N-1的chebyshev多项式的系数
    else                                % 奇数情况
        for j = 1:M                    
            S(i,j) = A(2*j-1,2*i-1);   
        end
        S_compare(i) = A(N,2*i-1);    
    end
end
% 通过S和S_compare系数比对求出电流幅度
for k = 1:M
    i = M-k+1;
    if rem(N,2)==0                      % 偶数
        I(i) = (S_compare(i)*x0^(2*i-1) -...
            I*S(i,:)')/S(i,i);
    else                                % 奇数
        I(i) = (S_compare(i)*x0^(2*(i-1)) -...
            I*S(i,:)')/S(i,i);
    end
end
I = I/max(I);                         % 对I归一化
if rem(N,2)==0
    I_final = [fliplr(I),I];          % 最终的单元排列（左右对称）
else
    I_final = [fliplr(I),I(2:end)];
end
sprintf('天线单元归一化电流幅度：')
sprintf('%.3f  ',I_final)
%% 获得最终阵列方向图S_P
theta_rad = 0:0.01:pi;
theta = theta_rad*180/pi;
u = pi*d/lamuda*( cos(theta_rad)- cos(theta0));
S_P = zeros(1,length(theta_rad));       % 最终方向图
for k = 1:M
    if rem(N,2)==0
        S_P = S_P + I(k)*cos((2*k-1)*u);% 偶数
    else
        S_P = S_P + I(k)*cos(2*(k-1)*u);% 奇数
    end
end
S_P_abs = abs(S_P);                     % 对S_P取绝对值
S_PdB = 20*log10(S_P_abs/max(S_P_abs)); % 对S_P取dB值
%% 绘图
H = -ones(1,length(S_P_abs))*25;        % 根据预先设置的主副瓣比得到的参考曲线
% 直角坐标系
figure('NumberTitle', 'off', 'Name', 'S Parameter (abs)-Plot');
plot(theta,S_P_abs,'b','LineWidth',1.5)
xlabel('theta(°)')
ylabel('|S| ')
title('chebyshev低副瓣阵列直角坐标图')
figure('NumberTitle', 'off', 'Name', 'S Parameter (dB)-Plot');
plot(theta,H,'r--','LineWidth',1.5)
hold on
plot(theta,S_PdB,'b','LineWidth',1.5)
xlabel('theta(°)')
ylabel('|S| dB')
title('chebyshev低副瓣阵列直角坐标图')
legend('预设副瓣参考曲线','方向图')
% 极坐标系
figure('NumberTitle', 'off', 'Name', 'S Parameter (dB)-Polar');
polarplot(theta_rad,H,'r--','LineWidth',1.5)
hold on
polarplot(theta_rad,S_PdB,'b','LineWidth',1.5)
thetalim([0 180]);
rmin = S_PdB(1,1);
rmax = max(S_PdB);
rlim([-50 rmax]);
title('chebyshev低副瓣阵列极坐标图')
legend('预设副瓣参考曲线RdB','方向图(dB)')