Jun 30, 2015

System parameter estimation with Simulink

In previous blog entries, I created custom Simscape components for TEC and forced air heatsink for a thermal control system.  I built an open-loop system model with a current source, and now I want to find reasonable values for the model's parameters using openloop data.  While I've done this exercise many times in the past, I always did it within Matlab itself rather than in Simulink.  In an effort to catch up to the modern Matlab workflow, I will document my new work procedure.

Create input/output ports for the model

I already have my model (creatively called "take1") open, but it is currently missing output signal.  I dragged an out1 port from the Simulink Library Browser --> Simulink --> Sinks, and changed the port property to show port name and number.  Similarly, I dragged an in1 port from the Library Browser and connected to the ideal current source in my model.

Bring the openloop data into Matlab workspace

My data is already in a CSV file in my Matlab folder (~/Documents/MATLAB), so I can easily bring it into the workspace by right clicking on the file in the Current Folder view  --> Import Data... --> Import.  Now my Matlab Workspace is choke-full of column vectors.

I believe the csvread() command would work just as well.  In there, I identify the vectors I want to use for this experiment:
  • input1: TECValue, commanded value to the current generating FW
  • output1: CenterTemp
  • ouput2: I have several thermistor readings, which I should average, like this:

    >> center = mean([Sensor5 Sensor6 Sensor9 Sensor10], 2);
    >> center = center(4215:6000); % Just deal with relevant data

  • output3: HeatsinkTemp
The data was captured at 2 Hz, but the time was not recorded.  So I just make it up with this Matlab command:

>> time = (0 : 0.5 : 0.5*(size(center,1)-1))';

Note the transpose (') operator to keep the time as a column vector, just like the other openloop data columns.  I suppose I could have defined a sampling frequency as a constant, but let's just keep moving for now.

Create an experiment

In the Simulink model window, select Analysis > Parameter Estimation.  This action opens a new session, Parameter Estimation - take1, in the Parameter Estimation tool.  Click "New Experiment" in the toolbar to bring up "Edit Experiment: Exp" wizard--"Exp" is the default name of the experiment.
  • Select all output ports I created in the model, and supply the openloop data
    • output 1: [time, CenterTemp(4215:6000)]
    • output 2: [time, center]
    • output 3: [time, HeatsinkTemp] I ignored the heatsink temperature because it is actively controlled.
  • In the "Inputs" category, select only some of the rows: [time, TECValue(4215:6000)]
  • Select Initial States. Optional?
  • Select parameters: Simulink found all parameters in my model: those that were left in symbolic form.  I checked all of them.  If necessary (usually a good idea), the parameters should be constrained to sensible bounds. when expanded, the parameter tab shows the minimum and maximum bounds.  For almost all physical parameters, the minimum should be some small value.

Clean up data

To view the imported data, click  Add Plot on the Parameter Estimation tab and select the experiment name--"Exp" here.

Since the open loop data covers many different cases, here is a way to restrict the data range: in the "EXPERIMENTAL PLOT" tab, click "Extract Data" button, and specify the "Start/End Time", as shown below:
Then I clicked "Save As", which created another experiment "Exp1".  I renamed Exp1 to "Exp TEC".  In the EXPERIMENT PLOT tab above, there are other ways to clean up the experiment data:
  • Remove outlier, offset
  • Scale, filter data

Estimate parameters

When I click on the "Estimate" button shown below, the GUI runs optimization loop, and adjusts the parameters until the score (sum of the squared error) stops converging.
The EsitmatedParams "Scales Value" shows that my initial guess for one of the parameters was particularly poor.  Going back to the Experiments browser --> Edit --> Parameters, the parameter values at the end of the estimation iteration are shown, and I can just read them off.

Clicking on "Plot Model Response" button above overlays the simulation result against the experimental data.
The fit is not all that great, suggesting that my simple model is missing at least 1 first order dynamics between the TEC and the out1, and maybe the TEC data I have is somewhat off.  Rather than complicate the model further, I let the parameter estimator vary the TEC parameters to obtain a better fit, as shown below:

If I had more openloop data, I could validate the model + parameter against it.

Next step

I have to design a feedback and feed-forward controller for this system.   The feedback controller will be relatively easy, but the feed-forward tends to have lots of heuristics and tuning knobs.

Jun 29, 2015

Modeling a forced air heatsink in Simscape

In a previous blog post, I simulated a temperature controlled device that uses TEC to move heat in either direction.  Another common practice in thermal control is to thermally glue a heatsink on the "cold side" of the TEC.  When cooling an area thermally glued to the "hot side" of the TEC, the "cold side" gets hot, and reduces the cooling efficiency; if the "cold side" continues to get hotter, the TEC will go into a runaway condition.  The solution is to increase the effective surface area of the TEC cold side with a heatsink.  But as you see in your own laptop, natural convection occurring on the larger surface area alone cannot move the many Watts the modern TEC can move and need the help of forced convection cooling; air cooling is cheaper than liquid cooling.

Forced convection constant is a non-linear function of the fan speed

Regardless of whether forced or unforced, all heatsinks have thermal mass (capacitance).  So an easy way to think about the problem at a high level is to consider the change in the thermal resistance from the heatsink fins to the ambient (air if air cooled, or liquid if liquid cooled).  If I plot Tsink - Tambient (the temperature of the heatsink relative to the ambient) with DC fan at max speed, I may see a classical RC circuit cooling behavior, like this:
The expotential decay portion is the exp(-t/RC) curve, so I can plot ln(Tsink - Tambient) and read off the negative slope to estimate the 1/RC.  Intuitively, higher DC fan speed yields higher 1/RC--faster cooling, BUT 1/RC does NOT scale linearly with the fan speed.

Heatsink convection vs. air speed

A simplified model of convection is Q = h (Tsink - Tambient), and the problem boils down to estimating the convection constant for various conditions.  I remember the scary integrals in the heat transfer class, which are often done in FEA these days.  I didn't do well in fluid mechanics or heat transfer, but was able to follow along a heat transfer textbook: Engineering Heat Transfer 2nd Edition, by Rathore and Kapuno.  On page 563, there is a heatsink fan configuration that looks similar to mine (array of cylinderical rods facing the air), so I started there.  For forced air convection, the h term is Nu ∙ kf / D, where
  • kf is the thermal conduction coefficient of air film around heat fin.  Typically, we would just look up the value at an average temperature of (Tsink + Tambient)/2.  Between 26.3 ~ 30E-3 W/m∙K for 300~350 K.
  • D is the rod diameter, say 2 mm.
  • Nu is the Nusselt number that depends on Reynolds number (which in turn depends on the forced air speed).  For my configuration, the textbook says it is 0.27 ∙ (u_air D / ν)^0.63 ∙ Pr^0.36 ∙ (Pr∞ / Prsink)^1/4, where
    • u_air is the bluk air speed through the heatsink,
    • ν is the air viscosity, ranging 15.89 to 20.92 m2/s for air between 300 to 350 K,
    • Pr is the Prandtl number: 0.707 at 300 K and 0.700 at 350 K.
    • The Nu therefore works out to something like 0.7E-3 u_air^0.63
Note that Nu depends the most on the air speed uair, namely, non-linearly as a power of 0.63.  With this Nu, you are supposed to solve for the temperature of the air leaving the heatsink fins: Texit = Tsink – (Tsink – Tambient) exp[-π D N h / (ρ u_air N_T p_T C_p)], where:
  • N is the total number of fins
  • h was calculated above,
  • N_T is the number of rods in the transverse plane
  • ρ is air density, about 1 Kg/m3 around 300~350 K
  • p_T is the fin spacing (pitch)
  • C_p is the air's heat capacity (about 1 J/g∙K)
After substitution, you arrive at Texit = Tsink – (Tsink – Tambient) exp[Ca / u_air0.37], where C_a is depends on your situation.  What this equation says is that the exit temperature at the heatsink is almost the Tambient at all air speeds except the really slow case, in which case the exit temperature is almost the Tsink, driven by this example of the exponential term:
Intuitively, convection from a heatsink is NOT from Tsink, but rather a temperature that is closer to Tambient.  How close depends on the heatsink shape and air speed.

Finally, what you really want to know is the aggregate heat flow from the heatsink to the ambient air.  For the N number of rod case, Qconvect = N (h π D L) ∆T = π N ∙ Nu ∙ kf ∙ L ∙ ∆T, where:
  • ∆T is calculated above--practically speaking something like 90% of (Tsink - Tambient),
  • L is the rod height
Since Nu is proportional to (u_air D / ν)^0.63, Qconvect ~ some_number ∙ u_air^0.63 ∙ (Tsink – Tambient) for all but the very low air speeds.

Now that you understand, keep the simplest 1 equation

If the textbook was completely applicable for my case, I should be able to predict the 1/RC slope for all forced air cooling based on the airspeed--which according to the DC fan manufacturer scales linearly with the fan speed.  But instead of u_air^0.63, I found a factor of u_air^0.8 to be a better fit.  So instead of losing sleep over the discrepancy, I will just use the power of 1/0.8 = 1.25 when calculating the desired fan speed.  That is, my forced air heatsink is simply characterized by 
Qconvect(duty) = [h@100% * duty^0.8 + h@0%] (Tsink - Tambient),
where 0 <= duty <= 1 ; that is, consisting of the forced and natural convection terms.

Custom Simscape component for the forced air heatsink

I explained the steps to creating a custom Simscape component in a previous blog.  Here is the Simscape code for the above heat flux equation:

component Heatsink
% Heatsink: 2
% If the duty cycle is s, where 0 <= s <= 1,
%   q = Csink (h100 * s^0.8 + h0) (H - C)
%
% where the following parameters are treated as constants for this simple
% model:
%
% * h100: convection constant at max fan duty cycle [W/K]
%
% * h0: natural convection constant [W/K]

nodes %conserving ports
H = foundation.thermal.thermal; % H:top 
C = foundation.thermal.thermal; % C:bottom
end

parameters
h100 = {0.032, 'W/K'};
h0 = {0.0019, 'W/K'};
end
inputs
s = {0, '1'};% s:bottom
end

variables(Access=private)
q = { 0, 'W' };  % heat flux
end

function setup
if h0 <= 0
pm_error('simscape:GreaterThanZero', 'h0')
end
if h100 < 0
pm_error('simscape:GreaterThanZero', 'h100')
end
end

branches % Relationship between the Through variables
q : H.Q -> C.Q;%In heatsink, heat flows hot to cold
end

equations % Relationship between the Across variables
q == (h100 * s^0.8 + h0) * (H.T - C.T);
end % equation
end %component

And while I may yet get tired of such frivolities, I got a picture of a heatsink from Bing Images and put it right next to the source code during the compilation of my simscape library, so that this custom component may show up as a picture.

Simulate cooling a hot block with the heatsink

As a sanity test, set the initial temperature of the heatsink at 100 C, when turn the fan on at 50% duty cycle.  Here is the Simscape block diagram for it:
And here is the 10000 second simulation result:
As you can see, this heatsink is not all that powerful.

Jun 28, 2015

Simscape: look ma, (almost) no equations!

I studied control engineering in school, but like many people, I found only a few occasion to directly use what I learned in school.  If I had attended this lecture live while I was at Caltech, I might have been more inspired and finished my Ph.D in control and mechanical engineering, but what is the use of thinking about "what might have been" 20 years earlier?

What little control engineering I got to do in my otherwise mostly an embedded software engineering career has been in environment control: temperature, humidity, air flow, etc.  While I still dabble in other control applications like power, motor, pneumatic, there is a universal need for this HVAC-like applications in biotech industry, and therefore pays the bills.  The control engineering method I learned in school starts with simple schematic-like diagram of the device under control, and the simplest possible relevant governing equations: for example temperature and heat flow in thermal domain, which works out to a bunch of RC-circuit like nodes.  While a single first order system is a control engineering 101 lecture material, the real world is rarely that simple: if you put a few such first order systems together, sufficiently complex behavior will result.  To even describe such system, the governing equation balloons out due to the number of states, and manually manipulating the equation before deriving a LTI (linear time invariant) system equation gets unwieldy until all the equation can be beautifully consolidated to the final state-space form: dX/dt = A X + B U; Y = C X.  But since the controller is implemented in a digital form, the 1st order differential equation is then promptly converted to the difference equation: X[k+1] = A X[k] + B U[k].  The linear algebraic formalism allowed advanced MIMO (multi-input-multi-output) controllers that can control the system from incomplete/noisy measurements.  For example, this is a block diagram of closed loop feedback regulator (tries to achieve the commanded value of the output) with an estimating observer:
A standard way to deal with noisy measurement is Kalman filter.

While great in theory, I found that complex systems (systems with say 5 states or more) are difficult to control using these LTI controllers, and made simplifications reduce the number of states before I start working with the state space representation in my ancient copy of Matlab controls toolbox, running on Windows XP.  But thanks to MSFT's finally killing the Windows XP, I just upgraded to Matlab and Simulink 2015a, and found many new things, including the Simscape--which lets you DRAW the system model (look ma, no equations!) and simulate the response.  It almost feels like cheating, but since I am not in school any more, I will explore in this blog entry how control engineering has become slightly easier since I left school.

TEC (thermoelectric cooler) Simscape model

In biotech, TEC is the temperature cooling device of choice, because the heat transfer requirement is relatively low and its non-moving part nature.  I was hoping that Simscape already includes the TEC model--if not in the fundamental models, perhaps in the SimPower or SimElectronics, but alas, it did NOT.  I ran into an article where the author built a Simscape model [Caroff, T., et al. "Transient cooling of power electronic devices using thermoelectric coolers coupled with phase change materials." 2013 19th International Workshop on Thermal Investigations of ICs and Systems (THERMINIC)] but did NOT release that Simscape model to the public domain.  So I wound up creating a custom Simscape model--not exactly a gentle introduction into Simscape.

Peltier governing equation

A TEC can be driven with either positive or negative current (with respect to the factory indicated ground terminal--usually colored black); it will move heat from one side of the TEC to the other depending on the polarity of the current.  In dealing with TEC, I found it helpful to adopt a notion that the side that gets hot when positive current is applied is the "hot side", and the heat flow direction is from positive in the cold to the hot side (yes, it's counter-intuitive, but you will see it is notation-wise simpler).
Q is a non-linear function consisting of 3 terms:
  1. Peltier effect, driven by the current, which is non-symmetric across current = 0
  2. Self-heating due to current flow (resistive heating), driven by the current
  3. Conduction from the hot to the cold side
The governing equation for the TEC's heat moving capacity is:

Q = Sm Th Itec + 0.5 Rm Itec^2 + Km (Tc - Th) for Itec >= 0
    = Sm Tc Itec + 0.5 Rm Itec^2 + Km (Tc - Th) otherwise

Staring at the equation, I am relieved that Q is continuous at Itec = 0.  The thermal gradient is positive from Tc to Th, as I had established in the convention above.  It also makes sense that the electrical self-heating should generate heat inside the TEC.  To cool (Q < 0), Itec must be < 0, and has to overcome the resistive heating and the thermal gradient (because the Tc is usually > Th when cooling).  Clearly, the TEC efficiency will suffer drastically if the "cold side" is not adequately cooled with a heatsink.  Because the self heating grows as quadratically with the current, a TEC will eventually reach a runaway condition if it pulls a lot of current for extended period.

The Peltier effect also works the other way, and produces a "back-EMF"-like voltage across the TEC's positive and negative terminals (V+ - V-): V = Sm (Th - Tc) + Rm Itec, which can be rearranged to Itec = [V - Sm (Th - Tc)] / Rm.

A custom Simscape component for the TEC

Creating my first Simscape package

The Simscape documentation said:
Simscape files must reside in a +package directory on the MATLAB path.
After a fresh install, Matlab ~/Documents/MATLAB and put that in the path.  So I just created my own Simscape package folder there:

~/Documents/MATLAB$ mkdir +henryscape

Here, I created a new file TEC.ssc and wrote the following code:
component TEC
% Thermoelectric cooler: 2.0
The first comment line after "component TEC" is the "user friendly name" to be displayed in the model view.  "2.0" after the colon 
% Move heat with electrical work.
% More comments not shown in the interest of space

The rest of the comment shows up as the component explanation.
nodes %conserving ports
p = foundation.electrical.electrical; % +:top
n = foundation.electrical.electrical; % -:bottom
Th = foundation.thermal.thermal; % H:top 
Tc = foundation.thermal.thermal; % C:bottom
end
These ports show up in the block diagram.  The comments immediately after the domain declaration controls the displayed name and location of the ports in the block diagram
parameters
S = {0.08, 'V/K'};  % Siebeck constant
K = {2, 'J/(K*s)'}; % thermal conductivity; should use [W]?
R = {2.2, 'Ohm'};   % electrical resistance
end
Parametrizing the component's property enhances reuse.  In my TEC component's case, the 3 parameters were captured in the governing equation given above.  The default values are for the TEC I happen to be using, but easily changeable when I double-click on the component in the model editor.  I hesitated a bit when picking the unit for the thermal conductivity: Watts per Kelvin [W/K] would have been more concise, but I stuck with [J/(s*K)] because I saw the DC motor component example use [J/s] when it could have used [W].
variables(Access=private)
i = { 0, 'A' };  % Current
v = { 0, 'V' };  % Voltage
qK = { 0, 'J/s' };% conductive heat flux
 qR = { 0, 'J/s' };% resistive heat generation
 qS = { 0, 'J/s' };% Siebeck heat flux
end

I actually don't know yet the reason for choosing different access level; just copied from a Simscape component example.  Could I have used 'W' instead of 'J/s'?  Yes.
function setup
if S <= 0
pm_error('simscape:GreaterThanZero', 'S')
end
if K <= 0
pm_error('simscape:GreaterThanZero', 'K')
end
if R <= 0
pm_error('simscape:GreaterThanZero', 'R')
end
end

Parameter validation.
branches
i : p.i -> n.i;
qK : Th.Q -> Tc.Q;
 qR : * -> Th.Q;
 qR : *-> Tc.Q;

 qS : Th.Q -> Tc.Q;
end

qK (conduction always goes from hot to cold).  Keep the sign convention this way and also keep your sanity.

Resistive heat is generated from inside the TEC and flows outside.

Peltier heat flows from the hot side to the cold.

branches specifies the relationship between the Through variables
equations
v == p.v - n.v;
v == S * (Th.T - Tc.T) + R * i;
qK == K * (Th.T - Tc.T);
 qR == 0.5 * R * i * i;

if ge(Th.T, Tc.T)
qS == S * Th.T * i;
else
qS == S * Tc.T * i;
end
end % equations

Equation specifies the relationship between the Across variables.

qR is always NEGATIVE so that heat always flow from the TEC to the outside.  But this is OPPOSITE of the above definition.

qS is based on the HOTTER side temp.  It is NON-linear.

"let" keyword would allow me to describe TEC heat flow with intermediate terms for comprehension.  end keyword is required at the end of a let-in-end statement.
end %componentI miss the curly braces...

Compiling my custom package

I also wanted a TEC icon for the block, so I found a TEC image (below) on the web and put it right next to TEC.ssc.

I like this picture because it shows the heat "going up" as I've written the equation.  To actually generate the package, I ran this command in the Matlab shell:

>> cd ~/Documents/MATLAB
>> ssc_build henryscape
Generating Simulink library 'henryscape_lib' in the current directory '/home/henry/Documents/MATLAB' ...

After this, I find a file henryscape_lib.slx in ~/Documents/MATLAB folder, and when I double-clicked on the lib file in Matlab's folder, I saw the component appear in the Simulink libarry viewer:
Note that the + and the H nodes are on top, and the - and C nodes are on the bottom, as intended in the source code.

Using the TEC component in a Simscape model

I create a new Simscape model with this Matlab command:

>> ssc_new

This populates a new model with components that every Simscape model needs, as shown below:
Note the model is using the recommended solver for Simscape: ode23t.

I then dragged the TEC component from the library viewer shown in the previous section into this model, and then added an ideal current source controlled by a square wave, and measured the temperature on the hot side of the TEC, using this model:
The resulting plot of the Th - Tc (the temperature difference between the TEC hot and the cold side) is shown here, overlaid against the +/- 1A current being driven from the current source:
Even with an ideal current source, the Peltier effect prevents the current from changing discontinuously.

Next step

I will connect this TEC component to the part that I am trying to heat/cool, and also to the heatsink that sits under the TEC.  I saw a Simscape tutorial that fits parametrized model to openloop data.  In that scenario, some of the parameters are left in symbolic form that can be changed from Matlab script.  For example, a thermal capacitance (mass) of the part I am trying to cool is called wallMass in the thermal mass block property, as shown below:
Ultimately, I will throw in a real-time controller that acts on the measured temperature of the part that I am trying to heat/cool.

Jun 24, 2015

Bare metal code to read ADC on Zynq

In a previous blog entry, I studied interacting with the on-chip ADC peripheral using the Xilinx XADC device driver.  I was able to to understand device driver, and get results out, but I prefer real-time HW control, with minimum code between the HW and the data collection.  Xenomai patch will yield a nearly real-time behavior (as I reported in a real-time jitter test of various Linux options), and even though Linux device driver framework is a wonderful SW (literally; whenever I stare at the Linux kernel code, I think: "how did they pull all of this off?"), I just don't want to run all that code just read a few bits, when the HW itself is relatively simple.

Aside from my minimalist preference, there is a concrete benefit to bare metal code: certifiability; when I was in the defense industry for a few years after school, the software certification cost was running at about $100~$150 per LOC.  I am sure it is more now.

Status (so you don't have to read the whole thing)

  • Able to read the dedicated analog input and 15 other auxiliary channels in channel sequencing mode, kicked off from a timer interrupt (to turn off XADC when not in use).
  • Have a thermistor circuit that can connect to the Zedboard.
  • Reading out the XADC seems to take too long, even though it's just a bunch of memory mapped register read/write operations.  I think it may be due to the XADC wizard performing a read over the DRP FIFO internally...

XADC peripheral

To recap, the 2 XADC peripherals on Zynq is (as described in the Zynq TRM) a 12-bit 1M sample/s HW, with up to 17 external analog channels each.  There is a not-so-high-quality on-chip reference voltage is good enough for the internal temperature and voltage sensors, but the TRM recommended an external 1.25 V reference IC.  The HW has these potentially useful features:
  1. Remembering the maximum/minimum values in dedicated registers--until the next XADC reset.
  2. Interrupting based on threshold breach--so that the SW filtering can be bypassed.  The maskable alarm is latched in the XADCIF_INT_STS register, and is cleared by writing to the same register.
  3. Averaging in HW: if the system's anti-aliasing filter BW is higher than the application's sampling BW, sampling the ADC fast and then averaging in HW may be the only option.
Q: why does UG480 say: "auxiliary analog inputs do not require any user-specificied constraints or pin locations...  All configuration is automatic when the analog inputs are connected to the top level of the design."?

ADC sampling is divided into 3 phases:
  1. Acquisition: the capacitor in the ADC HW is charged from the input source.  Depending on the input impedance, this may take longer than the default 4 ADCCLK.
  2. Conversion: the charge in the capacitor is converted to digital value, and written out to the appropriate register.
  3. Readout: SW should read out from the register--unless a streaming arrangement is made ahead of time.
In sequencing mode, the acquisition for the next channel is overlapped with the conversion of the current channel (the pipeline concept).  If the acquisition requires more time, the conversion for the next channel can be made to start 10 ADCCLKs later by asserting the ACQ in the control register.

Control interface

The peripheral can be accessed through 2 bus (you have to choose):
  • PS-XADC: directly from the processor, through the APB bus.  Commands are serialized through a 15-word deep FIFO, and against through a read FIFO on the way back.  This is a lower speed interface, but more convenient than setting up the AXI channel to the XADC. 
  • PS-to AXI ADC: through the AXI bus.
It appears that to enable the PS-XADC interface, the XADCIF_CFG[ENABLE] bit has to be asserted, to enable the interface arbiter to choose the PS-XADC (instead of PL-JTAG).  The serial interface is clocked by the PCAP_2x clock from the PS clocking subsystem (nominally 200 MHz), but divided by 4 * XADCIF_CFG[TCLK_RATE], to guarantee the serializer clock rate to maximum of 50 MHz. The interface is complicated slightly by its SPI-like behavior: the HW pushes out data to the read FIFO (read through the XADCIF_RDFIFO register) only when the write FIFO is written to, so that a NOOP command (XADCIF_CMD[29:26] = 4'b0000) is necessary to read the value of a single read.  In a continuous read case, such NOOP would not be required.  I can see that concretely in the XSDK BSP generated examples.

XSDK BSP generated PS XADC example

Q: why is the INTR_ID in the generated code 39?
A: BSP thinks I am using sysmon feature to talk to the XADC.  xparameters.h:

#define XPS_SYSMON_INT_ID 39

It also thinks my base address is different than what I was assigned!  XADC inteferface configuration is apparently not the same as the AXI 4 Lite interface to the XADC wizard IP.

#define XPAR_PS7_XADC_0_BASEADDR 0xF8007100

This register is shown in the Zynq TRM, XADCIF_CFG, in section B.16, Device Configuration Interface.  What is confusing about this at first is that this is NOT the direct interface to the XADC HW itself, but rather the SPI-like FIFO to push the commands to.  To control the XADC from the CPU, several writes to the device configuration register blocks were required:
  1. Write a magic value (0x757BDF0D) to the unlock register (offset 0x34 in the configuration register base address of 0xF807000).  This write is a simple memory mapped register write.
  2. And then the example writes to the XADCIF_CFG register (0xF8007100), turning on the following bits:
    1. ENABLE
    2. CFIFOTH = 0xF
    3. DFIFOTH = 0xF
  3. Clear the miscellaneous control (MCTL), to release XADC reset
To self-test the XADC:
  1. Reset XADC by bouncing the reset bit in MCTL register just discussed above.
  2. Write some values into the alarm threshold (like VCCINT_UPPER, which is at 0x51 in the XADC internal register banks) and read it back.  This tests the command FIFO to the XADC.  Internally, this is what happens during write/read:
    1. Write
      1. 32 bit message is formatted into this 32-bit JTAG message
        1. MSB: write (0x08) or read (0x04)
        2. register address left shifted by 16 bits (address bits are position [25:16])
        3. 16-bit data in the 2 LSB
      2. Write into the command FIFO (@ 0xF8007110).
      3. Read the Read FIFO (@ 0xF8007114) after any write since for each write, one location of Read FIFO gets updated.  Throw the read value away.
    2. Read
      1. Create a dummy data (0) to write into the write FIFO: 0x04000000 | (reg offset << 16)
      2. Write into the CMDFIFO
      3. Read from the RDFIFO TWICE, and keep only the 2nd read
To configure the XADC for real now,
  1. Set sequencer mode to the default mode (value 0 to CFG1[15:12]), to prevent alarms from tripping.
  2. Disable all alarms (value 0 to CFG1[11:0]).
  3. Restore the on-chip temperature and voltage alarms.
  4. Register SCU interrupt handler.  This interrupt handler reads the XADC int status register @ 0xF8007104 (which only tells you about the FIFO threshold, overtemp, or alarm)--these are not what I am interested in.
  5. Change the sequencer mode again.
  6. Read internal registers, like temperature (0x00) and Vccint (0x1)
While it IS possible to control the XADC HW directly through the built-in control path in the PS, it appears that going through the memory mapped XADC wizard register will be simpler (next).

XSDK BSP generate sysmon example

To reset the XADC wizard IP (I guess Xilinx calls it sysmon), example writes 0xA to the XADC wizard base address.  Self test consists of writing and reading back the Vccint alarm threshold--over memory mapped 32-bit register access.  To configure the sequencer, write safe (default) sequencer mode to the CFR1 (configuration register 1) bits [15:12] @ offset 0x304 from the XADC wizard base address.  This is apparently necessary before enabling channels (by writing to the group of 8 registers @ offset 0x320 ~ 0x33C).  The format of these channel sequence registers are NOT given in PG091 (XADC wizard IP documentation), but rather in UG480 (XADC HW documentation), Chapter 4, Automatic Channel Sequencer.  To add even more confusion, sequence 8 and 9 configuration registers are NOT even continuous with the rest, but in what PG091 calls "test register" block: @ offset 0x318 and 0x31C, and therefore the bits of the register are not even documented.

XADC clock rate is configurable by changing the clock divisor in XADC wizard configuration register 2, @ offset 0x308 (documented in UG480, table 3-6).

Config register @ offset 0x300 MUX bit must be asserted--with an associated input channel--to use the external mux feature.

The BSP generated example only checks for EOS (end of sequence), and reads out all channels.  When reading the converted value, one generally reads from blocks at offset 0x200 through 0x2B8.  What is NOT obvious from the documentation that is clarified in the example code is that while the XADC input is shared for different channels, the conversion results are stored in the respective readout registers according to the channel.  Curiously, the example does NOT right shift the 32 bit readout value, even though the documentation says that the data is MSB justified.

I/O

As you can see in this image, the XADC can become an IO pin hog.  If I wanted to sample 16 channels, I would have to sacrifice 32 I/O pins (UG480: "all analog input channels are differential and require two package balls...  The analog inputs of the ADC use a differential sampling scheme to reduce the effects of common-mode noise signals")!
Vp/Vn pins should be grounded if not used.  Normally, the Vivado XADC wizard will only expose used pins, according to the channel selection mode.  But it seems that the XADC wizard always instantiates the Vp/Vn pins.

A way to reduce the I/O pin count is to use an external multiplexer.  For example, when in simultaneous mode, two 3-bit multiplexer chip (something I have to solder outside Zynq chip) can select among the 8 pairs of input channels, as you can see below:
XADC wizard configuration does not seem to expose this feature, because I cannot select a pair of inputs into the XADC wizard module...

Zedboard XADC

That the Zedboard is almost the "Cadillac of the Zynq eval boards" can be seen when you see the external reference voltage available for ADC, and that the analog ground is decoupled from the digital ground with a ferrite bead.

Reference voltage

ADCs need reference voltage; sometimes for the external circuit, but often for internal digitization implementation as well.  XADC reference voltage requires 1.25 V, which should be placed as close as possible to the reference pins and connected directly to the V REFP input, using the decoupling capacitors recommended in the reference IC data sheet.

Zedboard is generating the 1.25V Vref, as you can see here:
This Vref is available for circuits that need high quality reference voltage (like thermistors) on the Zedboard's XADC header (pin 11):
If I want to use some other reference voltage than the Zedboard generated one above, I can just lift R186 below, and connect pin 11 above to the desired reference voltage.
For XADC on Zynq, Vref pin should be tied to ground if internal reference is to be used.

The DXP/N (XADC-DXP/N: pins 7/12; according to Zynq TRM Table 2-13 PL Pin Summary, these are temperature sensing diode pins) are bit of an anachronism; according to a Xilinx forum discussion, they should not be used any more.

XADC input header on Zedboard

As shown in the XADC header pin assignment above, there are 3 differential pairs of ADC input accessible on the XADC header Zedboard (DXP/N are completely independent of XADC):
  1. XADC-VP/N: dedicated analog input pair on pins 2/1.  There is only 1 such dedicated analog input on the XADC.  When not used, these should be shorted to ground.
  2. XADC-AD0_P/N: pins 3/6
  3. XADC-AD8_P/N: pins 8/7
Note that the pairs 0 and 8 can be SIMULTANEOUSLY sampled if using the simultaneous selection mode in the XADC wizard.  These are are low-pass filtered on Zedboard with R=100 Ohm and C=1nF as shown below, so that the filter bandwidth is 10 MHz.
This is hardly sufficient for low frequency signals like those from an accelerometer, for which 100 Hz BW is more appropriate and can be accomplished by using a 10 kOhm resistor and a 1 uF capacitor combination.  This is easily possible with an extra 1 uF capacitor across the P/N pins, and a 1 kOhm resistor on the positive pin.

An FMC connector breakout board allows access to all other differential input pairs, all of them on bank 35 (VADJ).

Analog ground

When I looked at the accelerometer signal on a scope, I connected the analog ground to the Zedboard's digital ground.  It is NOT advisable to use the digital ground as an analog ground reference for XADC, because the ground "shakes" in response to the high frequency digital currents.  In an effort to improve the ADC performance, a dedicated supply and ground reference is provided on the Zedboard: A ferrite bead filters out the high frequency noise from the command ground is fed to the Zynq's analog ground input pin GNDADC_0, as you can see below.
The jumper J12 lets you bypass this ferrite bead, although I cannot imagine why you would want to do that.  It is this ground that I want to use for the analog ground of the accelerometer, so I should short JP12's pin 1 to the common analog ground output of the accelerometer AND the negative pins of the XADC's differential input pins XADC-VN, XADC-AD0_N, XADC-AD8_N in the XADC header.

I cannot directly connect the accelerometer output's 3.3V to the XADC positive pins, because the maximum voltage difference between the XADC P/N pair is 1 V; I have to use a voltage divider instead; something like 2.5 kOhm and 1 kOhm should work.

Including XADC IP in a Zynq HW design

While Zedboard itself is ready for XADC, the XADC is an optional IP in Zynq; the "Ubuntu on Zedboard" Zynq reference HW design I have been using for the bare metal, hard-real-time SW development till now LACKS the XADC block.  To include it, open the system.bd (block diagram) in a Vivado project--> click the "Add IP" icon (the one with the + symbol) on the left toolbar -->  Enter "xadc" in the search window --> double-click on XADC Wizard --> double-click on the resulting IP block to configure it.  To understand the configuration options for the XADC Wizard, read Xilinx document PG091, rather than UG480 which explains the HW level primitive that the wizard encapsulates.  The following setup targets human-machine interface application that uses accelerometer and mic, and 1 thermistor for ambient temperature sensing.
  • Basic tab
    • Leave the interface as AXI4Lite, to control XADC from the SW through memory mapped registers
    • Startup channel: sequencer, with MUX, with Vp/n and all Vaux pairs enabled (except for vaxu4 pair, which cannot be selected for some reason, as shown below)
      • Single channel mode is inefficient if many channels need to be sampled, because setting up a channel, and then waiting for the conversion takes many clock cycles (26 ADCCLKs for signal acquisition and conversion in a continuous mode, but commanding the XADC over the shared and slow AXI4Lite bus will take many more cycles).  The biggest problem with the single channel mode is that a channel is selected (fixed) in the wizard, so there is no chance to monitor multiple channels without using an external mux.
      • In all products I worked on, there are MANY ADC channels to monitor, so either a sequencer or simultaneous (sequencing) mode are used.
        • In simultaneous mode, 8 pairs channels can be converted simultaneously.  The pairs are hard coded in HW to 0/8, 1/9, ..., 7/15.  Note that Vp/Vn pair cannot be sampled in this mode.  In XADC wizard IP configuration window, selecting simultaneous selection does NOT add another mux input, adding to the confusion.
        • Independent ADC is necessary if you need to monitor the Zynq die voltage/temperature constantly, even while monitoring external channels.
    • Temp bus if for monitoring the DRAM internal temperature, which I am going to ignore for now.
    • Timing mode: event mode vs. continuous; since temperature does not change that fast, there is no need to continuously sample the channels at maximum throughput.  Actually, I've used continuous sampling in the past mostly for convenience.
      • Event mode trigger: convst_in.  Since the IP acts on the OR of the convst_in and the CONVST register (which I will control), I conected convst_in to a Const IP module to pull it up. 
  • ADC setup
    • I enabled external mux even though I only want to sample 1 channel (external thermistor) for now, because a real application will have many channels. I will use the dedicated analog input channel Vp/n.
  • Alarms tab: turn off all temperature and voltage alarms, to reduce the pin count of the core. 
    • Note: a simple control application implements its own alarm in SW, but some safety critical application may want a SW independent alarm trip (that can inhibit the circuit that is driving significant current/voltage).
When using the mux, the SW cannot control which channel to sequence: the XADC picks the next MUX channel, as you can see from this screenshot of UG480:
My understanding is that the core remembers the current channel, and advances to the next one.  XADC can (in simultaneous sampling mode) sample 2 simultaneous inputs, with 2 muxes.

The resulting XADC interface, with only the pins I connected (see explanation below) are shown:
Since the interrupt goes high for any of the interrupt conditions that can be enabled through the interrupt enable register, the eoc_out, alarm_out, eos_out seem redundant.  But because I am only interested in the eoc for now (or eos, if I use the channel sequencer in the future), I connect eoc_out to the only free pin left in the interrupt concatenator (sys_concat_intc[2]).  As explained in a previous blog, the Zynq interrupt ID for this interrupt line would be 89.

The existing "Ubuntu on Zedboard" already contains an AXI interconnect for the AXI Lite interfaces.  I keep adding more IPs to this multiplexer, as shown in a previous blog.  I just add another AXI4Lite port, and connect that to the s_axi_lite interface shown above.  As before, I let Vivado assign register address for this new module, and obtain 0x43C10000 as the base address.

Although I will only sample 1 channel in this investigation and therefore do NOT need a MUX, I will bring out the MUX to the top level for posterity's sake.  muxaddr_out should be used with an external 5 bit mux, to switch between channels.  While this approach saves pins, this presents a problem to a multi-axis device such as an accelerometer, whose 3 axis should ideally be sampled simultaneously.  Even with the simultaneous sequencing mode, XADC cannot sample more than 2 channels at the same time.  To expose muxaddr_out bus to the top level, right click on the bus in the system diagram view --> make external --> rename to XADC_mux.  Similarly, make the ADC input channel Vp_Vn external (and optionally, rename the interface to XADC_in).

Conversion can be started by driving const_in high, but I left it disconnected because I will start the conversion from the SW.

After generating the HDL wrapper (right click on the block design --> Sources --> Hierarchy view --> system.bd --> Generate HDL wrapper), these are the new interface wires to the system:

  input XADC_in_v_n;
  input XADC_in_v_p;
  output [4:0]XADC_mux;

These wires should be brought out to the top level Verilog file verbatim, so the block diagram's pins can be constrained.  Therefore, the top level module will just copy these wires verbatim as well.  Finally, they have to be constrained.  This is what PG091 said about constraining the XADC input pins: "VP/VN and 16 VAUXP/VAUXN pin pairs do not need LOC constraints to be specified in XDC.  VP/VN is a dedicated input and VAUXP/VAUXN I/Os are dual mode I/O for 7 series FPGAs.
Vivado tool performs placement of these analog inputs automatically, but VP/VN and 16 VAUXP/VAUXN pin pairs need analog I/O standard constraints for the implementation."  These are the Zynq XADC VP/N pins I read off the Zedboard schematic:
  • VP: L11
  • VN: M12
Since these reserved pins are in bank 0, where all other pins seems to be connected to 3.3 V, I guessed at the IO standard.  Contrary to the documentation, I found that I had to locate the Vp/n pins for implementation (DRC: design rule check) to pass.

set_property  -dict {PACKAGE_PIN  L11   IOSTANDARD LVCMOS33} [get_ports XADC_in_v_p];
set_property  -dict {PACKAGE_PIN  M12   IOSTANDARD LVCMOS33} [get_ports XADC_in_v_n];

The mux output can be placed anywhere, so I chose the unused OLED pins.
  • XADC_mux[0] --> V20
  • XADC_mux[1] --> U20
  • XADC_mux[2] --> V19
  • XADC_mux[3] --> V18
  • XADC_mux[4] --> AB22
VGA R, G, B channels mix multiple bits into analog value as you can see below, so the resistors will have to be removed if the mux will actually be used in the future.

set_property  -dict {PACKAGE_PIN  V20   IOSTANDARD LVCMOS33} [get_ports XADC_mux[0]];
set_property  -dict {PACKAGE_PIN  U20   IOSTANDARD LVCMOS33} [get_ports XADC_mux[1]];
set_property  -dict {PACKAGE_PIN  V19   IOSTANDARD LVCMOS33} [get_ports XADC_mux[2]];
set_property  -dict {PACKAGE_PIN  V18   IOSTANDARD LVCMOS33} [get_ports XADC_mux[3]];
set_property  -dict {PACKAGE_PIN  AB22  IOSTANDARD LVCMOS33} [get_ports XADC_mux[4]];

After a successful bitstream generation and export (of the hardware definition and the bitstream), the memory mapped registers for xadc_wiz_0 appears at 0x43C10000, which the SW can then access.

A thermistor circuit driven by 1.25 Vref

For this experiment, I will use a 10 KOhm thermistor (this is the resistance at 20 C), which has an inverse relationship between the measured temperature and the resistance.  The input to the XADC Vp should be the high voltage side of the thermistor and the Vn should be the analog ground, in a voltage divider configuration, like this:
The 5 KOhm value shown above is just an example.  To properly size the voltage divider resistor, we have to know the expected operating temperature, and more importantly, the nominal resistance of the thermistor at those 2 extremes.  The Vp-Vn into the XADC will be the greatest when the thermistor resistance is high--that is, when the temperature at themistor is low.  Let's say we want to measure between 0 C to 50 C.  An NTC themistor I have the following resistance values:
  • @ 0 C,  32336 Ohm
  • @ 50 C, 3635 Ohm
  • @ 100 C, 700 Ohm
 If I want Vp-Vn to be 1.0 V at 0 C, then I have to solve for the voltage divider equation: 32336/(32336 + R) = 1.0/1.25 => 32336 * 1.25 = 32336 + R => R = 32336 * (0.25) = 8084.  Using a 10 kOhm resistor, these are the Vp-Vn values I would see:
  • @ 0 C, Vref * 32336 / (32336 + 10k) = 0.954 V
  • @ 50 C, Vref * 3633 / (3633 + 10k) = 0.333 V
The ln(Rthermistor) is NON-linear, and the slope flattens out at higher temperature; on the cold side.  At the highest temperature, the thermistor resistance is roughly 3 kOhm, so putting a 1 uF capacitor across Vp-Vn would give me a low-pass filter with bandwidth of 1/(3k * 1E-6) ~ 330 Hz, which would reject a lot of high frequency noise.  1 uF capacitor would do an even better job of high frequency noise rejection.
The low-pass resistor 100 kOhm introduces a bias term (on the order of 1/10th) during transient.  A larger resistor (and a smaller capacitor) would reduce the error.

Bare metal C++ to read XADC wizard

Most ADC examples you will find on the web uses continuous sequencing, which is easier to program than "sequence only when needed" but wastes power.  In this code, I kick off a channel sequencing (sampling multiple channels in sequential order) in a system tick timer.  Since the XADC HW should be completely deterministic, the data should be available at a deterministic time WRT when it is needed (extra delay between ADC sampling and control value computation is a seldom discussed but a subtle contributor to a control loop instability).

The XADC wizard HW interface is completely over the AX4Lite shared bus, shown above.  Using the base register address assigned by Vivado, a generic macro to access the 32-bit register is:

#define XADC_WIZARD(offset) (*(volatile uint32_t*)(0x43C10000 + offset))

POST (power-on-self-test)

It is a good practice to perform some POST in FW, by explicitly resetting the HW to the default state and performing a basic "wiggle toe" test.

    XADC_WIZARD(0) = 0xA;//reset is active for 16 clock cycles; see PG091, SRR
    //POST the XADC wizard
    XADC_WIZARD(0x340) = 0x55;//This is the test value the BSP example used
    Q_ASSERT(0x55 == XADC_WIZARD(0x340));
    XADC_WIZARD(0) = 0xA;//reset is active for 16 clock cycles; see PG091, SRR

In case I am writing to the wrong register that does NOT remember the values I wrote, the assert will park the FW in an error state, so that it will be obvious to all.

Initialize the XADC for 16 channels

To configure the sequencer, it must first be put into the default (safe) mode:

#define XADC_SEQUENCE_SAFE      (0 << 12 | 0xFFF)
#define XADC_SEQUENCE_OFF       (3 << 12 | 0xFFF)
#define XADC_SEQUENCE_CONTINOUS (2 << 12 | 0xFFF)
#define XADC_SEQUENCE_1PASS     (1 << 12 | 0xFFF)

    //Disable channel sequencer before configuring the sequencer
    XADC_WIZARD(0x304) = XADC_SEQUENCE_SAFE;

Select the 16 channels.  Given the Zynq's power and flexibility, I hope to read nearly all of the 16 channels.

    XADC_WIZARD(0x320) = 1 << 11;//Vp/Vn pair
    XADC_WIZARD(0x324) = ~(1 << 4);//Enable all Vaux pairs, except Vaux4
    //Leave averaging, input-mode (bi/unipolar), and acquisition time at default

Then enable interrupts:

    XADC_WIZARD(0x60) = XADC_WIZARD(0x60);//Clear any pending interrupts
    XADC_WIZARD(0x5C) = 1 << 31;//PG091, GIER (global interrupt enable) register
    XADC_WIZARD(0x68) = 1 << 4;//PG091, IPIER register, bit EOS

Note that I am only interested in the end of sequencing (vs. end of conversion for each channel).

Enable external mux, which is necessary to save I/O pin counts:

    //Enable external mux and connect to Vp/n
    XADC_WIZARD(0x300) = 0 << 12 // no averaging
    | 1 << 11 //enable mux
    | 0 << 9 //event driven sampling: doesn't work?
    | 3 << 0; // mux input channel = Vp/n

Finally, put the XADC in low power mode until needed:

    //Turn off all alarms and enable calibration; see UG480 Table 3-9
    XADC_WIZARD(0x304) = XADC_SEQUENCE_OFF;//Write to the config register 1

    XADC_WIZARD(0x308) = 3 << 4; //power down ADC to save power

Start ADC acquisition and conversion from a system tick handler

My system tick timer handles timeout driven processing, so I want that to be jitter-free.  Even though starting the sequence itself is deterministic and can potentially be done before I handle the timeout event in the FW, the timeout handling code for all my state machines may take more than the sequencing (of the 16 channels).  If so, the timeout handler may get interrupted by the completion of the ADC sequencing.  While I write the code to be re-entrant in general, avoiding multi-thread contention is a good practice in general.  If sampling ADC and then running a control algorithm is deemed higher priority, I can easily move the following code to BEFORE my timeout handler.

//process the system tick (Q_TIMEOUT_SIG)

//Start the ADC conversion AFTER systick processing is done, to
//keep the Q_TIMEOUT_SIG jitter-free
XADC_WIZARD(0x308) = 0;//Start the XADC clocks

//Start another pass through sequence; PG480
XADC_WIZARD(0x304) = XADC_SEQUENCE_1PASS;

It takes MANY clock cycles to acquire and convert the 16 channels, as you can see in the scope capture below, where the 1st pulse is the timeout interrupt handling (at the end of which I start the acquisition, as explained above):
Eyeballing from the scope capture timeline, it takes about 19 usec for acquisition and conversion to complete, and another 5 usec to read out the ADC values from the registers (code given below).

Handle the EOS (completion of ADC sequencing) interrupt

Please remember from the above discussion that in my HW design that uses an FPGA to CPU interrupt concatenator, the XADC interrupt was mapped to interrupt ID 89, which I enum to INT_ID_XADC.  The first thing I should do is acknowledge the interrupt (turn it off)

case INT_ID_XADC: {
uint32_t status = XADC_WIZARD(0x60);//PG091, IPISR register
XADC_WIZARD(0x60) = status;//acknowledge interrupt
XADC_WIZARD(0x304) = XADC_SEQUENCE_OFF;

If the interrupt is EOS, I can read out the data.

if(status & (1 << 4)) {//EOS (end of sequence)
XADC_val[0] = XADC_WIZARD(0x20C);//PG091 table 2-4, Vp/Vn
XADC_val[1] = XADC_WIZARD(0x240);
XADC_val[2] = XADC_WIZARD(0x244);
XADC_val[3] = XADC_WIZARD(0x248);
XADC_val[4] = XADC_WIZARD(0x24C);
XADC_val[5] = XADC_WIZARD(0x254);
XADC_val[6] = XADC_WIZARD(0x258);
XADC_val[7] = XADC_WIZARD(0x25C);
XADC_val[8] = XADC_WIZARD(0x260);
XADC_val[9] = XADC_WIZARD(0x264);
XADC_val[10] = XADC_WIZARD(0x268);
XADC_val[11] = XADC_WIZARD(0x26C);
XADC_val[12] = XADC_WIZARD(0x270);
XADC_val[13] = XADC_WIZARD(0x274);
XADC_val[14] = XADC_WIZARD(0x278);
XADC_val[15] = XADC_WIZARD(0x27C);

XADC_WIZARD(0x304) = XADC_SEQUENCE_OFF;
XADC_WIZARD(0x308) = 3 << 4;//Stop the XADC clocks
}

It was not obvious at first why the sequencer could not be left in a single pass mode, until I found a cryptic passage in the datasheet while debugging: the sequencing starts when the mode CHANGES to single pass mode!  Stopping the clock is the PRINCIPLE method of saving power, but I have to wait until the sequencing is complete.

Result without any circuit connected to the Vp/n input (just reading noise)

This is what I see in the JTAG debugger:

XADC_val long unsigned int [16] [0x000020b7, 0x000020d2, 0x000020c3, 0x000020d1, 0x000020c7, 0x000020d0, 0x000020d5, 0x000020d1, 0x000020bd, 0x000020ca, 0x000020b5, 0x000020c6, 0x000020d7, 0x000020ca, 0x000020bd, 0x000020c9]

Since the XADC is a 12-bit ADC, the leading '2' is strange, and probably an invalid bit that should be thrown away.  The values change with every read, so I think the read values are indeed noise.

Note that the Xilinx documentation is WRONG: it said the data is MSB justified.  Clearly, the data is LSB justified!

Removing the XADC Linux device driver from the kernel

The  XADC device driver can be removed from the DTS file, but that will leave unused in the kernel.  For a thorough excision, I changed the kernel config like this:

# CONFIG_XILINX_XADC is not set