DDA on FPGA
A modern Analog Computer

The Digital Differential Analyzer (DDA) is a device to directly compute the solution of differential equations. The DDA is a simulation of an analog computer, and is inherently parallel in operation. It can be mapped easily into an FPGA by defining mathematical integrators, adders, multipliers and other operations, then wiring them together to produce the desired dynamics. The system described here has an maximum integrator update rate of higher than 50 MHz, so it is possible to accurately simulate audio-rate signals.

For more background on analog computers see:

• Building an analog computer (Rice)
• Analog Computing (Vaxman)
• Handbook of Analog Computation (EAI)
• Fundamentals of the Analog Computer (Howe)
• Educational Analog Computer (Systron-Donner)
• Page from my freshman Physics lab notebook showing coupled oscillator analysis on analog computer

On this page you can follow these links to examples of DDAs applied to linear and nonlinear systems:

• Second order system impulse response
• Second order system driven
• Second order system with 27-bit resolution
• Second order system controlled by a NiosII cpu
• Second order system with modularized integrator and NiosII cpu

A summary of the second order system controlled by NiosII appeared in Circuit Cellar Magazine (manuscript).

Second order system (damped spring-mass oscillator):
As an example, consider the linear, second-order differential equation resulting from a damped spring-mass system:

d2x/dt2 = -k/m*x-d/m*(dx/dt)

where k is the spring constant, d the damping coefficient, m the mass, and x the displacement. We will simulate this by converting the second-order system into a coupled first-order system. If we let v1=x and v2=dx/dt then the second order equation is equivalent to

dv1/dt = v2
dv2/dt = -k/m*v1-d/m*v2

These equations can be solved by wiring together two integrators, two multipliers and an adder as shown below. In the past this would have been done by using operational amplifiers to compute each mathematical operation. Each integrator must be supplied with an initial condition.



Converting this diagram to Verilog, the top-level module verilog code defines the 18-bit, signed, state variables and a clock divider variable (count). The clocked section resets and updates the state variables. The combinatorial statements compute the Euler approximation to the F(t,V(n)). The separate multiply module ensures that the multiplies will be instantiated as hardware multipliers. The Audio_DAC_ADC module was modifed to allow either ADC-to-DAC passthru or to connect the computation output to the DAC, depending on the position of SW17. SW17 up connects the computation.

/state variables
reg signed [17:0] v1, v2 ;
wire signed [17:0] v1new, v2new ;
//signed mult output
wire signed [17:0] v1xK_M, v2xD_M ;
// the clock divider
reg [4:0] count;


//Update state variables of simulation of spring- mass
	always @ (posedge CLOCK_50) 
	begin
		count <= count + 1; 
		if (KEY[3]==0) //reset
		begin	
			v1 <= 32'h10000 ; // 
			v2 <= 32'h00000 ;
			//count <= 0;
		end
		else if (count==0)
		begin 
			v1 <= v1new ;
			v2 <= v2new ;
		end
	end
	
	// Compute new F(t,v) with dt = 2>>9
	// v1(n+1) = v1(n) + dt*v2(n)
	assign v1new = v1 + (v2>>>9);
	// v2(n+1) = v2(n) + dt*(-k/m*v1(n) - d/m*v2(n))
	signed_mult K_M(v1xK_M, v1, 18'h10000);
	signed_mult D_M(v2xD_M, v2, 18'h00800);
	assign v2new = v2 - ((v1xK_M + v2xD_M)>>>9);
	
module signed_mult (out, a, b);
	output 		[17:0]	out;
	input 	signed	[17:0] 	a;
	input 	signed	[17:0] 	b;
	wire	signed	[17:0]	out;
	wire 	signed	[35:0]	mult_out;
	assign mult_out = a * b;
	assign out = {mult_out[35], mult_out[32:16]};
endmodule

// dt = 2>>9
// v1(n+1) = v1(n) + dt*v2(n)
assign v1new = v1 + (v2>>>9)  ;
// v2(n+1) = v2(n) + dt*(-k/m*v1(n) - d/m*v2(n))
signed_mult K_M(v1xK_M, v1, 18'h1_0000);
signed_mult D_M(v2xD_M, v2, 18'h0_0800); //light damping
//scale the input so that it does not saturate at resonance
signed_mult Sine_gain(Sinput,{sine_out[15],sine_out[15],sine_out}, 18'h0_4000);
assign v2new = v2 - ((v1xK_M + v2xD_M )>>>9) + (Sinput>>>9) ; 

Two mpegs of the scope screen show the low-damping (band pass) case for fast and slow frequency steps. The top trace is the input, and the bottom is the filter output. A sudden increase in output is seen near the resonant frequency of 485 Hz. If you stop the video at a frequency below the resonant frequency, you will see that the input and output are in phase. Above the resonant frequency, the phase is 180 degrees. The whole project is zipped here.

//////////////////////////////////////////////////
//// signed mult of 3.24 format 2'comp////////////
//////////////////////////////////////////////////
module signed_mult (out, a, b);
	output 		[26:0]	out;
	input 	signed	[26:0] 	a;
	input 	signed	[26:0] 	b;
	wire	signed	[26:0]	out;
	wire 	signed	[53:0]	mult_out;
	assign mult_out = a * b;
	assign out = {mult_out[53], mult_out[50:24]};
endmodule

Second order system controlled by NiosII.
It is often useful to be able to control an analog computer using a digital computer. Such a system is called a hybrid computer. For example you might want to sequence through a set of input frequencies applied to the simulated system to build a Bode plot. A NiosII cpu was built on the FPGA with ports to (1) control the DDS frequency, (2) control the gain of the sinewave applied to a second order system, (3) control the analog reset and start the simulation, (4) record the amplitude and phase shift of the system under test. The top-level verilog module contains the DDA simulation of the second order system and the NiosII cpu. The GCC program running on the NiosII:

1. Holds the DDA in reset, and initializes the test frequency
2. Releases the DDA reset and waits 10 cycles of the test frequency
3. Waits until the output phase and magnitude reach steady-state
4. Prints the phase and magnitude
5. Iincrements the frequency by a factor
6. Repeats the above steps for a number of frequencies

An mpeg of the frequency scan shown on an oscilloscope shows the input on the top trace and output on the bottom trace. You can easily see the dramatic increase in amplitude near resonance. You can also detect the phase shift. The results printed by the NiosII are available as text for a frequency sweep. The amplitude is in output/input ratio, phase is degrees. Note that the amplitude is not very accurate at low amplitudes because the input gain was not modulated, which causes quantization errors. The results are combined into a magnitude and phase plot with this matlab routine. The vertical red line indicates the computed resonant frequency for the linear system. the red curves are the analytical solution of a second order system with the same natural frequency and damping. The phase should be -90 degrees at the red line.

The whole project is zipped here.

Second order system with modularized integrator
A clearer version of the above project can be made by factoring the integration code into a separate module. The top-level module seems more like wiring an analog computer. A bit of the top-level module shows the wiring.

	// wire the integrators
	// time step: dt = 2>>9
	// v1(n+1) = v1(n) + dt*v2(n)
	integrator int1(v1, v2, 0,9,AnalogClock,AnalogReset);
	
	// v2(n+1) = v2(n) + dt*(-k/m*v1(n) - d/m*v2(n))
	signed_mult K_M(v1xK_M, v1, 18'h1_0000); //Mult by k/m
	signed_mult D_M(v2xD_M, v2, 18'h0_0800); //Mult by d/m
	//scale the input so that it does not saturate at resonance
	signed_mult Sine_gain(Sinput,{sine_out[15],sine_out[15],sine_out}, {2'h0,Sgain});
	integrator int2(v2, (-v1xK_M-v2xD_M+Sinput), 0,9,AnalogClock,AnalogReset);

/////////////////////////////////////////////////
//// integrator /////////////////////////////////
/////////////////////////////////////////////////
module integrator(out,funct,InitialOut,dt,clk,reset);
	output [17:0] out; 				//the state variable V
	input signed [17:0] funct;      //the dV/dt function
	input [3:0] dt ;				// in units of SHIFT-right
	input clk, reset;
	input signed [17:0] InitialOut; //the initial state variable V
	
	wire signed	[17:0] out, v1new ;
	reg signed	[17:0] v1 ;
	
	always @ (posedge clk) 
	begin
		if (reset==0) //reset	
			v1 <= InitialOut ; // 
		else 
			v1 <= v1new ;	
	end
	assign v1new = v1 + (funct>>>dt) ;
	assign out = v1 ;
endmodule
//////////////////////////////////////////////////

References

• Ruben Guerrero-Rivera, et al, Programmable Logic Construction Kits for Hyper-Real-Time Neuronal Modeling, Neural Computation, Volume 18 , Issue 11 (November 2006) Pages: 2651 - 2679
• Eugene M. Izhikevich , Simple Model of Spiking Neurons, IEEE Transactions on Neural Networks (2003) 14:1569- 1572