na_gs.m

通信系统的模拟仿真程序

%绘制负倒描述函数图、乃氏图


disp('线性环节传函G(s)为：')
printsys(num,den,'s');

ngwin=figure('Position',[120 49 513 491],'Name','描述函数法','NumberTitle','off',...
      'MenuBar','none');

axes('position',[0 0 1 1]);
axis off;
dd=text(0.1,0.18,'实轴范围:');
set(dd,'fontsize',15,'color',[1 1 1],'fontname','宋体');
dd=text(0.1,0.08,'虚轴范围:');
set(dd,'fontsize',15,'color',[1 1 1],'fontname','宋体');

%设置绘图区域
axes('Position',[.15 .33 .8 .56]);

%绘制乃氏图，得到乃氏图数据矩阵gs_nqst、w矩阵w
[re_g,im_g,w]=nyquist(num,den);
gs_nqst=re_g+sqrt(-1)*im_g;
nyquist(num,den);
ylabel('虚     轴');
xlabel('实     轴');
title('描述函数法');
grid;

%得到当前坐标轴范围
axis0=axis;
xaxis=[axis0(1),axis0(2)];  %横坐标轴范围
yaxis=[axis0(3),axis0(4)];  %纵坐标轴范围


%绘制负倒描述函数图
disp(' ');
disp('非线性环节为:');

%幅值A矩阵为列矩阵，范围为(a,100*a),间隔为0.005
A=n_para(1)+0.0001:0.005:n_para(1)*100;
A=A';
len=length(A);   % A的长度

%  求描述函数矩阵N_A
if i==1
  %死区的描述函数，a为开关值，k为斜率
  a=n_para(1);
  k=n_para(2);
  disp(['  死区环节,其中a: ',num2str(a),' ; k: ',num2str(k),' .']);
  
  N_A=2*k/pi*(pi/2-asin(a*ones(len,1)./A)...
      -(a*ones(len,1)./A).*sqrt(1-(a*ones(len,1)./A).^2));

elseif i==2
  %饱和的描述函数，a为开关值，k为斜率
  a=n_para(1);
  k=n_para(2);
  disp(['  饱和环节,其中a: ',num2str(a),' ; k: ',num2str(k),' .']);
  
  N_A=2*k/pi*(asin(a*ones(len,1)./A)...
      +(a*ones(len,1)./A).*sqrt(1-(a*ones(len,1)./A).^2));

elseif i==3
  %间隙的描述函数，b为开关值，k为斜率
  b=n_para(1);
  k=n_para(2);
  disp(['  间隙环节,其中b: ',num2str(b),' ; k: ',num2str(k),' .']);

  N_A=k/pi*(pi/2+asin(1-2*b*ones(len,1)./A)...
            +2*(1-2*b*ones(len,1)./A).*sqrt((b*ones(len,1)./A).*(1-b*ones(len,1)./A)))...
      +sqrt(-1)*4*k*b*(ones(len,1)./(pi*A)).*(b*ones(len,1)./A-1);

elseif i==4
  %继电器的描述函数，h为上升值，m为继电器系数，M为开时的输出值
  h=n_para(1);
  m=n_para(2);
  M=n_para(3);
  disp(['  继电器环节,其中h: ',num2str(h),' ; m: ',num2str(m),' ; M: ',num2str(M),' .']);

  N_A=2*M*(ones(len,1)./(pi*A)).*(sqrt(1-(m*h)^2*ones(len,1)./A.^2)+sqrt(1-h^2*ones(len,1)./A.^2))...
      +sqrt(-1)*2*M*h*ones(len,1)./(pi*A.^2)*(m-1);

end
disp(' ');

% 求负倒描述函数矩阵f_na
f_na=-ones(len,1)./N_A;

%判断f_na的单调性
differ=diff(f_na);
acsend=differ>0;

 
if nonzeros(imag(f_na))==[]  % f_na为实数阵
   hold on;
   plot(f_na(1),0,'yo');     %在起始点出画一'o';表示不包括该点
   hold on;
   if acsend==1              %当f_na单调递增时，f_na曲线为红色
      plot(f_na,zeros(len,1),'r');
   elseif acsend==0          %当f_na单调递减时，f_na曲线为绿色
      plot(f_na,zeros(len,1),'g');
   else                      %当f_na不单调时，f_na曲线为黄色
      plot(f_na,zeros(len,1),'y');
   end;  
else                         % f_na为虚数阵                    
   hold on;
   plot(f_na(1),'yo');       %在起始点出画一'o';表示不包括该点
   hold on;
   if acsend==1              %当f_na单调递增时，f_na曲线为红色
      plot(f_na,'r');
   elseif acsend==0          %当f_na单调递减时，f_na曲线为绿色
      plot(f_na,'g');        
   else                      %当f_na不单调时，f_na曲线为黄色
      plot(f_na,'y');
   end;
end;

% 改变坐标轴范围时，判断输入是否正确，如正确，设置坐标轴范围；
%  否则，弹出错误信息，恢复原坐标轴范围
axis_com=['xaxis=eval(get(x_edit,''string''));yaxis=eval(get(y_edit,''string''));',...
          'if (xaxis(1)<xaxis(2))&(yaxis(1)<yaxis(2))',...
          '   axis([xaxis yaxis]);',...
          'else  ',...  
          '   errordlg(''坐标轴输入错误!'',''错误信息'',''on'');',...
          '   axis0=axis;',...
          '   xaxis=[axis0(1),axis0(2)];',... 
          '   yaxis=[axis0(3),axis0(4)];',... 
          '   set(x_edit,''string'',[''['',num2str(xaxis(1)),'' '',num2str(xaxis(2)),'']'']);',...
          '   set(y_edit,''string'',[''['',num2str(yaxis(1)),'' '',num2str(yaxis(2)),'']'']);',...          
          'end;'];

% 横轴范围输入框
x_edit=uicontrol(ngwin,'style','edit','units','normalized',...
       'position',[0.3,0.15,0.3,0.05],'string',['[',num2str(xaxis(1)),' ',num2str(xaxis(2)),']'],'back',[1 1 1],...
       'call','axis_com');
% 纵轴范围输入框
y_edit=uicontrol(ngwin,'style','edit','units','normalized',...
       'position',[0.3,0.05,0.3,0.05],'string',['[',num2str(yaxis(1)),' ',num2str(yaxis(2)),']'],'back',[1 1 1],...
       'call','axis_com');
% 查找输入点对应的A,w 按钮
uicontrol('style','push','units','normal','pos',[0.7 0.15 0.13 0.05],'string','find',...
         'call','find_aw')
% 关闭按钮
uicontrol('style','push','units','normal','pos',[0.7 0.05 0.13 0.05],'string','&Close',...
         'call','close')