# Simulating a simple RC circuit using an arduino and the circuit's equivalent discrete-time transfer function

Suppose we have an RC circuit with the continuous transfer function: $G(s)=\frac{1}{RCs+1}$ with R=1.6k, and C=2uF. Using the bilinear transformation(or better known as Tustin method) $s=\frac{2}{T_s}\frac{z-1}{z+1}$ with the sampling period being $T_s=0.5ms$, I arrive to the following discrete transfer function:

$G(z)=\frac{z+1}{13.8z-11.8}=\frac{X(z)}{U(z)}$

In terms of a difference equation(recurrence relation) the discrete system is thus described as $u_n+u_{n-1}=13.8x_n-11.8x_{n-1}$, where $x_n$ is the current output variable, $x_{n-1}$ is the previous output. Similarly $u_n$ is the current input variable, and $u_{n-1}$ is the previous input. The output of this system will obviously be discrete, unlike the real response of the RC circuit.

The code, I've written for the arduino (Uno model) is here:

#include "Wire.h"
#define PCF8591 (0x90 >> 1)

void AnalogOut(uint8_t value)
{
sei(); //enable interrupts
Wire.beginTransmission(PCF8591); // turn on the PCF8591
Wire.write(0x40); // control byte
Wire.write(value);
Wire.endTransmission();
}

void setup()
{
pinMode(A0, INPUT);
Wire.begin();
// timer2 initialization
noInterrupts();
TCCR2A = 0;
TCCR2B = 0;
TCCR2B = (0<<CS22) | (1<<CS20)| (1<<CS21);        //and the prescaler for a   0.5ms sampling time
TIMSK2 |= (1 << TOIE2);
interrupts();
}

uint8_t x_n_1=0;   //x[n-1]=0; initial  value of the output

uint8_t u_n;       //u[n]  current value of the input
uint8_t x_n;       //x[n] current value for the ouptut
uint8_t u_n_1;     //u[n-1] previous value for the input
//interrupt routine
ISR(TIMER2_OVF_vect)
{
x_n=uint8_t(0.07246*u_n+0.07246*u_n_1+0.855*x_n_1);   //x[n]=1/13.8(u[n]+u[n-1]+11.8x[n], also converting it to a 8bite integer
AnalogOut(x_n);                       //writing to the PC8591 d/a converter
x_n_1=x_n  ;                         //x[n-1] for the next sampling interval becomes the current x[n]
u_n_1=u_n;                                  //u[n-1] for the next sampling interval becomes the current u[n]

}
void loop()
{
u_n_1=analogRead(A0);  //getting an initial value for u[n-1]
while(1);

}


Note that I have used an interrupt routine for generating a 0.5ms delay. Also, the part of the code in the loop() function is better explained next:

 TCNT2=preload; //the timer starts ticking at its inital value
u_n_1=analogRead(A0); //we are doing this a/d conversion before even one interrupt has happened.


Here is the circuit I put together in proteus:

In proteus, I am sending $y=3+\sin(314t) volts$ to the A0 analog input. Also, I set the refference voltage for the a/d covnerter to 5V (in our actual laboratory we can set Vref using a potentiometer). Here is what I got on the oscilloscope in proteus when running the simulation: I am not sure whether this would work on the actual device since we'll be doing that in our lab in a couple of days as a part of a Digital control systems course. Is the crappy signal only a result of poor simulation or is there another error? I have tried generating a (discrete) sawtooth wave using this same circuit(but without no inputs or any transfer function, just plain signal output) and it worked well.

• You got something wrong in your calculations. I get: $u_n=0.135 x_n+0.865 u_{n-1}$ Dec 28 '16 at 19:59
• Looks like an arbitrary waveform with 64 levels to me and 40 time quantums per sine. Dec 28 '16 at 20:13
• @VladimirCravero I must have done the calculation a hundred times and I always arrive to the same result. Are you sure you aren't using another discretization method(zero order hold, IIR...)? Dec 28 '16 at 20:18
• Well I have done the calculations in another way, and yes I get your numbers. What I find odd is that for an RC filter I would not expect the output to depend on the old input, but only on the current input and the current and old output. Dec 28 '16 at 20:20
• What range of values do you expect for $u_n$?
– Chu
Dec 28 '16 at 20:55

I think it's problem of scaling the coefficients from your first-order filter equation. So, you can manipulate the values with proper precision, instead of simply converting to uint8_t the result of the sums / multiplications, as you are doing. Since at the end you will restrict the output of the filter to 8 bits (DAC) and also to avoid using 32 bit variables, I suggest to define a scale factor of 256. In this way, it will be necessary to use only 16-bit variables (uint16_t) for MAC operations. The analogRead() function returns a number from 0 to 1023 (10 bit ADC). For reasons similar to those already discussed, I will limit the result from digital conversion to 8 bits (division by 4, or right shift by 2 bits). So:

uint16_t u_n;


The original equation:

$$x_n=0.855072x_{n-1}+0.072464\left ( u_n+u_{n-1} \right )$$

Will be converted to:

$$256x_n= 219x_{n-1}+19\left ( u_n+u_{n-1} \right )$$

and implemented in code, with rounding, as:

x_n = (219*x_n_1 + 19*(u_n + u_n_1) + 128) >> 8;


In the resulting plot ($T_s=0.5 ms$) it's possible to notice that the system, when excited approximately at the frequency of the pole (50 Hz), exhibits a response with magnitude 0.707 (-3 dB) of the input and phase lag of 45 degrees. It's worth remembering that, for example, a 4V voltage is equivalent to a digital value of 4/5 * 256 (aprox. 204).

UPDATE:

You may notice in the previous plot the effect of discretization with insufficient precision, for example the presence of a DC offset error. Thus, by changing the scaling factor to 1024 and using 32-bit variables, it's possible to minimize the discretization error:

uint32_t u_n;

$$1024x_n= 876x_{n-1}+74\left ( u_n+u_{n-1} \right )$$
x_n = (876*x_n_1 + 74*(u_n + u_n_1) + 512) >> 10;