An Illustration of Central Limit Theorem Using R

Feb 4th, 2018 | Comments

The Central Limit Theorem is one of those spooky math results that probably keep mathematicians awake on full moon nights. This post will try and illustrate the central results of the theorem so that some of us slumbering mortals might also lose some sleep occasionally.

In essence the theorem states that given a bunch of numbers (aka the Population), if you randomly select small - but big enough - groups of numbers (aka sample) from it and add the numbers in these samples together; the group of numbers produced as a result of this addition will always be distributed in a specific pattern called the normal distribution. The addition component is very important, it’s not enough that you select samples from the population, you must add the data within the samples for the CLT to hold. This also means that whenever you are dealing with some data that could possibly be expressed as the result of randomly summing other quantities you can apply the CLT.

The cool part here is that even though you don’t know anything about the original population, you can find it’s mean and variance if you take a set of large enough samples.

Now if instead of just adding the numbers inside the samples, we calculate their average or mean, some more interesting results can be obtained. Lets also introduce some notation to make the textual statements clearer.

Population mean = $ \mu_P $, variance = $ \sigma_P $ (aka sample distribution)

Set of sample means = $\{ \mu_0 , \mu_1 , ... , \mu_N \}$ (aka sampling distribution)

Mean of set of sample means = $\mu_X$, variance = $ \sigma_X$

The mean ($ \mu_X $) of the set of sample means will approach the actual mean of the population
The variance ($ \sigma_X $) of this mean is approximately equal to the population variance divided by N ($ \sigma_P/N $), where N is the size of each sample (i.e. the count of numbers in each sample)

More about N

The sample size N is rather important, it must be reasonably big for the CLT to hold. If it’s too small, then there’s not going to be any soup for you. The right size of N depends on the original sample population. If the original population is a decent symmetric one then - according to the math gurus - a value of N <= 30 is good enough. However if the original population is skewed, then you might need higher values of N. For example imagine an impulse function - a value of K at one point and 0 elsewhere. You might need a large sample from the population to end up with a sample that has K inside it.

Proof or maybe not

Since the proof of a pudding is in the coding and not in the deriving, let try some experiments.

The below script calculates the Population and Sample statistics based on CLT, and plots a histogram of the sample means for a given population function. Feel free to play around with different settings for N and P. I use N = 40 as the threshold in the script, but 30 should also yield similar results.

[Central Limit Theorem] [lang : r] (central_limit.r) download

# Central Limit Theorem Illustration

# THEOREM
# ==========
# Tells us that the sampling distribution of the sample mean is, at least approximately, 
# normally distributed, regardless of the distribution of the underlying random sample.
# The maean of this distribution of samples means is EQUAL to the population mean
# The variance of this (new) mean of distributions is EQUAL to the (population variance / N)
# [The error with respect to the population mean decreases as N increases]
# where N is the sample size
# NOTE:
# (1) The larger the sample size N, the smaller the variance of the sample mean.
# (2) If the underlying distribution is skewed, then you need a larger sample size, typically N > 30
#     for the results to hold
# (3) This is only applicable to the Mean, variances follow the Chi-squared distribution
#
# https://onlinecourses.science.psu.edu/stat414/node/177
# EXTRA:
# The sampling distribution of the sample variance is a chi-squared distribution with degree of 
# freedom equals to N−1, where N is the sample size
# The chi-squared distribution looks skewed compared to normal distribution
# http://onlinestatbook.com/2/chi_square/distribution.html
#

# population size
P = 50000
# sample size --> Increase this to getter closer to ideal results
N = 40
# sample count --> This is also important in practise
C = 1000



## ----- Setup the population ----
# Select any probability function and convert it to discrete form with P points

# (1) Exponential function - change the factor from 1 - 100 to increase skew
# Sequential sampling of the PD function
pd = exp(seq(from = 0.0, to = 10.0, by = 10.0/P))
# Random sampling of the PD function
#pd = exp(runif(P, 0.0, 10.0))

# (2) Ramp function
#pd = seq(from = 0.0, to = 10.0, by = 10.0/P)

# (3) Sine function
pd = sin(seq(from = 0.0, to = 2*pi, by = 2*pi/P))

# (3) Random function
pd = runif(P, 0, 10^12)

# Allocate a matrix to hold the sample data
s = matrix(0,C,N)

## ----- Take C samples of same size from the population ----
for (i in 1:C) {
  # Create a matrix of random indices with sample size
  idx = matrix(sample(c(1:P), size = N, replace = TRUE))
  # Select random samples from our probability function
  s[i,] = pd[idx]
}

# Now s[] contains C samples from the population pd
# Lets do some statistics on the samples

par(mfrow=c(3,1))
plot(pd,main="Population function")
## MEAN
# prob TRUE converts the counts into probability and is necessary to draw the density function overlay
hist(rowMeans(s), plot = TRUE, main = sprintf("Histogram of sample means (N=%i)", N), prob=TRUE)
lines(density(rowMeans(s)), col="blue")

## Lets calculate the variance of all the samples
v_s = apply(s, 1, var)
hist(v_s, plot = TRUE, main = "Histogram of sample variances", prob=TRUE)
lines(density(v_s), col="red")

print("==========================\n")
buf = sprintf("Population:\t Mean = %f, Variance = %f, Variance/N = %f\n", mean(pd), var(pd), var(pd)/(N))
cat(buf)
buf = sprintf("Sample Means:\t Mean = %f, Variance = %f\n", mean(rowMeans(s)), var(rowMeans(s)))
cat(buf)
print("==========================\n")

A sinusoidal population

The population numbers are extracted from one period of the sin() function. Running the script produces the following result

Population:	 Mean = -0.000000, Variance = 0.500000, Variance/N = 0.012500
Sample Means:	 Mean = -0.002520, Variance = 0.012633

Since the function is nice and symmetric, a value of N=40 is more than enough to result in a sample mean that’s almost equal to the population mean, and the sample mean variance is also almost equal to that predicted by the (known) population variance.

An exponential population

This population is extracted from the exp() function.

Population:	 Mean = 2202.722807, Variance = 19411026.891557, Variance/N = 485275.672289
Sample Means:	 Mean = 2178.442988, Variance = 467492.546106

The results are not as perfect as before, but it’s still pretty close. Increasing the parameter of the exponential function from 10 to 100 will result in a plot that’s clearly not normal. Increasing N to 200 brings us closer to CLT expectations.

A random function

To illustrate that the CLT does not depend on the original population distribution, we set the population to be a random function. One can see that for a random data, N = 40 is pretty good.

Other considerations

Sampling with and without replacement

The examples are all based on samples taken with replacement, meaning that every sample is drawn from a population of P = 50000, and the samples taken are not deleted from the population. If you reduce the size of the population as you take samples, then the CLT formula needs a correction in the variance formula as specified here. However this can be ignored if the population size is much larger (20x) of the sample size.

Random data

It’s important that you take random data for each of the samples, otherwise you won’t get the desired results. For instance imagine that you keep on repeating the same sample of data, in this case you may never reach close to the population mean until N becomes almost as large as the population. In practice this is a difficult requirement to get right because you might unknowingly introduce biases. For example if you do an online survey as a means of sampling, you might only be targeting a certain segment (eg: tech savvy people) of the population which might invalidate your results with respect to other segments of the population.

Conclusion

That concludes the introduction to CLT. A lot of statistics tests introduce CLT using probability concepts, however CLT - as you can see - is a property of the sums of random numbers and does not require a background of probability.

All the plots I have generated show sampling distribution of variances (in addition to the mean) - this does not follow the normal distribution but instead follows the chi-squared distribution with N degrees of freedom, though to be honest, I can’t visually tell a huge difference.

References

Fun With Ecology - Richness and Diversity

Jan 31st, 2018 | Comments

Recently I had the opportunity to play around with some ecological data for statistical analysis using the R programming language. As I didn’t know anything about ecology (or R), I got to learn a bunch of new stuff, and realized that it’s actually very interesting. This post will talk about what I learnt and also serve as an introductory tutorial for someone else getting into stuff like this. However, do note that I’m not an ecologist or a data scientist, so the usual disclaimers apply - if you discover any errors please point them out.

This is the first article and will talk about the basics and species richness and diversity. Further articles will deal with Multivariate and Univariate analysis.

Software of the trade

A fair amount of proprietary software packages for ecological data analysis seem to be available, with PC-ORD, Canoco and PRIMER supposedly being the most popular. A larger list is available here. Luckily for us there are several FOC R packages that can do most of what the proprietary packages do, with greater flexibility. In the R universe, the most important package is vegan.

There is another wrapper package called BiodiversityR that wraps vegan and other stuff to present a GUI interface that allows you to quickly do any fancy test you want. It also comes with an excellent user manual that’s accessible to muggles. Given this, the best workflow is to experiment with the BiodiversityRGUI() and then eventually write your own scripts that natively access vegan APIs for better control over the plotting and parameters for the final publication worthy results.

With regards to species diversity and richness estimation SpadeR is another useful package. It comes with an online version that users who don’t know how to use R can use. EstimateS is standard proprietary software specializing in diversity and richness estimation.

Inputs

Samples and the Species / Environmental Matrix

Generally after you’ve designed your experiments and collected data in the forest or sewer, identified the species you discovered and measured the environmental conditions you’ll be left with two tables of data. The most important one is called the species matrix and this tells you what species you discovered, in what quantity (aka abundance) every time you took a “measurement”. Each “measurement” is called a sample in literature, and typically you’re going to have samples separated across location and/or time. The term sample was a bit confusing to me as coming from a signal processing background I thought it meant each individual species sample. However it doesn’t, and ecologically speaking if you go down to your neighbor’s pond today and come back with a fistful of fish, then that’s one sample. If you do that again tomorrow, that’s another sample.

The species matrix is by convention a series of species samples, with each sample assigned to a new row. The columns are the abundance of each species. All software programs expect such an arrangement by default.

Now in addition to the abundance information you are probably also interested in knowing about how other “factors” (or independent variables) effect the abundance or existence of various species. These “factors” could be environmental (like temperature, pressure..) or not (like location, your collection method..), quantitative (anything that can be assigned a numeric value) or not (like the name of the gardener at each site). These are all recorded in the environmental matrix, and you need to have an entry for each sample in the corresponding species matrix.

Example Analysis

We use some contrived data to illustrate the concepts.

The Accumulation table is:

A	B	C	D	E	F	H	I	J	K
91	1	1	1	1	1	1	1	1	1
90	1	1	1	1	1	1	1	1	0
10	10	10	10	10	10	10	10	10	10
20	20	20	20	20	20	20	20	20	20
1	20	11	28	0	0	20	20	0	0
0	0	11	28	0	0	20	20	20	1
1	1	1	1	1	1	1	1	1	1
1	1	1	1	1	1	1	1	1	1

The corresponding “environment” table:

Type
Uneven
Uneven
Even
Even
Random
Random
Low
Low

The Uneven set of samples show a skew with species A being dominant and the rest of the species even at 1 individual each. The Even dataset has an even abundance, the Random set has varying abundance and the Low dataset has the smallest total number of individuals.

The following commands load the test abundance and environment data

# Load the BiodiversityR UI
library(BiodiversityR)
BiodiversityRGUI()

# Load data
spa.dummy = read.csv('C:/Users/gov/Google Drive/Zoology/species-dummy.csv', header=TRUE)
env.dummy = read.csv('C:/Users/gov/Google Drive/Zoology/env-dummy.csv', header=TRUE)

Species Richness and Diversity

Species richness is a scalar number to illustrate the variety of species available in an area. Introducing this measure immediately creates a problematic question, is a sample with 9 species of 1 individual each + 1 species of 100 individuals equally rich as one with 10 species of 10 individuals each? The solution is to answer yes to this question, and introduce another measure called diversity that takes into account the relative abundance in addition to the richness.

Richness

When calculating richness with small samples of data, you run into two main problems

There is a good chance that you haven’t sampled enough to detect all the present species
There is likely to be a wide variation in number of individuals you have in the different samples

The second problem can be solved by using rarefaction curves.

Rarefaction

Rarefaction is an easy graphical approach which involves plotting the species accumulation curve backwards. To plot the curve you take the pooled data and randomly remove an individual / sample and mark the resulting species richness; the process is repeated until you have 0 individuals or samples.

As BiodiversityR does not seem to have a GUI option to plot individual based rarefaction instead of site/sample based, I’m using the vegan commands directly as so:

# Pool the two samples per community
spa.dummy.pooled = matrix(data = c(colSums(spa.dummy[1:2,]), colSums(spa.dummy[3:4,]), colSums(spa.dummy[5:6,]), colSums(spa.dummy[7:8,])), nrow = 4, byrow = TRUE)

# Calculate lowest number of individuals towards which we must rarefy
raremax = min(rowSums(spa.dummy.pooled))

# Plot
rarecurve(spa.dummy.pooled, step = 20, sample = raremax, col = c("red", "blue", "green", "orange"), xlab = "Individuals", cex = 0.6)
legend("bottomright", legend = c("Uneven", "Even", "Random", "Low"),
	   lty=c(1,1,1,1), # gives the legend appropriate symbols (lines)
	   lwd=c(2.5,2.5, 2.5, 2.5),
	   col=c("red", "blue", "green","orange"), # gives the legend lines the correct color and width
	   cex=0.75 # reduce size of legend box
	   )

To compare communities, you check the species richness for all communities at the point corresponding to the lowest abundance (the vertical line in the plot below). The horizontal black lines correspond to the respective rarefied richness.

Non-parametric richness estimators

Although rarefaction is simple, it’s accuracy is not so obvious. The curve does not predict the eventual species richness, and it’s possible that with more samples the curves may intersect and change in rank; and sampling down means that the lowest sample quantity should be good enough for an accurate deduction. (Apparently this is possible with EstimateS, but I wouldn’t know)

Another class of methods rely on estimating the eventual species richness based on the data you have. Many of them are based on a technique called Good-Turing frequency estimation, which correct the observed richness with estimates of undetected species calculated by looking at the number of singletons/doubletons and low abundance species. The idea being that the existence of singletons etc allude to the fact that you don’t have enough samples and there’s a good chance of discovering more species, where as if you have large abundances for all the species concerned, then your existing data is probably good enough.

There are a bunch of non-parametric methods available called Chao1, Chao2, ACE. All of these can be tried out online using SpadeR and in BiodiversityR. In vegan the specpool() function helps with this.

Calculating the Chao diversity index gives the following result, showing that we expect some new species to be discovered at the Uneven plot while the richness for the Even plot is as calculated from the data without any extra correction.

Diversity.1 <- diversitycomp(spa.dummy, y=env.dummy, factor1='Type', index='chao' , method='pooled', sortit=FALSE, digits=6)
Diversity.1

Type     n    chao
  Even   2 10.0000
  Low    2 10.0000
  Random 2  9.0000
  Uneven 2 10.0278

Diversity

Next we turn to diversity which takes into account the relative distribution of species (aka evenness) along with the total number of species. The diversity calculated from our data is the diversity at individual spots/samples and is called alpha diversity.

A popular measure of diversity is the Shannon Entropy function invented by Claude Shannon in his Mathematical Theory of Communication. This function measures entropy which is the uncertainty of the process under consideration. If you apply this to an abundance sample, what you get is the uncertainty in species identity - i.e. if you take a random individual from the sample, how confidently / easily can you guess that it belongs to a particular species X. The function attains it’s highest value when all species are equally probable (i.e. the species distribution is perfectly even).

$\begin{align} H & = -\sum_{i=1}^N p_{i} \ln{(p_{i})} \end{align}$

When the diversity is evenly distributed, the equation reduces to

$\begin{align} H & = -\ln{(1/N)} \end{align}$

Because of the log here you have to be a bit careful when comparing communities using this index directly, for example (using the base 10 logarithm) if you have two communities with evenly distributed species (for simplicity), one with 10 evenly distributed species and the other with 100, the two communities will have an entropy of 1 and 2 respectively, even though intuitively you would expect the second community to be 10x more diverse than the first one. To fix this problem, one can transform the Shannon Index by taking the inverse logarithm of the index. For a given diversity index this is equivalent to finding the diversity of a community with the same diversity index but with a perfectly evenly distributed population. This technique is applicable for other diversity indices too as is detailed here based on the earlier work by Hill here

Therefore to summarize, it’s helpful to transform the Shannon index to “true diversity” by

$\begin{align} D & = \exp{(H)} \end{align}$

This transformed value is also called the first order Hill Number.

The Shannon Diversity for our test dataset shows that - as expected - the Even and Low dataset is the same and higher compared to the other two datasets. BiodiversityR does not seem to have an option for Hill Numbers, so you’ll have to calculate it separately.

Diversity.1 <- diversitycomp(spa.dummy, y=env.dummy, factor1='Type', index='Shannon' , method='pooled', sortit=FALSE, digits=6)
Diversity.1

Type     n  Shannon
  Even   2 2.302585
  Low    2 2.302585
  Random 2 1.756506
  Uneven 2 0.480094

The classical Shannon diversity index suffers from the same problem as that of the simple richness measures - unknown species. A similar correction factor can be introduced as per Chao et al.. You can experiment with this using SpadeR.

Rank-abundance Curves

Another alternative to the above method of calculating a single numeric number to represent the diversity of a place is to analyze diversity visually. A simple rank abundance curve can show you both the richness and evenness - curves with more species extend longer and those with better evenness are more “horizontal” in nature.

The rank abundance curve of the test data clearly shows the unevenness of the Uneven and Random data points in contrast to the Even and Low datasets.

References

An excellent overview of richness and diversity measures along with details on plotting Diversity profile is in Anne Chao’s paper here and also here
Applications to microbial diversity is illustrated here.
Simplified version of Shannon’s paper - Information theory for intelligent people
The books - Measuring Biological Diversity and Biological Diversity - frontiers in measurement and assessment by Anne E. Magurran provide an in-depth treatment of biodiversity indicators

Using GCC and Clang With Eclipse on Windows

Sep 13th, 2016 | Comments

Running GCC or Clang on Windows has never been easier thanks to the MSYS2 project. Setting up the Eclipse CDT environment to use these compilers is also very easy. This post will illustrate exactly how easy it is.

Installing MSYS2

Head over to the MSYS2 website and follow steps 1 - 7. At the end of it you should have a working MSYS2 shell.
Next you need to install GCC + other utilities (MAKE, GDB, BINUTILS, … etc) for windows. The Mingw-64 project provides versions of GCC that are capable of running natively on Windows and building Windows 64/32 bit target binaries, and MSYS2 has great support for mingw-64. Download the mingw64 packages using the MSYS2 package manager pacman.

There are a couple of ways to go about this - you can either download the individual packages (gcc, gdb …) or you can download all the relevant packages in one shot using the mingw-w64-x86_64-toolchain or mingw-w64-i686-toolchain depending on what you want to use. You can search the available packages using the command pacman -Ss string_to_search.

# This should install the toolchain into the /msys64/mingw64 folder
pacman -S mingw-w64-x86_64-toolchain
# Install make if it's not already installed. Eclipse will use make to run it's autogenerated makefiles. It goes into the /usr/bin directory
pacman -S make

Once this is done you will have a new shell in the /msys64 folder called mingw64.exe or mingw32.exe (if you installed the i686 toolchain). Runnnig this shell and typing gcc -v will give you some information on what you have installed. This shell can also be used to directly compile any makefile projects you have. Next we talk about integration with the Eclipse IDE.

GCC with Eclipse

Assuming you have already installed the Eclipse IDE for C/C++ developers. Eclipse has some built-in support for Mingw projects which depends on the presence of some environment variables like MINGW_HOME, MSYS_HOME etc. However the method I outline here is not related to this and is (imho) simpler and more robust.

In eclipse go to File | New | C Project | and select either Cross GCC or MinGW GCC project. Then in the project properties go to C/C++ Build | Environment

Select the option to “Replace native environment with specified one”. This is a really useful option that allows you to customize the environment for this build without having any dependencies on the system PATH.
Set the PATH variable to have access to the relevant gcc and make binaries. MSYS2 installs make in the /usr/bin directory and the gcc executables are in either mingw64/bin or mingw32/bin.
Define a TMP variable to point to the system TMP environment variable by using the syntax ${env_var:TMP}

Now you have configured the build environment. All that’s remaining is to setup the compiler and linker options. Go to C/C++ Build | Settings | Tool Settings | Cross GCC Compiler and set the command to gcc. Similarly set the Cross GCC Linker command to gcc. The Cross GCC Assembler command should be as.

Thats it! Select apply and build the project.

If you wan to use GDB, you need to create an Eclipse debug configuration.

Clang with Eclipse

MSYS2 provides packages for Clang also (mingw-w64-x86_64-clang and mingw-w64-i686-clang). There seem to be some dependencies with mingw, so you are better off first installing mingw first followed by clang.

Thereafter the settings are exactly the same as above except that this time instead of gcc you type clang for the compiler and linker commands. The assembler is llvm-as. lldb is not yet part of the MSYS2 package, so you can’t run the debugger yet; but you might be able to run gdb.

C Tips and Tricks[0]

Sep 13th, 2016 | Comments

These posts document some of the C concepts that I often forget or confound and some idioms that I’ve encountered. It is intended as a quick reference (for myself). Most of it is not original work, some relevant sources can be found in the inline references, others might be burried in my subconscious.

Arrays

Whenever an array name appears in an expression except as the operand of unary & or the sizeof operator, the value of the array name becomes a pointer to the first element of the array.

int a[10];
int (*p_a)[10] = &a;	// &a is a pointer to an array of int with bounds 10, i.e. p_a = &a[0];
int *p_int = a;			// a is a pointer to int, i.e. p_int = &a[0]

int b[10][20];
int (*p_b)[20] = b;		// b is a pointer to an array of 20 ints, i.e. p_b = &b[0];

Arrays as function parameters

You can declare functions with array parameters, this is exactly equivalent to using the more common pointer syntax. The type qualifiers following the [ act as if they followed the * when the parameter’s type is rewritten as a pointer type. You can also use the static keyword within the [] to specify a minimum size for the caller passed array pointer. Thus:

int foo(int a[n]);  // is the same as foo(int* a)
int f(float array[const restrict]); // can be read as int f(float *const restrict array);
void foo(int a[static 10]); //  actual argument corresponding to a must be an array with at least 10 elements

Multidimensional arrays (Variable Length)

In C90 multidimensional arrays can be used with the convinent [i][j][k] indexing syntax only if the total size of N-1 dimensions (bounds a.k.a stride) are known at compile time. You don’t need to know the size of the last dimension as the stride of this dimension is equal to the sizeof(array_element_type). C99 created VLAs to overcome this restriction. Arrays whose size is determined at runtime are called Variable Length Arrays (VLA).

// C90 multidimensional arrays (using a single dimension array)
float a[m*n*sizeof(float)];
for (i = 0; i < m; ++i)
  for (j = 0; j < n; ++j)
	a[i*n + j] = 1.0;

// C99 multidimensional indexing
float a[m][n];
int i, j;
for (i = 0; i < m; ++i)
  for (j = 0; j < n; ++j)
	a[i][j] = 1.0;

References:

Randy Meyers has a series of informative articles in the now defunct C/C++ Users Journal that discuss VLAs and other topics summarized on this page in more detail. The journal arcives are available here. Dr.Dobb’s journal (also defunct) has copies of the articles too. VLA Part 1 VLA Part 2 VLA Part 3 VLA Part 4

C99 C99 Rationale C11 Final Draft

Strings

String literals v/s arrays

A string in C is simply an array of chars terminated by a '\0'. Initialized strings can be created in two ways - one is to statially allocate them in memory and hold a pointer to address where the string is stored, and the other is to store the string as an array of chars. In the first case they are constant and it’s undefined behavior to modify them, while in the second case their contents can be modified [C99 N1256 draft 6.7.8/32 “Initialization”].

/* pChar is a pointer to the actual lccation of the string, so we require an additional 4 bytes or more to store the pointer compared to string[] */
char* pChar = "I'm a string literal, I can be modified but no guarantees on the result"; // This always goes into the data segment
char string[] = "I'm a string array, can I be modified?"; // This could go into the stack or data segment depending on where it's defined

Anonymous string literals can also be used as in printf("Another literal"), in this case the pointer to the literal is implicitly passed to printf() on the stack. This by default allocation of debugging strings into data memory eat up a lot of useful RAM. Some targets have a way to relocate strings into program memory, for example AVR has the PROGMEM keyword which the Arduino platform extended to the F(x) macro.

Another trick is to use the C99 compound literal sytax to convert the string literal into a char array and push it onto the stack. However this trick only works for short length string literals, larger literals are stored in the data segment and then coped onto the stack at runtime by your smart ass compiler. The boundary was 18 bytes on an embedded MCU I tested with.

// AVR specific redirection into flash
#define PSTR(s) ((const PROGMEM char *)(s))

// Convert a string literal into an unnamed char array
printf("%s", (char []){"This is a string array stored on the stack"});
// A macro to do this
#define F(x)	(char []){x}
// A string with extended ascii chars
char ex = "An ex char to end it \xB0";

Jagged char arrays

C does not natively support jagged arrays but you can declare an array of variable length strings like:

char* pChar2D[] = {"A long long long string", "short", "something else"};
printf(pChar2D[0]); // A long long long string

How does it work? Relook at the declreation of pChar2D- it is an array of pointers to char, which means that when you say pChar2D[0] you end up getting the address of the string literal "A long long long string". Not surprisingly, this is exactly how jagged arrays of other types are programatically dealt with usually. i.e.

int r0[] = malloc(sizeof(int) * 10);
int r1[] = malloc(sizeof(int) * 2);;
int *jagged[] = { r0, r1 };
printf("%d %d, jagged[0][9], jagged[1][1]);

If you are unsure about a C declerations, they can be looked up at cdecl.org

Designated Initializers

Initializing individual elements with arrays and structs can be done with designated initializers. Elements that are not specifically initialized are automatically initialized to 0. This is clearly better than the old method of having to initialize members individually which lends garbage values to uninitialized elements. More details here

// Initializing an array
int a[100] = {[0]=5, [99]=100};

// Struct
struct S1 {
  int i;
  float f;
  int a[2];
};

struct S1 x = {
  .f=3.1,
  .i=2,
  .a[1]=9
};

Compound literals

Just like how like languages like C# allow you to create and manipulate anonymous objects, C99 allows you to do the same thing with C types like int, struct, union, array etc. These are called compound literals. You are allowed to take the address of these constants and pass them around as function parameters. The basic syntax is (<type>){<initializer-list-expression>}.

(float){3.14};
(struct POINT){10, 10};
(char []){"Compound Literal Array"};
(const char[]){"A CONSTANT compound literal array"};

Drawing Anti-aliased Circular Points Using OpenGL/WebGL

May 5th, 2016 | Comments

This article explores points in OpenGL, and I use WebGL to illustrate the concepts. The demos may not render in Internet Explorer as I use new ECMA Script 6 syntax for strings (storing the shader code) which apparently is not supported by IE.

Drawing a Point

If you discount the WebGL boilerplate code, drawing a point is easy. Attributes to the vertex shader are a vec3 point position, vec4 color of the point and a float point size. The shader updates the OpenGL built-in variables gl_Position and gl_PointSize from this and passes the color as a varying to the fragment shader. The fragment shader simply outputs the color it got from the vertex shader. Note that the maximum Point Size for WebGL is only guaranteed to be at least equal to 1 pixel, though in most cases it’s more than one.

[Our simple vertex shader] [lang : c] (vs-points.js) download

﻿var vs_points = `
attribute vec3 vPosition;
attribute vec4 vRgbaColor;
attribute float vPointSize;
varying  vec4 color;

void
main()
{
    gl_Position =  vec4(vPosition, 1.0);
	//gl_PointSize = 10.0;
	gl_PointSize = vPointSize;
    color = vRgbaColor;
	//color = vec4(1.0,0.0,0.0,1.0);
}`;

[Our simple fragment shader] [lang : c] (fs_point_default.js) download

﻿var fs_point_default = `
precision mediump float;
varying  vec4 color;

void
main()
{
    float alpha = 1.0;
    gl_FragColor = color * (alpha);
}`;

I will ignore the boilerplate code as there’s nothing special about it. You can figure it out by just going through the points.js script. The relevant high-level code is shown below:

  var point = [0.0, 0.0, 0.0];
  var color = [1.0, 0.0, 0.0, 1.0];
  var size = [400.0];

  var canvas = document.getElementById("point-canvas");
  initGL(canvas);
  shaderPointDefault = initShaders(0);
  initBuffers(point, color, size);
  gl.clearColor(0.0, 0.0, 0.0, 1.0);
  gl.enable(gl.DEPTH_TEST);
  drawScene(shaderPointDefault);

The result of this is - a nice big square!

Drawing a disk

To get a circle out of this square we need to first understand the square better. An important parameter is gl_PointCoord. This is available at the Fragment Shader and tells the location of the current fragment inside the current point. The idea is that depending on the value of gl_PointSize, a point will be actually composed of many pixel fragments and this parameter allows us to identify each of these fragments. Using this information we can draw any figure inside a “point”.

Now gl_PointCoord has an X and Y coordinate and varies from $ [0,1] $ just like a 2D texture. Recall that the equation of a circle with origin at the center and radius 1 unit is

$\begin{align} x^2 + y^2 = 1 \end{align}$

This is the same as the dot product of the 2D vector with itself. Since GLSL provides us a dot() function, we can draw a circle in the square by coloring the region in space where the dot product of a vector with itself is $ <= 1 $:

[Extract a circle from a point] [lang : c] (fs_point_alias.js) download

var fs_point_alias = `
precision mediump float;
varying  vec4 color;

void
main()
{
    float r = 0.0, delta = 0.0, alpha = 1.0;
    vec2 cxy = 2.0 * gl_PointCoord - 1.0;
    r = dot(cxy, cxy);
    if (r > 1.0) {
        discard;
    }
    gl_FragColor = color * (alpha);
}`;

All pixels outside the radius of 1 are discarded. The result is:

Since we discard all pixels outside the radius, you can see that the edges are rough and aliasing is pretty bad. Use the CNTL + keys to zoom into the web page to see the distortion better.

Anti-aliasing the disk

Instead of abruptly discarding the pixels at the fragments where $ r > 1 $, we need to smooth out the colors at the edge pixels. One way to do this is to use the smoothstep() function provided by GLSL that allows you to interpolate the value of a step function in the region where the independent variable lies between the min and max values. In the case of a normal step function that we indirectly used before, the min and max was $ [1,1] $ respectively - anything greater than 1 is clipped to 1 and anything less than 1 is clipped to 0. To use smoothstep, we need to find an $ \epsilon $ so that we can define a range of $ [1+\epsilon, 1-\epsilon] $ overwhich the color of the fragments can be smoothed.

The GLSL fwidth() function can give us this information. It returns the maximum change in a given fragment shader variable in the neighbourhood of the current pixel (8 surrounding pixels). So fwidth(r) gives us the $ \epsilon $ and we can have a shader as:

[Our simple vertex shader] [lang : c] (fs_point_anti_alias.js) download

var fs_point_anti_alias = `
#ifdef GL_OES_standard_derivatives
#extension GL_OES_standard_derivatives : enable
#endif

precision mediump float;
varying  vec4 color;

void main()
{
    float r = 0.0, delta = 0.0, alpha = 1.0;
    vec2 cxy = 2.0 * gl_PointCoord - 1.0;
    r = dot(cxy, cxy);
#ifdef GL_OES_standard_derivatives
    delta = fwidth(r);
    alpha = 1.0 - smoothstep(1.0 - delta, 1.0 + delta, r);
#endif

gl_FragColor = color * alpha;

}`;

And now we finally have a lovely circle.

A demo to put it together

The below demo draws such point-circles of randomly varying sizes, color, alpha and position. The idea for the demo is from the Gamedueno2 book The rendering code uses callbacks registered with the browser requestAnimationFrame API insead of a timer based rendering. Use the sliders to vary the maximum size of the circles and their total number and see the impact on the fps. The point size slider is clipped to the max point size obtained by querying ALIASED_POINT_SIZE_RANGE.

Max Size:

Count:

fps:

References

Anti-aliased circle using fwidth()
Shader anti-aliasing
Using requestAnimationFrame
Perhaps the best WebGL tutorials out there.
WebGL regression test using gl_PointCoord

I2C Wrappers for the Arduino Wire Library

Feb 22nd, 2016 | Comments

Motivation

The current (v1.6.7) Arduino Wire library doesn’t seem to have a convinent API to read/write from/to I2C devices. All the pieces are available but what’s lacking is a wrapper around Wire that makes it straight forward to talk to I2C slaves. I wrote one similar to the PIC code I’ve used in the past, and it’s posted here.

Initialization

The Wire interface can be initialized by calling begin(). The default speed is 100kbps, which can be increased by calling SetClock().

  Wire.begin();
  Wire.setClock(400000L); // Set clock to 400kbps

Read

Reading from a slave typically consists of first sending the register ID you want to read out from followed by a repeated start condition [using endTransmission(false)]. Then you set the read bit and clock out the data [via requestFrom()].

void i2c_read(byte seven_bit_address, byte reg, byte* pData, byte len ) {
  byte i = 0;
  Wire.beginTransmission(seven_bit_address);
  Wire.write(reg);
  Wire.endTransmission(false);     // send a repreated start
  Wire.requestFrom(seven_bit_address, len, true); // read N bytes and send a stop bit
  while (Wire.available() && (i < len)) // Copy the already read out data into our buffer
  {
    *pData++ = Wire.read();
    i++;
  }
}

Write

Write is self explanatory.

void i2c_write(byte seven_bit_address, byte reg, byte* pData, byte len ) {

  Wire.beginTransmission(seven_bit_address);
  Wire.write(reg);
  for (int i = 0; i < len; i++) {
    Wire.write(*pData++);
  }
  Wire.endTransmission();     // send a stop
}

// An even simpler interface to write a byte
void i2c_write(byte seven_bit_address, byte reg, byte value ) {
  Wire.beginTransmission(seven_bit_address);
  Wire.write(reg);
  Wire.write(value);
  Wire.endTransmission();     // send a stop
}

Example Usage

The code below shows how to read/write some data from/to the MAX30100 Pulse Oximetry sensor [I2C address = 0x57].

  Serial.begin(230400);
  // Initialize
  Wire.begin();
  Wire.setClock(400000L);

  // Do stuff..

  // Read
  byte id;
  i2c_read(0x57, 0xFF, &id, 1);         // Part ID register = 0xFF
  Serial.println("Part ID:");
  Serial.print(id, HEX);                // expected value is 0x11

  // Write
  byte data;
  data = (byte)0xA << 4 | (byte)0xB;    // Set the LED config register to 0xAB
  i2c_write(0x57, 0x09, &data, 1);      // LED config register = 0x09 
// OR i2c_write(0x57, 0x09, data);

  // Read back and check
  i2c_read(0x57, 0x09, &data, 1);
  Serial.println("LED Config = ");
  Serial.print(data, HEX);

Real-time 2D Plotting on Windows Using OpenGL

Jan 10th, 2015 | Comments

Introduction - Why on earth?

There are plenty of plotting libraries for Windows using WPF / Win Forms like Oxyplot which can be configured to work with real time data. But rather than using yet another framework, I thought it would be a great idea to do some plotting of my own while also learning something about OpenGL and modern GPUs which I knew nothing about. This post is a result of that endeavour so please bear in mind that if you just want to draw a graph in a WPF application, you don’t really have to worry about OpenGL; it can be done quite easily with a Canvas element, or GDI if you are using Win Forms.

Now that all that’s clear, lets get into the lovely details.

Example Application

The example application plots real time data received from an accelerometer that’s mounted on an embedded device. The data is streamed to the PC via a Serial COM port which happens to actually be over Bluetooth. But any COM interface (USB-CDC / RS232..) will suffice. The accelerometer chip used is from ST Microelectronics - LIS3DH

A video of the application running on my PC is here:

The project includes a Simulation Mode which allows you to test the plotting without actually using a sensor, the plot just loops through some pre-recorded data. Just select the “Simulate” CheckBox and the “Start” button to start the simulation.

Download the example code here. Which you should be able to run immediately if you have Visual Studio 2012 or above.

(Modern) OpenGL in 10 minutes

I assume that like me you don’t know anything about OpenGL. If you do, then you wouldn’t be reading this anyway right? There are plenty of good tutorials on OpenGL written by those who know a lot more than me so I shall link to some that I found useful. You only need to know the absolute basics to be able to do a 2D plot, so once you reach the point where you can draw a Triangle, you know (one point) more than required and you can return back to this post. Most of the tutorials are in C++, but don’t worry about that just concentrate on understanding the principles.

Assuming that you have glanced through one or more of the tutorials above, I’ll briefly summarize the concepts relevant to the plotting scenario.

OpenGL is basically a specification that lays out rules on how software can interact with a GPU chip (you can actually think of it as a GPU driver specification). It tries to strike a balance between being low level enough to be able to utilize all of the GPUs bells and whistles while being high level enough to be portable across different GPU chips. When GPU chip maker X releases a new chip, he/she makes sure that the drivers that they release for that chip contain APIs that are compatible with the latest OpenGL. specification. On Windows, drivers also comply with the Direct3D specification.
The OpenGL interface specifies a lot of flags which encapsulate the State or Mode of the graphics engine. Typically you change whatever flags you need, then pass the vertices and other properties of the object you want to draw to the GPU and then issue the commands to draw the objects, that will be modified according to the current state (color, alpha etc).
The GPU actually has a lot of little processors called Shaders that run at different stages of the drawing pipeline, and can run in parallel. OpenGL allows you to program these shaders via a C like language (GLSL). Generally you will create a separate shader program for each stage in the drawing pipeline. At minimum you need two shader programs - one for the Vertex shader (affects vertex properties) and another for the Fragment shader (affects pixel properties).
Shader programs are compiled and linked somewhere during the run-time of your application, usually at initialization. You can switch the programs as and when you wish too. There are two kinds of variables that you can pass from your main application to a shader program
- Attributes: These represent variables that is uniquely defined for every data point that’s input to that shader, for example every pixel in a fragment shader has a unique color/alpha value.
- Uniform variables - These variables are constant for a particular draw call. Lets say you want to scale all the vertices by a factor of 2, you would store 2 in a uniform variable.
The output of a Shader program is usually the setting of some pre-defined OpenGL variable (prefixed by gl ) or output some pre-defined value. For example the vertex shader program must minimally set the variable gl_Position for all vertices which denotes the position in space of each and every vertex that you draw. And the fragment shader must return the color for each pixel.
The fundamental data structure in OpenGL is an Array (of float[]). These arrays store everything from Vertices to colors of pixels.
Inside the application you will typically store vertices, colors etc inside Buffers/Buffer Objects. Usually you have one buffer in your main application for each attribute in the shader program. You can use an Array Object to encapsulate (bind) all the Buffer Objects related to a particular scene. So when you want to draw something you just have to bind the respective array object to the OpenGL driver and then ask it to draw the scene.
When you want to draw stuff, you simply have to transfer all of these attribute buffers to the GPU and ask it to draw the image by telling the GPU the relationship between the buffers you sent and the attributes in it’s shader programs.

Wrappers for OpenGL on Windows

There are many wrappers that allow you to access the OpenGL API on Windows via C++ or C# like GLUT, GLFW, OpenTK, SharpGL etc. Except for GLUT, all of them can be used with C#/.NET. In this example I’m using SharpGL as it has built-in support for integrating OpenGL context within a WPF application. To use SharpGL just install the package with NuGet (Install-Package SharpGL.WPF) along with GlmNet which is useful for Math stuff. SharpGL also has Visual Studio project templates that you can use to create example projects.

Even though we use SharpGL for this example, as it’s basically just a thin wrapper over OpenGL it is easy to port the code to other frameworks like OpenTK.

Application Design

These are mainly 3 components to the application - the WPF User Interface that includes all the OpenGL plotting magic, a serial communication thread that collects data from the embedded board and a small module to interpret the data received from the accelerometer into X/Y/Z “points”. We shall look at each of these in turn

User Interface and Plotting

The UI consists of two areas - one occupied by the SharpGL control where we will do the plotting and another where we will have some simple buttons to ask the device to Start/Stop communication and a ComboBox to select the COM Port we want to communicate over. Data binding can be used to fill the combobox with the available serial ports by binding the contents with the GetPortNames method in System.IO.Ports :

<ObjectDataProvider ObjectType="{x:Type ports:SerialPort}" MethodName="GetPortNames" x:Key="portNames"/>
<ComboBox  ItemsSource="{Binding Source={StaticResource portNames}}" SelectedIndex="0"  x:Name="cmbxPorts" />

Since each data point from the accelerometer has 3 axes, the display surface available to the SharpGL control should be split into 3 regions. Also as far as OpenGL is concerned the entire 2D display region ranges from the X and Y values of [-1,1] so we need to convert every acceleration point into this range (normalization) before we push the data onto the VBO.

Vertex and Fragment Shaders

The Vertex Shader takes in two 3D vectors called vPosition and vColor for every buffer object and updates the variables gl_Position and color based on the information in vPosition and vColor. color is then passed on to the Fragment Shader which is used to generate the color for every pixel.

// The Vertex Shader 
// (1) takes in a 3D point and adds 'w' to it 
// (2) takes in a 3D color variable and adds Alpha to it and sets a custom variable called 'color'
#version 330

in vec3 vPosition;
in  vec3 vColor;
out vec4 color;

void main()
{

    gl_Position =  vec4(vPosition, 1.0);
    color = vec4( vColor, 1.0);
}


// The Fragment Shader outputs takes in color and outputs color

#version 330

in vec4 color;
out vec4 outputColor;

void main()
{
    outputColor = color;
}

Plotting Data

Plotting the background grid

Since the the grid is just 4 parallel lines that divide the screen into equal portions, it’s easy to draw first. This is done in the method void CreateAndDrawGrid(OpenGL GL) First of all we need to figure out the vertices to be drawn which just consists of 8 points to give 4 lines to split the figure into 3 equal regions. The points are - [(-1,-1), (1,-1)],[(-1,-1/3), (1,-1/3) …], done in the snippet

for (int i = 0; i < (C_NUM_Y_DIV + 1) * 2; i = i + 2)
{
    xgrid[i].x = -1f; xgrid[i + 1].x = 1f;
    xgrid[i].y = C_X_MARGIN_MIN + C_Y_STEP * i / 2; xgrid[i + 1].y = C_X_MARGIN_MIN + C_Y_STEP * i / 2;
    xgrid[i].z = 0f; xgrid[i + 1].z = 0f;
}

Next we need to set the colors or each of the vertices, which is set to White. If you select two different colors then the drawn line will be a color gradient between the colors.

Once the vertices and color buffers are set up, we just need to set up an Vertex Array Object (VAO) containing the two VBOs for vertices and colors, pass this VAO to OpenGL and ask it to draw a figure between these points. That’s what the rest of the code in the method does.

////  Create the vertex array object.
vertexBufferArray = new VertexBufferArray();
vertexBufferArray.Create(GL);
vertexBufferArray.Bind(GL);

//  Create a vertex buffer for the vertex data.
var vertexDataBuffer = new VertexBuffer();
vertexDataBuffer.Create(GL);
vertexDataBuffer.Bind(GL);
vertexDataBuffer.SetData(GL, attribute_vpos, xgrid, false, 3);

//  Now do the same for the colour data.
var colourDataBuffer = new VertexBuffer();
colourDataBuffer.Create(GL);
colourDataBuffer.Bind(GL);
colourDataBuffer.SetData(GL, attribute_vcol, coldata, false, 3);

//  Bind the shader, set the matrices.
shaderProgram.Bind(GL);
shaderProgram.SetUniformMatrix4(GL, "modelview", mviewdata.to_array());

GL.DrawArrays(OpenGL.GL_LINES, 0, 8);

//  Unbind our vertex array and shader.
vertexBufferArray.Unbind(GL);
shaderProgram.Unbind(GL);

Plotting the raw data

Plotting the raw data is quite similar to what was done for the grid, only difference is that the vertices are created by mapping the accelerometer data into the [-1,1] space. The rest of the stuff - VAO / VBO is the same.

The method openGLControl_OpenGLDraw(object sender, OpenGLEventArgs args) is called automatically by the control, so every time control enters here we check if there is a new data point available to plot. If yes, that’s added to the VBO and the whole scene is redrawn. The two BOs used are: vec3[][] _data for vertices of the acceleration points and vec3[][] _dataColor that stores the colors for each vertex.

Rendering is done in the method void CreateAndPlotData(OpenGL GL, float[] raw_data, int figureID) which takes in the new data available and the subplot ID. It converts the points into the [-1,1] space and draws lines through them. Conversion is done by first normalizing the data into the range [-1/3,1/3] and then shifting the data by a fixed factor to align it into one of the subplots. i.e.

for (int i = 0; i < point.Length; i++)
{

    // 1. Normalize the point.x to range [-.33,.33] so we can fit in 3 figures in the range of [1,-1]
    _data[figNum][offset + i].y = 2 * (point[i] - Accelerometer.C_MIN_LSB) / (Accelerometer.C_MAX_LSB - Accelerometer.C_MIN_LSB) - 1; // [1,-1]

    _data[figNum][offset + i].y = _data[figNum][offset + i].y * C_Y_STEP / 2; // [-.33,.33]

    if (figNum == 0)
	_data[figNum][offset + i].y = _data[figNum][offset + i].y + C_Y_STEP + factor * C_Y_STEP;
    else if (figNum == 1)
	_data[figNum][offset + i].y += factor * C_Y_STEP;
    else if (figNum == 2)
	_data[figNum][offset + i].y = _data[figNum][offset + i].y - C_Y_STEP - factor * C_Y_STEP;
}

You can avoid having to reload some data by observing that in the BOs for the raw data, two things are fixed and can be loaded at start-up -

The colors for each of the sub-plots are fixed (to Red/Green/Blue).
The X and Z values of all the samples that will be shown on screen. The X values are limited by the total number of points we want to show on screen - C_SAMPLES_PER_FRAME

Serial Communication in Windows

The raw data comes over the COM Port and is received in the class called SerialPortAdapter. It runs a background thread that constantly tries to read bytes from the serial port using the SerialPort.Read() API. If you are running .NET 4.5 you can also read data asynchronously by accessing the SerialPort.BaseStream class, which is also shown in the example. Both methods worked fine for me.

Once one or more bytes are received, it’s pushed into a Circular Buffer and then the new bytes are passed to the Accelerometer class which can interpret the data received. The circular buffer in SerialPortAdapter discards bytes only if they have been successfully processed by Accelerometer. This is important because the Accelerometer actually needs at least 6 contiguous bytes to get a valid acceleration point.

Streaming can be started and stopped by sending a special character to the device.

Parsing Accelerometer Data

The Accelerometer class provides two important interfaces -

Provides an API to interpret the byte stream received over the serial port. Converts these bytes into acceleration points and stores them in a Queue.
Defines the Maximum and Minimum value of acceleration points possible. This allows the Drawing logic to scale the screen appropriately.

If you plan to use any other sensor, you just need to reimplement these things and everything else should scale automatically. The rest of the code could be split into another class called Sensor for a more modular design, if required.

For this example, the received stream consists of a continuous stream of acceleration points and an acceleration point contains 6 bytes, 2 bytes each per axis. The bytes are arranged in Little Endian format. So (X,Y,Z) of (0x0123,0x145,0x167) would be received as the stream: 0x23, 0x01,0x45,0x01,0x67,0x01

Final words

The code can be extended easily to use other sensors with other formats by modifying the parser and the constants in Accelerometer and some of the rendering code. You can also receive data from multiple serial ports by simply creating more instances of SerialPortAdapter and hooking them to the appropriate stream parser.

This example assumes that no bytes are lost in transmission, if a byte is lost then all the subsequent acceleration values could be misinterpreted. So to make the plotting robust to such failures, you would need to send data in frames and introduce headers, length bytes etc.

To make it convinent to use .NET colors and vec3[] arrays directly, I have added two Extension methods to the SharpGL DLL which are defined in SharpGLEx they define overloads for the OpenGL.Color() and VertexBuffer.SetData() methods.

References

Delegates in C#, Objective C and C - Exchanging Data Between Threads / Objects

Oct 14th, 2013 | Comments

Introduction

A common problem in Object Oriented programming systems is how to pass data between Objects asynchronously. If your system is simple and procedural, communication is easy, the Parent Object can pass data to the Child Object as arguments and the Child returns the result. A simple method call will suffice. Things are not so straightforward if the objects are running asynchronously on separate threads and you want to pass data between them. Lets say you want to update your User Interface with the status of some task running in another object. This post looks at some of the ways to solve this problem in .NET C# and Objective C.

First we will look at asynchronous communication from the Child Object to the Parent Object, since this the most common use case. Data exchange from Parent to Child can usually be accomplished procedurally (via a method call on the Child Object) since generally you have a some UI Event or Parent Object Event which spawns the start of some background process.

Definitions and Example Scenario

Let’s think of the problem in terms of a Parent and Child object. By Parent Object I mean any object that instantiates another object called the Child Object. A good example of a Parent is the UI Object - an instance of UIViewController in Objective C or WindowsForm in .NET C#. A child object is any object that is instantiated in the ViewController or Form.

Let’s assume that the parent is a text box, and the child is a class that is communicating with a autonomous Robot. Let’s assume that the Robot sends us a message every 1 second about its whereabouts (over the internet or some other technology). We can easily simulate this Robot by having a periodic timer run in the Child class that spits out some data every 1s. I’ve chosen a timer for simplicity, it could very well be another thread. Our problem is how to display the information the Robot sends us on the UI (and have the user control the Robot).

Traditional Solution in C

If you were programming in C, this problem could be solved by using Callback methods. Callbacks are the standard pattern in event-driven Embedded Systems to pass data from lower layers into application layers. The UI file would define a method that would be called by the Child thread when it receives a message from the Robot. The disadvantage here is that we have ended up tightly coupling the UI with the Robot communication. If we change the UI and don’t want to display the message anymore, we’ll have to modify the Robot communication implementation to not trigger the callback anymore.

To make the implementation more flexible, you can use a higher level of indirection via function pointers. Instead of having a fixed callback function, define a callback function pointer or better yet a table of function pointers. For example:

// Callback functions receive a byte array pointer which points to the start of received data and length denotes the length of valid data in the array
#define C_MAX_RX_CALLBACKS  3u
typedef CALLBACKPTR (void)*(UINT8*rx_ptr, UINT8 length);
CALLBACKPTR RxCallBackArray[] = {CallBackMethod1, CallBackMethod2, CallBackMethod3 };

// Somewhere in RobotCommunication.c
if(Data_Received == TRUE){
    // Notify the application
    for(i=0;i<C_MAX_RX_CALLBACKS;i++)
        if(RxCallBackArray != NULL){
            (RxCallBackArray[i])(rx_data_ptr,length ); // Call the indexed function
        }
    }

}

The advantage here is that CallBackMethod1 could be deleted later on a UI change, and all you have to do is change the CallBackMethod1 entry to NULL. Also with an Array of callbacks, you can send the data to multiple application modules, say the UI module and a filtering module and so on.

Object Oriented Solution

The solution in Object Oriented languages is similar to that of C, in that it makes use of function (“method” to be more exact) pointers, but is much more powerful and easier to implement. At the heart of it is a thing called “Delegate”.

Delegates in C-Sharp

In C# a delegate is a build-in reference type, just like class, string and others. It specifically indicates a reference to a (surprise, surprise) method (or function in C). So it is nothing but a “Function Pointer Type” with some cool features that we will see soon. It can be used to declare any (static or instance) method - with any return value and any number of parameters.

Coming back to the Robot example, what we need to do is create a property with delegete type in the Robot Class. Parent classes that wish to be notified on the reception of a message, must assign specific methods to this property. In the Robot Class we just call all the methods assigned to this property when we receive some data that other classes may be interested in.

// ------------- Class : Robot Communication --------------------
// 1. Declare the delegate to be implemented in the parent class
    public delegate void dataReceivedHandler(byte[] array);
    public dataReceivedHandler dataReceivedDelegate;
// ........	
// 2. Pass the received data to the delegate(s)
    if (dataReceivedDelegate != null)
    {
        byte[] data_to_send = {1,2,3,4};
        dataReceivedDelegate(data_to_send);
    }
// ------------- END --------------------
// ------------- Class : Parent (Form) --------------------
// 1. Assign the delegate to a specific method
	Robot = new RoboComm();
    Robot.dataReceivedDelegate += ProcessDataFromBot;
// 2. Process the data in the method
    private void ProcessDataFromBot(byte[] data)
    {
		// Do something with the data..
		
		// Since the delegate is called on the calee thread, UI updates can only be done using Invoke
        this.Invoke((MethodInvoker)delegate
        {
            listBox_Logger.Items.Add(DateTime.Now.ToString() + "Got some bytes");
        });
    }
// ------------- END --------------------

There’s a lot more you can do with delegates in C# - you can obviously assign and remove delegates at run-time, the methods assigned to a delegate may be static or an instance methods; if at the parent side you do not wish to go through the trouble of defining a method you can use an Anonymous Delegate - which is a lambda expression containing the code you want to execute in the delegate.

// Anonymous delegate version
Robot.dataReceivedDelegate += delegate(byte[] data){
    this.Invoke((MethodInvoker)delegate
    {
        listBox_Logger.Items.Add(DateTime.Now.ToString() + "And some Anon bytes");
        listBox_Logger.SelectedIndex = listBox_Logger.Items.Count - 1;
    });
};

If a delegate declaration specifies a return type and multiple methods are added to a delegate instance (multicast delegate), an invocation of the delegate instance returns the return value of the last method referenced. The return values of the other methods cannot be retrieved (unless explicitly stored somewhere in addition to being returned).

Delegates can be invoked from a class that did not declare the delegate. This may seem unnecessary but is important when you are inheriting a class from another, this allows the child class to trigger the delegate if required.

Events in C-Sharp

An Event is also a delegate, but it is restricted that it can only be called (triggered) from the class where it is declared, so inheriting subclasses cannot be trigger an event. You can work around this by declaring a protected method in the base class containing the trigger, that is invoked by the child class.

Also the .NET framework imposes some more restrictions on the delegate types that can be used as Events - they should “take two parameters, an “object source” parameter indicating the source of the event, and an “e” parameter that encapsulates any additional information about the event. The type of the “e” parameter should derive from the EventArgs class. For events that do not pass any additional information, the .NET Framework has already defined an appropriate delegate type: EventHandler.”

In other words it’s simpler to use events if you don’t have any data to pass, else just use delegates.

Synchronous and Asynchronous delegates

All the delegates we have seen so far are synchronous delegates. i.e. they are invoked synchronously by the callee object. The Robot object must wait for the execution of the delegates to finish before it can continue processing the next message. There may be situations where you want to do an asynchronous call on the delegates, C# has a couple of different ways for implementing this, some more equal than others. From .NET framework 4.0 and above you should use the Task Parallel Library (TPL) for such tasks. It provides a clean and feature rich interface quite similar in functionality to the Grand Central Dispatch (GCD) in Cocoa that we will see later. Using TPL APIs you can queue lambda functions or any delegates for execution at a later point of time. You can also easily receive a return value from this delegate.

The syntax to schedule a callback to be executed asynchronously on one of the system thread pools is:

byte[] test = { 1, 2, 3, 4 };
Task<bool> task = Task.Run(() =>
{
  // This is a delegate function
    return (asyncTplDataReceivedDelegate(test));
});

// This callback will execute once the task completes
// You can use some flags to set different continuations in case the task throws an exception and/or completes successfully
task.ContinueWith(_ =>
{
    // A lambda function that will be called automatically when the task finishes execution
    Console.WriteLine(" The result in TPL callback is: {0}.", task.Result);
});

// At any point, accessing the task.Result property causes the execution to block until the task is complete and the result is available.

If for some reason you cannot use TPL, then you can use BeginInvoke to achieve the same result; except that it’s a very ugly and messy way of going about it. Once upon a time, it was the only option though so I’ve described next it since I had gone thtough the trouble of figuring it out. Please skip it if you can :)

(Optional) Asynhronous delegates using BeginInvoke (Obsolete, use TPL if you can)

the BeginInvoke method (defined in the CLR for each delegate type) is used for this. Note that this is only applicable for delegates that have only one assigned target method.

When calling BeginInvoke on the delegate, you pass the delegate parameters + two additional parameters. One of them is a callback method which - if it exists - will be executed after the asynchronous method executes. The second is state information that can be passed to the callback, it can be any object. Once a BeginInvoke call is executed, the delegate will be executed on a new thread from the thread pool. BeginInvoke will also return an object implementing the IAsyncResult interface which contains information about the current state of the execution of the delegate method. Then the code execution continues.

How do you know that the delegate has completed execution? You have to use EndInvoke on the delegate. For every BeginInvoke there must be a corresponding EndInvoke. You have to pass the IAsyncResult object returned by BeginInvoke to it and any other ref or out parameters you expect to get back from the asynchronous delegate. Note that the call to EndInvoke is blocking, and only returns when the async operation has ended or an exception is generated.

This gives you 3 ways to get the result from the asynchronous call -

Wait-till-done Just call the EndInvoke once you have completed doing your tasks, and wait till the async operation completes.

  byte[] test = { 1, 2, 3, 4 };
  // Wait - till - done pattern
  IAsyncResult iar = asyncDataProcessDelegate.BeginInvoke(test, null, null); // Start async
  Console.WriteLine("Doing stuff");
  bool result = asyncDataProcessDelegate.EndInvoke(iar); // Wait for end and get result
  Console.WriteLine("After EndInvoke: {0}", result);

Polling method Poll the property AsyncResult.IsCompleted to know if the async operation has completed or not.

  IAsyncResult iar = asyncDataProcessDelegate.BeginInvoke(test, null, null); // Start async
  while (!iar.IsCompleted)
  {
      Console.WriteLine("Doing stuff");
  }
  bool result = asyncDataProcessDelegate.EndInvoke(iar); // get result and cleanup

Callback method This presents the cleanest Async pattern in that there is no waiting involved - you pass a callback to be invoked when the Async operation completes, and call EndInvoke from within this callback. You can also pass an arbitrary object to the callback method via the state parameter. The callback must be of type AsyncCallback.

// Option 2 - use a callback, no waiting!!
IAsyncResult iar2 = asyncDataProcessDelegate.BeginInvoke(test, new AsyncCallback(DelegateExOverCb), null);

(Optional) Redesign using Task based Asynchronous Pattern (TAP)

The use of TPL and async/await keyword suggests a redesign of the system based on the Task Based Asynchronous Pattern. What it means that instead of thinking in terms of callbacks, you can actually write your code as if you were dealing with a synchronous system while still not stalling any of the tasks. So you can write the code as if you were polling the Robot for any new data while not starving the parent object. Note that this is only supported from .NET 4.5, and therefore cannot run on a Windows XP system.

For example you could write something like this:

// On the main form
bool _isCancelled = false;
async private void button1_Click(object sender, EventArgs e)
{
    while (!_isCancelled)
    {
	// Get new data
	byte[] data = await Robot.ReadBytesFromRobotAsync();
	// Update the UI
	Console.WriteLine("Got data from Robot: " + ByteArrayToString(data));
	listBox_Logger_TAP.Items.Add(DateTime.Now.ToString() + " " + ByteArrayToString(data));
	listBox_Logger_TAP.SelectedIndex = listBox_Logger_TAP.Items.Count - 1;
    }
}

// Inside the RobotComm class
async public Task<byte[]> ReadBytesFromRobotAsync()
    {
  TaskCompletionSource<byte[]> dataRxTcs = new TaskCompletionSource<byte[]>();
  Console.WriteLine("Awaiting begins ...");
  await Task.Delay(C_COMM_INTERVAL).ContinueWith(_ => { dataRxTcs.SetResult(new byte[] { 3, 4, 5, 6, 7 }); });
  Console.WriteLine("Awaiting ends ...");
  return dataRxTcs.Task.Result;
    }

Note how the design is simplified as we don’t have to keep track of any delegate pointers and the flow of the program is very obvious. In practice of course, everything after the await keyword will be executed asynchronously just as if it were a delegate.

Example Project for C

A Visual Studio project illustrating the above concepts simulating the Robot is here It includes code illustrating all of the above approaches.

Delegates in Objective - C

Objective-C uses delegates extensively to implement the callback pattern, just as C#. However there are some differences in the implementation and some more bells and whistles. So in Objective-C a delegate is NOT a function pointer; but rather an object pointer that points to the object that implements the callback (delegate) methods. Obj-C defines a keyword called protocol which can be used to declare the callback functions that a delegate to a class must implement. By default these functions are mandatory, meaning that all objects that want to be a delegate of a particular class must implement all the callback functions that the class expects. Optional functions must be declared using the optional keyword.

For our example, we have a View Controller class that is the parent and wishes to be the delegate of the RobotSim class. RobotSim class will call the View Controller’s implementation of the protocol.

  // RobotSim.h
  // In RobotSim declare the “protocol”. Any delegate of RobotSim must implement all the callbacks declared here (unless they are optional). 
  @protocol RobotCommDelegate
  // Required APIs go here:

  - (void) handleRxedDataSync: (NSArray*)array; // We'll encapsulate the received data in an NSArray for simplicity. This is not the most efficient idea.

  @optional
  // Optional APIs go here:
  @end

  // Now we need a property to hold the pointer to our delegate
  @interface RobotSim : NSObject

  // A property to hold the pointer to my delegate
  @property (nonatomic, weak) id <RobotCommDelegate> delegate;

  @end

  // RobotSim.m
  // In response to some event, trigger the delegate
  if ([delegate respondsToSelector:@selector(newFrameReceived)]) {
  [delegate handleRxedDataSync:data];
  }


  // ViewController.h
  // Tell the compiler that VC wants to be a delegate of RobotSim by using the <> brackets
  #import "RobotSim.h"

  @interface ViewController : UIViewController <RobotCommDelegate>

  @property (nonatomic, strong) RobotSim myRobot;

  @end

  // ViewController.m
  // Assign the delegate pointer
  myRobot.delegate = self;

  // Implement the callback which is part of the protocol
  - (void) handleRxedDataSync: (NSArray*)array
  {

  // blah blah blah
  }

That’s all there is to it. An interesting corollary of the delegate being an object pointer is that, the delegate class can cheat and call any function defined in the parent class if required (even if it’s not part of the protocol), and can also use this reference to access public properties of the parent class.

So this is an example of a “Synchronous delegate”. The callback will execute in the delegate object’s context, and it’s a blocking call as far as the delegate thread is concerned.

Asynchronous delegates using the Grand Central Dispatch (GCD)

Apple has made it really easy to live asynchronously with the GCD. To execute a function asynchronously, all you need to do is to pass it as a block to a GCD queue, and the framework takes care of almost everything else. If you want to be notified when the function has finished executing, there’s a GCD API for that too. Though there’s a lot to be said about GCD, we will just concentrate on what necessary to implement Asynchronous callbacks.

Steps to execute an asynchronous method are:

Make a Queue At some point in your application create a separate queue for the method to run, so that you don’t block the main thread. Or you can use one of the “Global Queues” via dispatch_get_global_queue(). Creating a custom queue is good if you want to run code you push onto it serially, FIFO order (not LIFO). If you use one of the global queues, then code runs concurrently and there is no guarantee that the block you pushed onto it first will execute first. The main (UI) queue can be accessed by calling dispatch_get_main_queue(), but you are generally better off not using it.
Pass the method to the Queue To execute the callback asynchronously, just use dispatch_async(), and pass it the Queue and Callback you want to execute.

      dispatch_queue_t myQ;
      myQ = dispatch_queue_create("com.example.myQ", NULL);

      // callback
      dispatch_async(myQ, ^{

      if ([delegate respondsToSelector:@selector(newFrameReceived)]) {
	  [delegate handleRxedDataSync:data];
      }

      });

Call the completion handler If you want to be notified asynchronously of when the delegate has finished processing the data, then you simply have to call the completion block after you executed the call on the delegate. it’s better to use the global queue

      dispatch_queue_t myQ;
      myQ = dispatch_queue_create("com.example.myQ", NULL);

      // callback
      dispatch_async(myQ, ^{

      dispatch_retain(myQ);

      if ([delegate respondsToSelector:@selector(newFrameReceived)]) {
	  [delegate handleRxedDataSync:data];
      }

      dispatch_async(dispatch_get_global_queue(DISPATCH_QUEUE_PRIORITY_DEFAULT, 0), ^{ [self completionHandler];}); // Delegate has completed its work!

      dispatch_release(myQ);

      });

And that’s all there’s to it!!

A minor note: if your delegate happens to be running on the main thread, AND you only have one parameter to pass via the delegate AND it does not return anything; then a simpler alternative without using GCD is to call performSelectorOnMainThread:

Delegate pointers and ARC / Garbage Collection

In the Objective-C example above, you might have noticed that the delegate property was declared as a weak reference. This is important as the Automatic Reference Counting scheme of memory management employed by Objective-C cannot figure out circular references. What has happened here is that the ViewController object holds a reference to RobotSim via the property myRobot, and myRobot in-turn holds a reference to ViewController via it’s delegate property, thus setting ViewControllers retain count to 2. So when you or the system tries to deallocate the ViewController object, it won’t work because it’s reference count is non-zero (1). If you declare the delegate as a weak reference, then the compiler will take care that the ViewControllers retain count will not be incremented when you do myRobot.delegate = self;.

For C#, since the CLR does not use ARC but something called mark and sweep garbage collection which is able to release circular references, so you don’t need to do anything special.

Example Objective C Project

Since the iPhone comes built-in with a lot of nice sensors, we will use this data to replace the hypothetical Robot. To Illustrate delegates, I’ve chosen to read the accelerometer + magnetometer + gyro data available from the Core Motion framework. And to demonstrate Asynchronous delegates, I’ve chosen to read data from the GPS via the Core Location framework. There’s a class called RobotSim which reads the accelerometer, gyro and magnetometer data every 0.5s via a timer. This data is passed on to the ViewController using a synchronous delegate. The GPS location data is posted asynchronously using the Grand Central Dispatch methods described above. The project can be downloaded here.

References:

Adding Bluetooth to Your Embedded Device

May 31st, 2013 | Comments

Introduction

So you are tired of using USB cables or even USB – UART converters and want to have some wireless action on your device. You think bluetooth might be the best option but are confused by all the layers on the stack and the various versions like BT2.0, BT4.0, BLE and so on. This post will try to shed some light on this and summarize the options available and how difficult or easy it can be to add BT to your device.

Intro to BT versions

At the moment the main BT versions common are BT 2.0, 2.1, 3.0 (High Speed) and 4.0 (latest version). Most Host devices at the moment support BT2.1 while newer devices like iPhone 4S have BT4.0. The range of a BT device depends on the “Class”, Class 1 being the maximum range and Class 3 the lowest. For indoor applications, with a range of 10m or less, Class 2 is the answer.

Now you can think of BT4.0 as being the sum of BT 2.0/2.1/3.0 (called “Classic Bluetooth”) and Bluetooth Low Energy (BLE). And the cool thing of BLE is that it is an entireley new and simpler protocol from the ground up, including the physical layer. It specifically targets applications where you don’t need to stream a lot of data, but you just want to send some small data periorically or so. The problem with BLE right now is that most Phones/PCs do not support it yet.

Another thing worthy of mention in the Classic BT vs BLE debate is that if you want to have a device that can talk to iPhone/iPad/iPod, then it is better to use BLE because you do not have to sign up for the Made For iPhone program to use BLE

BT Soc vs Module …BT2.0

Whatever BT version you want to add to your device, you have basically two options – use a BT SOC chip from Ti, CSR or someone else OR use a BT “module” from Roving Networks, Panasonic or someone else. The difference is that the SOC is just the chip with the RF front end and implementing the BT stack up to the HCI interface, while the module will come with an antenna and the entire BT stack up to some “Profile”, usually the Serial Port Profile (SPP). SPP is what we need to transmit data bytes between BT devices.

SOC

The SOC is cheap (1-3$) but it’s only the chip – so it is your responsibility to design the antenna, get the BT certifications and in the case of BT 2.0 – 3.0 write a significant amount of software (the BT stack from the HCI layer upwards to RFCOMM/SPP) before you can transmit a byte of data from your favourite mobile phone to your device. In other words, using BT 2.0 – 3.0 on your device by adding an SOC is not for the weak of heart.

If you do decide to go the SOC way, I have found an open source stack that actually works - I have tested it - on PIC32. You should check it out. It is entirely event based and build upon the existing Microchip USB drivers. The whole stack compiles into less than about 8KB ROM, though it does require a heap around 4kB.

Modules

The simpler option is to buy a ~15$ module from Roving Networks or Panasoic or someone else. These modules have the entire stack up to SPP built-in and also remove the burden of certification from you. All you have to do is plug in the module to your PCB and configure it with your MCU by sending some AT commands over the UART. You can almost immediately start working on your Application as if it were a serial port.

BT Soc vs Module …BLE

The SOC vs Module question plagues BLE too. But in this case the burden of implementing the stack is no longer upon you as the BLE stack is much simpler than Classic BT. Even the SOC chips come with a built-in stack but some of try to trip you by forcing you to buy IAR compilers and other usless stuff. All this is documented extensively here.

If you decide to use a module for BLE, then you would want to look at Murata or Bluegiga.

In the next posts I will post some of the code I have been using to talk to my Roving Networks BT2.0 enabled device from Windows, Linux and Android.

Integrating Chan’s FatFs Library With Microchip Applications Library’s USB Stack

May 31st, 2013 | Comments

Introduction

The Microchip Applications Library has a lot of great code that is easy to integrate and can bump up the feature set of your product in hours. So far I have had the opportunity to use some of the USB (MSD,HID,MS-Host) stacks along with their FAT library (FSIO). The USB stacks, along with their well written application notes are very good and quite easy to integrate, even if you know nothing about USB. Hovever, their FAT library is not the best out there. It takes up quite a lot of code + data space, is slower than Evolution (when there’s no external duress) and not extensible – for instance to support multiple volumes or long file names.

FatFs can solve all these problems. It seems to the handiwork of a one samurai army, who has done a lot of other cool projects too. In spite of doing all this stuff that I would gladly sell off my left testicle for, he actually refers to himself as a “mediocre” embedded systems engineer!! Anyhow, it’s also released with a free, open source, do-whatever-crap-you-can-with-it license.

Groundwork

This post will explain how to integrate FatFs to a Pic24 project. Sample code is also available for download. Most of the information contained here is just an aggregation of the things I found trawling the posts on Microchip forums.

Before going into the details, the files involved and their function should be evident. This is summarised in the table below. Note that both FatFs and FSIO have separated the FAT layer from the lower disk access layer. So you can build the system on different types of memory, furthermore in the case of FatFs your application can manipulate multiple memory types simultaneously. In this example I will only consider SD cards (over SPI) and USB Mass storage devices. You should also download the latest FatFs library and the sample

Function	MAL Files	FatFs Files
FAT File System	fsio.c/h + FS config files	ff.c/h + ffconf.h
SD card access	sd-spi.c/h	diskio.c/h

1. Integrating FatFs (with an SD card access)

In this scenario, your PIC is interfaced to an SD/MMC card via SPI and you have to set up a FAT system to access data on the SD card. You should not have to change anything to the ff.c/h files. But you have to add a call to the timer function and implement a function to return the timestamp for the files. Then you have to modify the diskio.c routines to match your particular SPI port. There are a few sections in Chan’s code marked as “Platform Dependent” that you need to modify.

2. Integrating FatFs with the MAL USB Library

Your device is a USB Mass Storage Device Something that might not be apparent on first glance is that for the case of a USB MS-device, you DO NOT need to implement a FAT system on the Device. The USB Host will take care of implementing the file system and/or formatting the device. So all you have to do is plug the MAL USB MS-Device routines to your diskio.c/h routines (see usbdevFatFs.c/h in the example code).
Your device is a USB Mass Storage Host: In this scenario, your PIC is a USB Mass Storage Host and wants to read/write data on a USB Mass Storage Device like a thumb drive. Now in this case, your application does have to implement a FAT (or some other) file system on TOP of the USB layers in order to read / write the thumb drive. The difference between the case 2.b and 1. is that the diskio.c routines are different and they access the USB BUS and the SPI BUS respectively.
Finally, What if your PIC talks to an SD card, and also acts as both a USB Mass Storage Device and Mass Storage Host? This is a case when you can use the multiple volume support in the FatFs (one volume is the SD card and the other is the USB MSD thumb drive). All you have to do is pass the volume number to the FAT routines, and you have to modify diskio.c to call the appropriate low level routines based on the passed in drive volume number (see the diskio.c file in the example code). The example code is based on the USB Mass Storage host+device usecase, it’s tested and it works!

ExampleCode

← Older Blog Archives

A	B	C	D	E	F	H	I	J	K
91	1	1	1	1	1	1	1	1	1
90	1	1	1	1	1	1	1	1	0
10	10	10	10	10	10	10	10	10	10
20	20	20	20	20	20	20	20	20	20
1	20	11	28	0	0	20	20	0	0
0	0	11	28	0	0	20	20	20	1
1	1	1	1	1	1	1	1	1	1
1	1	1	1	1	1	1	1	1	1

A	B	C	D	E	F	H	I	J	K
91	1	1	1	1	1	1	1	1	1
90	1	1	1	1	1	1	1	1	0
10	10	10	10	10	10	10	10	10	10
20	20	20	20	20	20	20	20	20	20
1	20	11	28	0	0	20	20	0	0
0	0	11	28	0	0	20	20	20	1
1	1	1	1	1	1	1	1	1	1
1	1	1	1	1	1	1	1	1	1

A	B	C	D	E	F	H	I	J	K
91	1	1	1	1	1	1	1	1	1
90	1	1	1	1	1	1	1	1	0
10	10	10	10	10	10	10	10	10	10
20	20	20	20	20	20	20	20	20	20
1	20	11	28	0	0	20	20	0	0
0	0	11	28	0	0	20	20	20	1
1	1	1	1	1	1	1	1	1	1
1	1	1	1	1	1	1	1	1	1