sa.m

这是一个用于语音信号处理的工具箱

% sa.m
% Simulated Annealing Algorithm

	ii1 = rem(floor(iseed1/177),177) + 2;
	j = rem(iseed1,177) + 2;
	k = rem(floor(iseed2/169),178) + 1;
	l = rem(iseed2,169);
	for ii=1:97,
		s = 0.0;
		t = 0.5;
		for jj=1:24,
			m = rem(rem(ii1*j,179)*k,179);
			ii1 = j;
			j = k;
			k = m;
			l = rem(53*l+1,169);
			if (rem(l*m,64) >= 32)
				s = s + t;
			end
			t = 0.5*t;
		end
		u(ii) = s;
	end
	rc = 362436.0 / 16777216.0;
	cd = 7654321.0 / 16777216.0;
	cm = 16777213.0 /16777216.0;
	I97 = 97;
	J97 = 33;


	NACC = 0;
	NOBDS = 0;
	Nfcnev = 0;
	for ii=1:Npara,
		xopt(ii) = xvect(ii);
		Nacp(ii) = 0;
	end
	for ii=1:NEPS,
		fstar(ii) = 1.0e20;
	end
	[f,area,leng]=fcn(Npara,xvect,currf,form_tar,area,leng);
	f=-f;
	Nfcnev = Nfcnev + 1;
	fopt = f;
	fstar(1) = f;
%       if (iprint >= 1)
%               prt2(f);
	quit=0;
while (Nfcnev < maxevl & quit ==0)

	nup = 0;
	nrej = 0;
	nnew = 0;
	ndown = 0;
	lnobds = 0;
	for m=1:NT,
	   for j=1:NS,
		for h=1:Npara,
		    for ii=1:Npara,
			if (ii==h)
%%% ranmar_start
				uni = u(I97) - u(J97);
				if (uni < 0.0)
					uni = uni + 1.0;
				end
				u(I97) = uni;
				I97 = I97 - 1;
				if (I97 == 0)
					I97 = 97;
				end
				J97 = J97 - 1;
				if (J97 == 0)
					J97 = 97;
				end
				rc = rc - cd;
				if (rc < 0.0)
					rc = rc + cm;
				end
				uni = uni - rc;
				if (uni < 0.0)
					uni = uni + 1.0;
				end
%%% ranmar_stop
			    xp_sa(ii)=xvect(ii)+(uni*2.0-1.0)*vm(ii);
			else
			    xp_sa(ii)=xvect(ii);
			end
			if (xp_sa(ii)<lb(ii)) | (xp_sa(ii)>ub(ii)) 
%%% ranmar_start
				uni = u(I97) - u(J97);
				if (uni < 0.0)
					uni = uni + 1.0;
				end
				u(I97) = uni;
				I97 = I97 - 1;
				if (I97 == 0)
					I97 = 97;
				end
				J97 = J97 - 1;
				if (J97 == 0)
					J97 = 97;
				end
				rc = rc - cd;
				if (rc < 0.0)
					rc = rc + cm;
				end
				uni = uni - rc;
				if (uni < 0.0)
					uni = uni + 1.0;
				end
%%% ranmar_stop
			    xp_sa(ii) = lb(ii) + (ub(ii) - lb(ii))*uni;
			    lnobds = lnobds + 1;
			    NOBDS = NOBDS + 1;
			end
		    end
		    [fp,area,leng] = fcn(Npara,xp_sa,currf,form_tar,area,leng);
		    fp=-fp;
		    Nfcnev = Nfcnev + 1;
          if (Nfcnev >= maxevl)
			fopt = -fopt;
			set(m4_1_h,'Value',xp_sa(1));
			set(m4_2_h,'String',round(xp_sa(1)));
			set(m4_3_h,'Value',xp_sa(2));
			set(m4_4_h,'String',round(xp_sa(2)));
			set(m4_5_h,'Value',xp_sa(3));
			set(m4_6_h,'String',round(xp_sa(3)));
			set(m4_7_h,'Value',xp_sa(4));
			set(m4_8_h,'String',round(xp_sa(4)));
			set(m4_9_h,'Value',xp_sa(5));
			set(m4_10_h,'String',round(xp_sa(5)));
			set(m4_11_h,'Value',xp_sa(6));
			set(m4_12_h,'String',round(xp_sa(6)));
			set(m4_13_h,'Value',xp_sa(7));
			set(m4_14_h,'String',round(xp_sa(7)));
			if (Npara==9 | Npara==12)
				set(m4_15_h,'Value',xp_sa(9));
				set(m4_16_h,'String',round(xp_sa(9)));
			else
				set(m4_15_h,'Value',0);
				set(m4_16_h,'String',0);
			end
			set(m4_17_h,'Value',xp_sa(8));
			set(m4_18_h,'String',round(xp_sa(8)));
			art_set(currf,1)=xp_sa(1);
			art_set(currf,2)=xp_sa(2);
			art_set(currf,3)=xp_sa(3);
			art_set(currf,4)=xp_sa(4);
			art_set(currf,5)=xp_sa(5);
			art_set(currf,6)=xp_sa(6);
			art_set(currf,7)=xp_sa(7);
			art_set(currf,9)=xp_sa(8);
			if Npara==8 
				art_set(currf,8)=artvelum;
			 	wh=1.34651;
				hk1=1.3431;
			        g1k=0.82081;
			elseif Npara==9
				art_set(currf,8) = xp_sa(9);
			 	wh=1.34651;
			        hk1=1.3431;
				g1k=0.82081;
			elseif Npara==11
				art_set(currf,8)=artvelum;
				art_set(currf,9) = 1.04651+0.0006*xp_sa(9);
				art_set(currf,10) = 1.2431+0.0002*xp_sa(10);
				art_set(currf,11) = 0.72+0.0002*xp_sa(11);
				wh = 1.04651+0.0006*xp_sa(9);
				hk1 = 1.2431+0.0002*xp_sa(10);
				g1k = 0.72+0.0002*xp_sa(11);
			elseif Npara==12
				art_set(currf,8) = xp_sa(9);
				art_set(currf,10) = 1.04651+0.0006*xp_sa(10);
				art_set(currf,11) = 1.2431+0.0002*xp_sa(11);
				art_set(currf,12) = 0.72+0.0002*xp_sa(12);
				wh = 1.04651+0.0006*xp_sa(10);
				hk1 = 1.2431+0.0002*xp_sa(11);
				g1k = 0.72+0.0002*xp_sa(12);
			end
			return;
		    end
		    if (fp >= f) 
			for ii=1:Npara,
			    xvect(ii) = xp_sa(ii);
			end
			f = fp;
			NACC = NACC + 1;
			Nacp(h) = Nacp(h) + 1;
			nup = nup + 1;
			if (fp > fopt) 
			    for ii=1:Npara,
				xopt(ii) = xp_sa(ii);
			    end
			    fopt = fp;
			    nnew = nnew + 1;
			    if iprint >=1
				[xxx,area,leng]=fcnok(Npara,xp_sa,currf,form_tar,area,leng);
            xvect_last=xvect;  %  This is to record the last 
            % xvect that is displayed with sagidraw.  This vector
            % will most likely have an error less than a vector
            % caluated for the loop Nfcnev >= maxevl.  The reason for 
            % this is that if the loop terminates early, then the vector
            % calculated may have an error greator than the previos vector.
         	 end
			end
		    else
			p = exprep((fp-f)/Temp);

				uni = u(I97) - u(J97);
				if (uni < 0.0)
					uni = uni + 1.0;
				end
				u(I97) = uni;
				I97 = I97 - 1;
				if (I97 == 0)
					I97 = 97;
				end
				J97 = J97 - 1;
				if (J97 == 0)
					J97 = 97;
				end
				rc = rc - cd;
				if (rc < 0.0)
					rc = rc + cm;
				end
				uni = uni - rc;
				if (uni < 0.0)
					uni = uni + 1.0;
				end

			pp = uni;
			if (pp < p) 
			    for j=1:Npara,
				xvect(j) = xp_sa(j);
			    end
			    f = fp;
			    NACC = NACC + 1;
			    Nacp(h) = Nacp(h) + 1;
			    ndown = ndown + 1;
			else 
			    nrej = nrej + 1;
			end
		    end
		end
	    end
	    for ii=1:Npara,
		ratio = Nacp(ii) / NS;
		if (ratio > 0.6)
		    vm(ii) = vm(ii)*(1.0+c_sa(ii)*(ratio - 0.6)/0.4);
		elseif (ratio < 0.4)
		    vm(ii) = vm(ii)/(1.0+c_sa(ii)*((0.4-ratio)/0.4));
		end
		if (vm(ii) > (ub(ii)-lb(ii)))
		    vm(ii) = ub(ii) - lb(ii);
		end
	    end
	    for ii=1:Npara,
		Nacp(ii) = 0;
	    end
	  end
	
%       if (iprint >= 1)
%           prt9(nup, ndown, nrej, lnobds, nnew);
	quit = 0;
	fstar(1) = f;
	if ((fopt-fstar(1)) <= eps_sa)
	    quit = 1;
	end
	for ii=1:NEPS,
	    if (abs(f - fstar(ii)) > eps_sa)
		quit = 0;
	    end
	end
   if (quit)
	    for ii=1:Npara,
		 	xvect(ii) = xopt(ii);
       end
     
	    fopt = -fopt;
		set(m4_1_h,'Value',xvect(1));
		set(m4_2_h,'String',round(xvect(1)));
		set(m4_3_h,'Value',xvect(2));
		set(m4_4_h,'String',round(xvect(2)));
		set(m4_5_h,'Value',xvect(3));
		set(m4_6_h,'String',round(xvect(3)));
		set(m4_7_h,'Value',xvect(4));
		set(m4_8_h,'String',round(xvect(4)));
		set(m4_9_h,'Value',xvect(5));
		set(m4_10_h,'String',round(xvect(5)));
		set(m4_11_h,'Value',xvect(6));
		set(m4_12_h,'String',round(xvect(6)));
		set(m4_13_h,'Value',xvect(7));
		set(m4_14_h,'String',round(xvect(7)));
		if (Npara==9 | Npara==12)
			set(m4_15_h,'Value',xvect(9));
			set(m4_16_h,'String',round(xvect(9)));
		else
			set(m4_15_h,'Value',0);
			set(m4_16_h,'String',0);
		end
		set(m4_17_h,'Value',xvect(8));
		set(m4_18_h,'String',round(xvect(8)));
		art_set(currf,1)=xvect(1);
		art_set(currf,2)=xvect(2);
		art_set(currf,3)=xvect(3);
		art_set(currf,4)=xvect(4);
		art_set(currf,5)=xvect(5);
		art_set(currf,6)=xvect(6);
		art_set(currf,7)=xvect(7);
		art_set(currf,9)=xvect(8);
		if Npara==8 
			art_set(currf,8)=artvelum;
		 	wh=1.34651;
		        hk1=1.3431;
			g1k=0.82081;
		elseif Npara==9
			art_set(currf,8) = xvect(9);
		 	wh=1.34651;
		        hk1=1.3431;
			g1k=0.82081;
		elseif Npara==11
			art_set(currf,8)=artvelum;
			art_set(currf,9) = 1.04651+0.0006*xvect(9);
			art_set(currf,10) = 1.2431+0.0002*xvect(10);
			art_set(currf,11) = 0.72+0.0002*xvect(11);
			wh = 1.04651+0.0006*xvect(9);
			hk1 = 1.2431+0.0002*xvect(10);
			g1k = 0.72+0.0002*xvect(11);
		elseif Npara==12
			art_set(currf,8) = xvect(9);
			art_set(currf,10) = 1.04651+0.0006*xvect(10);
			art_set(currf,11) = 1.2431+0.0002*xvect(11);
			art_set(currf,12) = 0.72+0.0002*xvect(12);
			wh = 1.04651+0.0006*xvect(10);
			hk1 = 1.2431+0.0002*xvect(11);
			g1k = 0.72+0.0002*xvect(12);
		end
		return;
	end
	Temp = rT*Temp;
	for ii=NEPS:-1:2,
	    fstar(ii) = fstar(ii-1);
	end
   f = fopt;
	for ii=1:Npara,   
	    xvect(ii) = xopt(ii);
   end 
 
end