find_aw.m

xnsy算法的matlab实现

%find a w
%

%变量information记录所选点在描述函数图上的信息
information=['                                                 '];
% x,y为鼠标输入点的坐标
[x,y]=ginput(1);

%变量d2gs为输入点到乃氏图的距离；
%变量d2gsmin点到乃氏图的最小距离；gs_i为最小距离点的下标
d2gs=abs(gs_nqst-(x+sqrt(-1)*y));
[d2gsmin,gs_i]=min(d2gs);

%求取精确的w值
if d2gsmin<0.4      
  % errw为w的最大误差
  errw=max(abs(w(gs_i)-w(gs_i-1)),abs(w(gs_i)-w(gs_i+1)));
  
  gs_nqst1=gs_nqst;
  gs_i1=gs_i;
  w1=w;
  %如果errw大于0.01时，将w1细分以后继续求乃氏图矩阵gs_nqst1
  while errw>=0.01
     w1=linspace(w1(gs_i1-1),w1(gs_i1+1),50);
     [re1,im1,w1]=nyquist(num,den,w1);   
     gs_nqst1=re1+sqrt(-1)*im1;
     d2gs1=abs(gs_nqst1-(x+sqrt(-1)*y));
     [d2gsmin1,gs_i1]=min(d2gs1);
     if gs_i1==length(w1)
        errw=abs(w1(gs_i1)-w1(gs_i1-1));
     else 
        errw=max(abs(w1(gs_i1)-w1(gs_i1-1)),abs(w1(gs_i1)-w1(gs_i1+1)));
     end;  
  end;
  %如果点到乃氏图的最小距离小于0.01时，可认为所选点在乃氏图上  
  if d2gsmin1<0.01
     information=['在乃氏图上 ',information];
  end;

end;

%求取A值
%如果f_na在实轴上且不单调时，需分为两个单调矩阵来分别求A值
[max_v,max_i]=max(f_na);
if (imag(f_na)==0)&(max_i~=1)&(max_i~=len)
   f_na1=f_na(1:max_i);
   A1=A(1:max_i);
   f_na2=f_na(max_i+1:len);
   A2=A(max_i+1:len);
%如果f_na不在实轴上或单调时，只需用一个单调矩阵来求A值
% 将第二个矩阵为空
else     
   f_na1=f_na;
   A1=A;
   f_na2=[];
   A2=[];
end;
   
%变量d2f_na1到f_na曲线的距离
d2f_na1=abs(f_na1-(x+sqrt(-1)*y));
[d2f_na1min,na1_i]=min(d2f_na1);

%如果点到f_na1图的最小距离小于0.02时，可认为所选点在f_na1图上  
if d2f_na1min<=0.02
  errA1=max(abs(A1(na1_i)-A1(na1_i-1)),abs(A1(na1_i)-A1(na1_i+1)));
  information=['在-1/N(A)图1上 ',information];
end;

%如果f_na在实轴上且不单调时，需对f_na2矩阵再次求A值
if (f_na2~=[])&(A2~=[])

   d2f_na2=abs(f_na2-(x+sqrt(-1)*y));
   [d2f_na2min,na2_i]=min(d2f_na2);

   if d2f_na2min<=0.02
      errA2=max(abs(A2(na2_i)-A2(na2_i-1)),abs(A2(na2_i)-A2(na2_i+1)));
      information=['在-1/N(A)图2上 ',information];
   end;
end;

disp(' ');
disp(['所选点为: ',num2str(x+sqrt(-1)*y)]);
hold on;
plot(x,y,'w+');

if findstr(information,'在-1/N(A)图2上 在-1/N(A)图1上 在乃氏图上 ')~=[]
   disp('该点为乃氏图与负倒描述函数图的交点!且是两个交点.');
   disp(['   其中w为:  ',num2str(w1(gs_i1)),' .']);
   disp(['       A1为:  ',num2str(A1(na1_i)),' .']);
   disp(['       A2为:  ',num2str(A2(na2_i)),' .']);
elseif findstr(information,'在-1/N(A)图')~=[]&findstr(information,'在乃氏图上 ')~=[]
   if findstr(information,'在-1/N(A)图1上 ')~=[]
      disp('该点为乃氏图与负倒描述函数图的交点!');
      disp(['   其中w为:  ',num2str(w1(gs_i1)),' .']);
      disp(['       A为:  ',num2str(A1(na1_i)),' .']);
   else
      disp('该点为乃氏图与负倒描述函数图的交点!');
      disp(['   其中w为:  ',num2str(w1(gs_i1)),' .']);
      disp(['       A为:  ',num2str(A2(na2_i)),' .']);
   end;
elseif findstr(information,'在乃氏图上 ')~=[]
   disp('该点在乃氏图上!');
   disp(['   其中w为:  ',num2str(w1(gs_i1)),' .']);
elseif findstr(information,'在-1/N(A)图')~=[]
   if findstr(information,'在-1/N(A)图2上 在-1/N(A)图1上 ')~=[]
      disp('该点为负倒描述函数图的点!且是两个点.');
      disp(['   其中A1为:  ',num2str(A1(na1_i)),' .']);
      disp(['       A2为:  ',num2str(A2(na2_i)),' .']);
   elseif findstr(information,'在-1/N(A)图1上 ')~=[]
      disp('该点为负倒描述函数图的点!');
      disp(['   其中A为:  ',num2str(A1(na1_i)),' .']);
   else
      disp('该点为负倒描述函数图的点!');
      disp(['   其中A为:  ',num2str(A2(na2_i)),' .']);
   end;
else
   disp('该点既不在乃氏图上,也不在负倒描述函数图上!');
end;