psrdf.m

伪速率脉冲调制器的稳定性分析

function psrdf
echo off,  clc
	He=0.03;	Ha=0.005;	Km=1;	Tm=2;	To=Tm; 	C=1; 	Nflag=1;	Dtmin=0.05;
  	He=0.05;	Ha=0.02;	Km=1;	Tm=2;	To=Tm; 	C=1; 	Nflag=1;	Dtmin=0.05;
  	He=0.04;	Ha=0.03;	Km=1;	Tm=2.0;	To=Tm; 	C=1; 	Nflag=1;	Dtmin=0.006;

   Wn=[0.05,0.2,0.8,3,10,20];
   Wn=[9.255,12.71,12.21,1,1.074,0.8505,3.251,9.262,1,1.073];
   Wn=[0.001 0.033 0.1 0.178 0.3086 0.5694 0.9815 1.8477 2 3]*2*pi;
   
   % some variance is define
   h=He-Ha;nw=length(Wn);
   
for vww=1:nw,
  w=Wn(vww);
     disp([' W= ',num2str(w),'   .... Working ....  Waiting.........'])
  a=1.02*He;  p2=1;  m1=1;
while 1,
  if a<=1.02*He,    [p1,p2,rea,ima,ys,ya,yc,a]=soldfjx(w,He,Ha,Tm,C,a,To,Km,Dtmin,m1);
  else,    
      [p1,p2,rea,ima,ys,ya,yc]=soldffz(w,He,Ha,Tm,C,a,To,Km,1,Dtmin);
  end
    if p1==0,break,end
       naa(m1,vww)=a;  ree(m1,vww)=rea; imm(m1,vww)=ima; yys(m1,vww)=ys;
       gnr(m1,vww)=ya; gncr(m1,vww)=yc; nnaa(vww)=m1;
    if m1>=20, break,end
    a=a+He*m1^2*sqrt(w)/200; m1=m1+1;
  clear  rea ima  ys  a00
end
end
[mm,nn]=size(naa);
    for kk1=1:nn,
	for kk2=1:mm,
	    if naa(kk2,kk1)==0,
	       naa(kk2,kk1)=naa(kk2-1,kk1);   ree(kk2,kk1)=ree(kk2-1,kk1);
	       imm(kk2,kk1)=imm(kk2-1,kk1);   yys(kk2,kk1)=yys(kk2-1,kk1);
	       gnr(kk2,kk1)=gnr(kk2-1,kk1);   gncr(kk2,kk1)=gncr(kk2-1,kk1);
	     end
	 end
     end
      naw=ree+sqrt(-1)*imm;
      rou=abs(naw); ndbc=-20*log10(rou); ndba=-20*log10(yys);
      nph=-angle(naw)*180/pi-180;ndbr=-20*log10(gnr);ndbcr=-20*log10(gncr);
%      
NM=2;
save psrdf Wn naa ndbc ndba ndbr ndbcr nph naw rou
%plot()
figure(1);
line1=line(nph(1,:),ndbc(1,:),'color',[1 0 0]);
line2=line(nph,ndbc,'color',[1 0 0]);
xlabel('Phase/Deg'),ylabel('Gain/dB')
figure(2);
line3=line(nph(1,:),ndba(1,:),'color',[0 0 1]);
line4=line(nph,ndba,'color',[0 0 1]);
xlabel('Phase/Deg'),ylabel('Gain/dB')
figure(3);
line5=line(nph(1,:),ndbr(1,:),'color',[0 1 1]);
line6=line(nph,ndbr,'color',[0 1 1]);
xlabel('Phase/Deg'),ylabel('Gain/dB')
figure(4);
line7=line(nph(1,:),ndbcr(1,:),'color',[0 0 0]);
line8=line(nph,ndbcr,'color',[0 0 0]);
xlabel('Phase/Deg'),ylabel('Gain/dB')


%	record('','Nonlinear NIDS','record.rep')
%	record(ndbc,'gain(dB)','record.rep')
%	record(nph,'phase(deg)','record.rep')
%	record(Wn,'Frequency','record.rep')

function [p1,p2,re,im,gna,gnr,gncr,a]=soldfjx(w,he,ha,tff,c,a,to,kf,dtmin,m1)
% SOLDFJX
%	 Nonlinear function equation method solve pulse parameter
    save ww
    h=he-ha;p2=1;  format short,format compact
    OPTIONS=fsolve('defaults');
    OPTIONS = OPTIMSET(OPTIONS,'TolFun',1.e-8,'MaxIter',800,'LargeScale','off');
    xx0=[pi/w/2.2,1.1*dtmin,-h/2,h/2,-h/2,he]';   % you may modify it
	if m1==1,  xy=fsolve('gtest28',xx0,OPTIONS);gtest28(xy);
	else,xx0(6)=[];xx0(5)=[]; xy=fsolve('gtest24',xx0,OPTIONS);gtest24(xy);
	end
   if length(xy)~=0,
	tpo=xy(1);tpc=xy(2)+tpo;p1=1;if m1==1,a=xy(6);end
	  disp('     w        a          dt         tonw     toffw')
	  disp([w a tpc-tpo tpo*w tpc*w])
	if xy(2)<=0,p1=0;return,end
   else,  p1=0;return
   end
 if xy(2)<dtmin ,
	 clear xy tpo tpc p1
     xx0=[pi/w/3,-h/2,h/2,-h/2,he]';          % yuo may modify it
	   if m1==1,xy=fsolve('gtest29',xx0,OPTIONS);
	   else,xx0(5)=[];xx0(4)=[]; xy=fsolve('gtest26',xx0,OPTIONS); end
	if  length(xy)~=0,
	   p1=1;tpo=xy(1);tpc=tpo+dtmin;if m1==1,a=xy(5);end
	   disp('     w        a          dt        tonw     toffw')
	   disp([w a tpc-tpo tpo*w tpc*w])
	 else,  p2=0;p1=0;return
	 end
 end
    dtw=(tpc-tpo)*w;    re=2*(cos(w*tpo)-cos(w*tpc))/pi/a;
    im=2*(sin(w*tpc)-sin(w*tpo))/pi/a;    gna=dtw/2/a;
    gnr=sqrt(2*dtw/pi)/a;    gncr=sqrt((sin(dtw/2))^2+(sin(3*dtw/2))^2/9)/a*4/pi;
    delete ww.mat
    return
    
    
 function [p1,p2,rea,ima,ys,gnr,gncr]=soldffz(w,he,ha,tff,c,a,to,kf,p3,dtmin)
% SOLDFFZ
%	    get psrm pulse parameters by simulation

   tt=2*pi./w;t2=tt;h=he-ha;af(1)=-h/2;
   aff=af(1); eps1=1; y=0;   pp=0; po=-dtmin;
   while abs(eps1)>he/1000,
      pp=pp+1;  if pp>20,break ,end
      t1=t2-tt;t2=t1;  clear n1 n2 tpc tpo yy1 yy2
      n1=0;n2=0;
      while t2 < tt,		 dt=tt/200;
         while t2 < tt,       t2=t2+dt;	    if t2>tt,t2=tt;end
            if y==0,af(2)=af(1)*exp(-(t2-t1)/tff);
            else, af(2)=af(1)*exp(-(t2-t1)/to)+kf*y*(1-exp(-(t2-t1)/to));
            end
            e=a*sin(w*t2)-af(2);
            if abs(e)>=he & y==0,
               if dt<h*to/kf/200 | dt<dtmin/200,
                  n1=n1+1;tpo(n1)=t2; y=c*sign(e); yy1(n1)=y;po=t2;  break
               end
               af(1)=af(1)*exp(-(t2-t1-dt)/tff);  t1=t2-dt; t2=t1; dt=dt/2;
            elseif (y==c & e<=ha) | (y==-c & e>=-ha),
               ss=t2-po;if ss<0,ss=ss+tt;end
               if ss>dtmin,
                  if dt<h*to/kf/500 | dt<dtmin/500,
                     n2=n2+1; tpc(n2)=t2; yy2(n2)=y;  y=0; break
                  end
               elseif ss<=dtmin,     t2=po+dtmin;if t2>tt,t2=t2-tt;end
                  n2=n2+1; tpc(n2)=t2; yy2(n2)=y;  y=0; break
               end
               af(1)=af(1)*exp(-(t2-t1-dt)/to)+kf*y*(1-exp(-(t2-t1-dt)/to));
               t1=t2-dt; t2=t1; dt=dt/2;
            end
         end
         t1=t2; af(1)=af(2);
      end
      eps1=-aff+af(2);  af(1)=af(2); aff=af(1);
   end
   
   % find pulse open and close time tpo tpc
   p1=1;p2=1;
   if n1~=0,
      n=length(yy1); im=0; re=0 ; ys=0; gnr=0;gncr=0;wdt=0;
      if n>8,p1=0;rea=nan;ima=nan;ys=nan;gnr=nan;gncr=nan;return,end  % more than 2 pulse
      if tpo(n1)>tpc(n2), tpc(n2+1)=tt; yy1(n2+1)=yy1(n1); end
      if tpo(1)>tpc(1),
         for k=n1:-1:1, tpo(k+1)=tpo(k); yy1(k)=yy2(k); end
         tpo(1)=0;
      end
      for k=1:n,
         re=re+yy1(k)*(cos(w*tpo(k))-cos(w*tpc(k)));
         im=im+yy1(k)*(sin(w*tpc(k))-sin(w*tpo(k)));
         wdt=wdt+w*(tpc(k)-tpo(k));
      end
      rea=re/a/pi; ima=im/a/pi;
      wdt=wdt/2; ys=wdt/2; ys=ys/a;
      gnr=sqrt(wdt*2/pi);gnr=gnr/a;
      gncr=sqrt((sin(wdt/2))^2+(sin(3*wdt/2))^2/9);
      gncr=gncr*4/pi/a;
      clear tpc tpo n1 n2 yy1 yy2
   end
   disp([p1,p2,rea,ima,ys,gnr,gncr])
   return
      
    
function yy=gtest24(xx)
%GTEST24
%	Nonlinear equations about genal- psrm
%	xx(1) : t1
%	xx(2) : dt
%	xx(3) : af(t1)
%	xx(4) : af(t1+dt)

  load ww
  yy=zeros(size(xx));
  yy(1)=a*sin(w*xx(1))-xx(3)-he;
  yy(2)=a*sin(w*(xx(1)+xx(2)))-xx(4)-ha;
  yy(3)=xx(4)-xx(3)*exp(-xx(2)/to)-kf*(1-exp(-xx(2)/to));
  yy(4)=xx(3)+xx(4)*exp((-pi/w+xx(2))/tff);
  return
  
function yy=gtest26(xx)
%GTEST26
%	Nonlinear equations about genal- psrm
%	     dt=dtmin
%	xx(1) : t1
%	xx(2) : af(t1)
%	xx(3) : af(t1+dt)

  load ww
  yy=zeros(size(xx));
  yy(1)=a*sin(w*xx(1))-xx(2)-he;
  yy(2)=xx(3)-xx(2)*exp(-dtmin/to)-kf*(1-exp(-dtmin/to));
  yy(3)=xx(2)+xx(3)*exp((-pi/w+dtmin)/tff);
  return
  
  function yy=gtest28(xx)
% GTEST28
%	 nonlinear equations about minimum amplititude input
%	 xx(1) : t1
%	 xx(2) : dt
%	 xx(3) : af(t1)
%	 xx(4) : af(t1+dt)
%	 xx(5) : af(0)
%	 xx(6) : amin

  load ww
  yy=zeros(size(xx));
  yy(1)=xx(6)*sin(w*xx(1))-xx(3)-he;
  yy(2)=xx(6)*sin(w*(xx(1)+xx(2)))-xx(4)-ha;
  yy(3)=xx(4)-xx(3)*exp(-xx(2)/to)-kf*c*(1-exp(-xx(2)/to));
  yy(4)=xx(3)+xx(4)*exp((-pi/w+xx(2))/tff);
  yy(5)=xx(3)-xx(5)*exp(-xx(1)/tff);
  yy(6)=xx(6)*w*cos(w*xx(1))+xx(3)/tff;  
  return
  
  function yy=gtest29(xx)
% GTEST29
%	 nonlinear equations about minimum amplititude input
%	     dt=dtmin
%	 xx(1) : t1
%	 xx(2) : af(t1)
%	 xx(3) : af(t1+dt)
%	 xx(4) : af(0)
%	 xx(5) : amin

  load ww
  yy=zeros(size(xx));
  yy(1)=xx(5)*sin(w*xx(1))-xx(2)-he;
  yy(2)=xx(3)-xx(2)*exp(-dtmin/to)-kf*c*(1-exp(-dtmin/to));
  yy(3)=xx(2)+xx(3)*exp((-pi/w+dtmin)/tff);
  yy(4)=xx(2)-xx(4)*exp(-xx(1)/tff);
  yy(5)=xx(5)*w*cos(w*xx(1))+xx(2)/tff;
  return