05ms612.txt
;***************************************************************************************************
; LISTING - COMPLETE CODE LISTING
;
; "Solving the loop-analysis puzzle," EDN Europe, March 2000, pg 55
; http://www.ednmag.com/ednmag/reg/2000/030200/05ms612.htm
;****************************************************************************************************

function $$loop(
   upstream_sig_name: string,
   dnstream_sig_name: string, 
   run_name : label string

),SEALED
{
   upstream_sig_name = $condition_sig_name(upstream_sig_name);
   if (upstream_sig_name == VOID) return;
   dnstream_sig_name = $condition_sig_name(dnstream_sig_name);
   if (dnstream_sig_name == VOID) return;

   	
   local first_harmonic = 250; //Hz
   local last_harmonic = 8000; //Hz
   local sim_start = .001; //1ms delay to let loop settle

   //calculate constants used in the calculations
   local num_freqs =1 + log(last_harmonic/first_harmonic)/log(2);
   local sim_end = 1/first_harmonic + sim_start;
   local sample_freq = 32*last_harmonic; 	// sample at least 32 samples of highest
 						//frequency
    
   //create local storage and file names
   local wf_index = $$get_aux_index();
   local aux_wf1_name = $strcat("aux@@Gain_",wf_index);
   local aux_wf2_name = $strcat("aux@@Phase_",wf_index);
   local Gains = $create_vector(num_freqs);
   local Phases = $create_vector(num_freqs);
   local Freqs = $create_vector(num_freqs);	
  
   //get transient data for loop input and output
   local upstream_sig_expr = $get_wvf_name(upstream_sig_name,VOID,run_name);
   local dnstream_sig_expr = $get_wvf_name(dnstream_sig_name,VOID,run_name);

   //calculate Fourier Transform:
   // F(jw)=integ(exp(-jwt)*f(t))dt
   // using Eulers Identity:
   // F(jw)=integ(cos(wt)*f(t))dt - j*integ(sin(wt)*f(t))dt
   // Setting A=integ(cos(wt)*f(t))dt and B=-integ(sin(wt)*f(t))dt
   // then |F(jw)| = sqrt(A**2+B**2)
   // and phase=arctan(B/A)
   //

   //variables used in the DFT calculation
   local freq_index; //index into frequency data
   local w;	//current radian frequency
   local f;	//current frequency
   local t;	//current time step
   local rough_t; //approximate time step
   local up;	//current upstream data point
   local dn;	//current downstream data point
   local last_t = 0; // last time step
   local dt;	//delta t
   local Aup=0, Adn=0, Bup=0, Bdn=0;
   local UpDc=0; local DnDc=0;

   //calculate DC value
   $clear_message();
   $message($strcat("Calculating DC values:,@note);
   for (t=sim_start; t<=sim_end; t=t+1/sample_freq) 
    	{
    	dt = t-last_t;
    	up = $get_signal_value(upstream_sig_expr, t);
    	dn = $get_signal_value(dnstream_sig_expr, t);
    	UpDc = UpDc + up*dt;
    	DnDc = DnDc + dn*dt;
    	last_t = t-sim_start;
	}
 
   UpDc=UpDc/(sim_end-sim_start);
   DnDc=DnDc/(sim_end-sim_start);	

   //calculate Fourier Transform
   
   for (freq_index=0; freq_index<num_freqs; freq_index=freq_index+1)
	{
	f = first_harmonic * pow(2,freq_index);
	$clear_message();
	$message($strcat("Processing Frequency: ", $format("%6.6f", f), "."),@note); 
	w = 2*pi*f;
	Aup=0; Adn=0; Bup=0; Bdn=0;
	last_t = 0;
	for (rough_t=sim_start; rough_t<=sim_end; rough_t=rough_t+1/sample_freq) 
	   {
	    local up_data =
$get_signal_transitions(upstream_sig_expr,@around,rough_t,,1,@horizontal); 
   	    if (up_data == UNDEFINED) return;
   	    local dn_data = 
$get_signal_transitions(dnstream_sig_expr,@around,rough_t,,1,@horizontal); 
     		 if (dn_data == UNDEFINED) return;
	    t = up_data[0][1] - sim_start; //get exact time value;
	    up = up_data[1][1]; //get the first datapoint after t
	    dn = dn_data[1][1];
	    dt = t - last_t; //use the real time point, not t
            last_t = t;
	    Aup = Aup + cos(w*t)*up*dt;
	    Bup = Bup + sin(w*t)*up*dt;
	    Adn = Adn + cos(w*t)*dn*dt;
	    Bdn = Bdn + sin(w*t)*dn*dt;
	    }

	local int_per = last_t; // *2 because area under cos*cos=.5*f 
	Aup=Aup*2/int_per; Bup=Bup*2/int_per; Adn=Adn*2/int_per; Bdn=Bdn*2/int_per;
	Freqs[freq_index] = f;
   	//one way to calculate gain and phase
   	//Gains[freq_index] = 20*log10(sqrt((Aup*Aup+Bup*Bup)/(Adn*Adn+Bdn*Bdn)));
   	//Phases[freq_index] = deg(atan((-Bup/Aup)/(-Bdn/Adn)));
	//another way to calculate gain and phase
	local Fup = $create_complex(Aup,-Bup);
	local Fdn = $create_complex(Adn,-Bdn);
	Gains[freq_index] = 20*log10($magnitude_complex($divide_complex(Fup,Fdn)));
	Phases[freq_index] = deg($phase_complex($divide_complex(Fup,Fdn)));
   	//test mode, for measuring gain and phase of perturbation
	//Gains[freq_index] = 20*log10($magnitude_complex(Fup)+1e-12);
	//Phases[freq_index] = deg($phase_complex(Fup));
	}


 
 //Faster Version

function $$loop_Fast(
   upstream_sig_name: string {default=$$get_selected_signal(1, @netpin)},
   dnstream_sig_name: string {default=$$get_selected_signal(2, @netpin)},
   run_name : optional string,
   first_harmonic: real {default = 125}, //Hz
   last_harmonic: real {default = 8000}, //Hz
   sim_start: real {default = .005} //1ms delay to let loop settle

),SEALED
{

   // condition signals 
   upstream_sig_name = $condition_sig_name(upstream_sig_name);
   if (upstream_sig_name == VOID) return;
   dnstream_sig_name = $condition_sig_name(dnstream_sig_name);
   if (dnstream_sig_name == VOID) return;

   local num_freqs = 1 + log(last_harmonic/first_harmonic)/log(2);
   local sim_end = 1/first_harmonic + sim_start;
   local samples = 1024; // sample at least 1024 samples per point
    
   //create local storage and file names
   local wf_index = $$get_aux_index();
   local aux_wf1_name = $strcat("aux@@Gain_",wf_index);
   local aux_wf2_name = $strcat("aux@@Phase_",wf_index);
   local Gains = $create_vector(num_freqs);
   local Phases = $create_vector(num_freqs);
   local Freqs = $create_vector(num_freqs);	
   local Aup=$create_vector(num_freqs);
 	local Adn=$create_vector(num_freqs); 
	local Bup=$create_vector(num_freqs); 
	local Bdn=$create_vector(num_freqs);
  
   //get transient data for loop input and output
   local upstream_sig_expr = $get_wvf_name(upstream_sig_name,VOID,run_name);
   local dnstream_sig_expr = $get_wvf_name(dnstream_sig_name,VOID,run_name);

   //calculate Fourier Transform:
   // F(jw)=integ(exp(-jwt)*f(t))dt
   // using Eulers Identity:
   // F(jw)=integ(cos(wt)*f(t))dt - j*integ(sin(wt)*f(t))dt
   // Setting A=integ(cos(wt)*f(t))dt and B=-integ(sin(wt)*f(t))dt
   // then |F(jw)| = sqrt(A**2+B**2)
   // and phase=arctan(B/A)
   //

   local freq_index; //index into frequency data
   local w;	//current radian frequency
   local f;	//current frequency
   local t;	//current time step
   local rough_t; //approximate time step
   local updt;	//current upstream data point
   local dndt;	//current downstream data point
   local last_t = sim_start; // last time step
	local dt;	//delta t
	local UpDc=0; local DnDc=0;

	//calculate DC value
 	$clear_message();
	$message("Calculating DC",@note);
	for (t=sim_start; t<=sim_end; t=t+1/(samples*first_harmonic))
	    {
	    dt = t-last_t;
	    local up = $get_signal_value(upstream_sig_expr, t);
	    local dn = $get_signal_value(dnstream_sig_expr, t);
	    UpDc = UpDc + up*dt;
	    DnDc = DnDc + dn*dt;
	    last_t = t;
	    }
  
	UpDc=UpDc/(sim_end-sim_start);
	DnDc=DnDc/(sim_end-sim_start);	

	for (freq_index=0; freq_index < num_freqs; freq_index=freq_index+1)
						{
						Aup[freq_index] = 0;
						Bup[freq_index] = 0;
						Adn[freq_index] = 0;
						Bdn[freq_index] = 0;
						}

   //calculate Fourier Transform
last_t = 0;
	for (rough_t=sim_start; rough_t<=sim_end; 
rough_t=rough_t+1/(samples*last_harmonic)) 
	       {
	        //get data
	        local up_data =
$get_signal_transitions(upstream_sig_expr,
@around,rough_t,,1,@horizontal); 
   	        if (up_data == UNDEFINED) return;
   	     local dn_data = 
$get_signal_transitions(dnstream_sig_expr,
@around,rough_t,,1,@horizontal); 
     		     if (dn_data == UNDEFINED) return;
		     t = up_data[0][1] - sim_start; //get exact time value;
	        dt = t - last_t; 
//   		  $clear_message();
//			  $message($strcat("Computing time point: ",
//						$format("%6.12f",t),
//						"."),@note); 
           last_t = t;
		     updt = up_data[1][1]*dt; //get the first datapoint after t
	 	     dndt = dn_data[1][1]*dt;

			  //perform summation
	        if (dt != 0)
		         {
				f = first_harmonic;
  				w = 2*pi*f;
				local cosine = cos(w*t);
				local sine = sin(w*t);
   				for (freq_index=0; freq_index < num_freqs;
freq_index=freq_index+1)
					{				
					Aup[freq_index] = Aup[freq_index] + cosine*updt;
					Bup[freq_index] = Bup[freq_index] + sine*updt;
					Adn[freq_index] = Adn[freq_index] + cosine*dndt;
					Bdn[freq_index] = Bdn[freq_index] + sine*dndt;

					sine = 2*sine*cosine;  //double angle rule
					cosine = 2*cosine*cosine -1;

					}

			      }
        }

       // divide the summations by 1/T
	    local int_per = last_t; // *2 because area under cos*cos=.5*f 
	    Aup=Aup*2/int_per; Bup=Bup*2/int_per;
 Adn=Adn*2/int_per; Bdn=Bdn*2/int_per;

  	    //calculate gain and phase and store results
   		for (freq_index=0; freq_index < num_freqs; freq_index=freq_index+1)
			{
		local Fup = $create_complex(Aup[freq_index],-Bup[freq_index]);
	   	local Fdn = $create_complex(Adn[freq_index],-Bdn[freq_index]);
	    	Gains[freq_index] = 
20*log10($magnitude_complex($divide_complex(Fup,Fdn)));
	    Phases[freq_index] = deg($phase_complex($divide_complex(Fup,Fdn)));   
	    Freqs[freq_index] = first_harmonic * pow(2,freq_index);
	    }

   //store the new curves to auxiliary waveforms
   $create_double_waveform(aux_wf1_name,[Freqs,Gains],"frequency",,@replace);
   $create_double_waveform(aux_wf2_name,[Freqs,Phases],"frequency",,@replace);

   //then chart them
   $$chart_result(@NEW,@stacked,"A",aux_wf1_name,,"Loop Gain",,"LOOP ANALYSIS");
   $change_axis_title("Loop Gain (dB)","Y:Unknown");
   $$chart_result(@SELECTED,@stacked,"A",aux_wf2_name,,"Loop Phase");
   $change_axis_title("Loop Phase (degrees)","Y2:Unknown");
   $keep_window_active();

}