Linear 3D triangulation

I came across this 3D linear triangular method in TheiaSFM:

bool TriangulateNView(const std::vector<Matrix3x4d>& poses,
                      const std::vector<Vector2d>& points,
Vector4d* triangulated_point)
  CHECK_EQ(poses.size(), points.size());

  Matrix4d design_matrix = Matrix4d::Zero();
  for (int i = 0; i < points.size(); i++) {
    const Vector3d norm_point = points[i].homogeneous().normalized();
    const Eigen::Matrix cost_term =
        poses[i].matrix() -
        norm_point * norm_point.transpose() * poses[i].matrix();
    design_matrix = design_matrix + cost_term.transpose() * cost_term;
  Eigen::SelfAdjointEigenSolver eigen_solver(design_matrix);
  *triangulated_point = eigen_solver.eigenvectors().col(0);
  return == Eigen::Success;

I was aware of the DLT (direct linear transform), but it didn't look like any formulation I've seen before. It's actually pretty neat. Let's say you're trying to find an unknown homogeneous point in 3D, \mathbf{X} = [X, Y, Z, 1]. What we have is N poses, P, represented as 3\times 4 matrices and the corresponding 2D coordinates represented as homogeneous points in \mathbb R^3. The 2D points are written as \mathbf{x} = [ x, y, 1].

Since we're triangulating the 3D point, and we have homogeneous coordinate (i.e. \alpha \mathbf{x} \equiv \mathbf{x}) then for all i we should have:
\alpha_i \mathbf{x}_i \approx P_i \mathbf X
given an scale factor \alpha.

Now let's pick apart the code above. Let's call design_matrix D and cost_term C. On line 12, we have:
\displaystyle D = \sum_{i=1}^{N} C_i^\top C_i
And line 15 we’re finding the eigenvector corresponding to the smallest eigenvalue of D (SelfAdjointSolver produces them in a sorted order), i.e.
\mathbf{X} \approx \displaystyle \underset{\mathbf{v}, |\mathbf{v}|=1}{\text{argmin}}\ \mathbf{v}^\top D \mathbf{v}

We can rewrite D = \mathbf{C}^\top\mathbf{C} where:
\mathbf{C} = \left[ \begin{matrix} C_1\\ C_2\\ \vdots\\ C_N\\ \end{matrix}\right], which substituting in above gives:
\mathbf{X} \approx \displaystyle \underset{\mathbf{v}, |\mathbf{v}|=1}{\text{argmin}}\ \|\mathbf{C v}\|_2^2,
which is of course the right singular vector corresponding to the smallest singular value of C. Using eigen decomposition is much more efficient the size is O(1), not O(N), but probably at the penalty of less numerical precision.

Either way we’re trying to find the approximate nullspace of \mathbf{C}, which means finding something that’s roughly in the null space of all the C_is. But why?

On lines 8–11, we have:
C_i = P_i - \mathbf{\hat{x}\hat{x}^\top}P_i,
and we’re claiming \mathbf{X} is about in the null space. Let’s see what happens when we multiply by it:
(P_i - \mathbf{\hat{x}\hat{x}^\top}P_i) \mathbf{X} = P_i \mathbf{X} -\mathbf{\hat{x}\hat{x}^\top}P_i \mathbf{X}\\
Now, substituring in the first equation we have all the way back at the top gives:
\approx \alpha \mathbf{x} - \alpha\mathbf{\hat{x}\hat{x}^\top x} = \alpha \mathbf{x} - \alpha\mathbf{\hat{x} |x|} = \alpha \mathbf{x} - \alpha\mathbf{x} = 0
taa daa!

So there you go! If there is no noise, \mathbf{X} is in the right null space of C_i and so the null space of \mathbf C and of course D. If there is noise, then it’s closest to being in the null space of all of C_i measured in a least squares sense.

Note though that it’s just an algebraic error, not a reprojection error so it will give slightly oddly scaled results which will give more weight to some points than others. It is however a rather elegant solution.

Neutrogena light mask part 3: about that LCD…

In Part 1 I hacked it to get 99 lives.
In Part 2 I did an excessively thorough analysis of the current limiting.
In Part 4 I found it keeps resetting on life 83.

The light mask comes with an LCD. I’ve always been curious about driving them, but never really taken the time to look. So, I figured I’d get my scope out and have a go. First off, it’s two seven segment displays. Rather handily, Neutrogena have left nice large test points for me. And if you look closely they even saved me the hassle of counting the number of them! There are 8. So, it must be a multiplexed display. Assuming it’s even vaguely symmetric, it’s got to be something like 5×3 or 4×4 (with some dead spots). So, time to break out the scope! First though, I have to solder some wires on:

Wiring the LCD into a breadboard by using the test points.

Except they don’t fit under the LCD. Plus if you connect the power up backwards, then it appears to stop working. Why knew? Fortunately I have one more spare. The failure mode seems to be that one digit is no longer working (the LCD works though—I rubbed some power supply leads across the zebra strip and they all lit up.). Weird.


Yes OK, so it’s actually working fine (which is weirder) it’s just that it displays “0”, not “00” because it’s made for end users who aren’t expecting everything to be in nicely formatted hex…

So, I’m vaguely aware of some things about LCDs, in no particular order:

  • They activate off DC, but that eventually destroys them so you should use AC.
  • Driving simple (non multiplexed) LCDs is easy.
  • Driving multiplexed ones is harder than you’d expect .
  • There’s a threshold voltage below which they don’t turn on.

And here’s what the waveform for the first 4 pins looks like:


I blame my scope. No way it’s a PEBKAC error. There are like, 4 levels and they go all over the place. It’s crazy.

OK, that is indeed harder than I’d expect. Reading around was a bit boring, confusing and badly written. There’s some information which indicates it’s possible to drive them off only 2 levels. This would make sense: if you want to make a particular row/column pair light up, then you drive those two out of phase. Everything else…

Oh I see.

So, I guess the other rows have to be driven in-phase with the column, and the other columns… hmm. OK, I guess that’s why they have multiple levels. If you have a 5V supply, and you drive the row/column out of phase, the one intersecting segment sees 10v pk-pk. If you drive the other rows and columns with an idle voltage (say 2.5V DC) then the segments see either 0V if neither the row/column is driven or 5V pk-pk if a row/column is driven and the column/row is idle.

Backing up, imagine driving a matrix of light bulbs (not LEDs because I don’t want to get bogged down in Charlieplexing) from switches:


A multiplexed display matrix with one rwo switch and two column switches closed.

Switches are either closed, which on a microcontroller means driving the pin as 1 or 0 depending on which rail the switch is connected to, or open which means tristated. For current driven things, like light emitters, it’s easy: tristated means no current flows. For voltage devices, it simply means “not actively driven”, so something needs to be done with it, i.e. bias it to some neutral voltage.

I have no idea which things are rows and which things are columns on the display. However, I do know that the voltages are 0, 1.6, 3.3 and 5V. I’d hazard a guess that 3v pk-pk (i.e. a neighbouring pair) won’t drive the display but 10V pk-pk will. Not sure about 6.7. Probably?

Well, I’ve got an Arduino and some voltage regulators. For a 5V drive, I can easily get 0, 2.5V and 5V by tristating an output and pulling it to a central value:


Pulling the pins to a 2.5V rail made with a potential divider and a TS358 opamp (costs 6p!)

I used 47k because it’s high ish (won’t load the Arduino or opamp much) and I have plenty of them. Anything in the range of about 1k to several meg would probably work OK.

I could use the 3.3V output as the central value but honestly that seems to be tempting fate and I don’t know if it can sink current. Instead, I’ll use the entertainingly cheap TS358CD op-amp. So, time to cut out the useless remains of the device and wire up the Arduino!

Tracks cut with a stanley knife

I can still fix this!

Also, the test pads aren’t all they’re cracked up to be: the joints need to be filed down. Even now I’m not sure I’m getting perfect contact between the LCD and the board.


I had to file the tops of the joints down very carefully (totally unnecessary as it transpired).

Anyway I’ve wired up pins 1-8 on the arduino (not 0!) to 1-8 on the LCD, more for my own sanity than anything else. And with a simple program:

void setup(){}

void loop() {

  pinMode(1, OUTPUT);
  digitalWrite(1, HIGH);
  digitalWrite(1, LOW);
  pinMode(1, INPUT);

I get this scope trace:SCR03

This had me going for quite a while. The pin set to input should be getting 2.5V.  But it’s not; it’s being pulled up. I looked up: the internal pullup is 20k. The voltage is consistent with that: 2.5 * 47 / (47 + 20) + 2.5 ≈ 4 ish.


Well, that had me going for a while. I went back and forth with trying to figure out how to turn the pullup off (it’s off by default), turning it on and seeing what happened, plugging and unplugging wires and all that. Turns out that I was using pin1 which on the Arduino is the TX pin for the serial port if you want it to be. That means it has an LED attached which is doing extra pullup effectively to almost exactly the right value.

Ho ho ho.

So looks like I’ll be using pins 2-9 instead and I won’t get to keep my sanity. But at least that works.  Also I realised after a, er, little debugging that the reason the device has screws next to the LCD is so that they push the board against the zebra strip ensuring good contact. I wonder if I should use those screws…

Anyhow, now I think what remains is to so a somewhat exhaustive search over all pairs of wires driving them in opposition, to see what happens when they’re activated.

static const int A=2, B=3;
void setup(){
  pinMode(B, OUTPUT);
  pinMode(A, OUTPUT);

void loop() {
  digitalWrite(A, 1);
  digitalWrite(B, 0);
  digitalWrite(A, 0);
  digitalWrite(B, 1);

That only took a few minutes: it’s only actually 28 combinations. Here are the results I’ve noted down, along with the matrix that’s implied by the numbers. I’ve written the numbers as pairs indicating the two pins that have been activated:

Oh actually! It’s even better! This proves that not only is 10v pk-pk sufficient to drive the segments (I was sure of this), but 5v pk-pk isn’t, which I wasn’t so sure about. That’s nice: no extra circuitry is required. So, what we have is a 4×4 matrix. What I’m going to do is drive each of the rows in turn, while driving all 4 columns simultaneously.

The mapping is very regular though and we actually have essentially two 4×2 matrices for the two digits. The plan: each digit will be represented by a 7 bit number, one bit for each segment. Then, a 4×2 matrix will be represented as 4 2 bit numbers. The next step is a little tedious and involves designing the segment pattern for each digit and figuring out the mapping of segments to columns for each row. I’ve gone for A-Z (excluding k,m,v,w,x), picking the lowercase letter in general when there’s no other criteria and of course 0-9:

And that’s most of it. All that remains is to write a program which drives the rows in sequence and sets the correct column pattern for each row. Also, I’m going to have a function that maps ASCII to the segment patterns so I can write stuff on the display. My choice of driving is to drive each row low in turn, then repeat with each row high in turn to make it AC. I did try row 1 high then low, then row 2 high then low etc too but it didn’t make much difference. Here’s he code:

void setup(){};

//Decode uppercase ASCII and numbers into
//a list of segments, stored as a bit pattern.
uint8_t ascii_to_code(uint8_t c)
  static const uint8_t letters[26] = {
  static const uint8_t numbers[10]={

  if(c >= 65 && c = 48 && c <= 57)
    return numbers[c-48];
    return 0;

//Decoders giving the first 2 column entries
//for the 4 rows, given the pattern for a
//single digit
uint8_t r5(uint8_t n){
 return (n&1)1;
uint8_t r4(uint8_t n){
 return (n&32)>>4 | (n&64)>>6;
uint8_t r3(uint8_t n){
 return (n&16)>>3 | (n&4)>>2;
uint8_t r2(uint8_t n){
 return (n&8)>>2;

//Set the column outputs as driven or inactive
//according to the input bit pattern. Polarity
//determins whether to drive high or low.
void set_columns(uint8_t n, bool polarity)
  for(uint8_t i=0; i < 4; i++)
      pinMode(i+6, OUTPUT);
      digitalWrite(i+6, polarity);
      pinMode(i+6, INPUT);

void display_digit(uint8_t left, uint8_t right)
  //Columns entries for both digits for
  //the 4 rows.
  uint8_t rows[4];
  rows[3] = r5(left)<<2 | r5(right);
  rows[2] = r4(left)<<2 | r4(right);
  rows[1] = r3(left)<<2 | r3(right);
  rows[0] = r2(left)<<2 | r2(right);  

  //Do positive columns/negative rows first
  //then repeat with the polarity inverted
  //Activate each row in turn and set the
  //column entries for that row.
  for(int p=0; p  100){

There was a lot of fiddling around, most of which did very little. Basically, driving a multiplexed display off 3 levels is pretty marginal. I found that often some digits would ghost on, and others would be very faint. I could increase the contrast by lengthening the amount of time it took to draw the whole display, by driving a row with many cycles then moving on to the next row, but it only got really good when it was far far too slow. I did find putting a slight pause in between driving rows did help. Removing it darkened everything including ghosted on digits, lengthening it lightened everything. The value I got was about the best one I could find.

Here’s what it looks like:

It’s not perfect, you can see the contrast is a bit poor and so “hELLo” comes out looking a bit like “hCCCo”. Nonetheless I’m moderately pleased, since it does kind of work. I have to be careful with the choice of symbol though because they’re not all that good.

Small hack of the day

Two things:

  1. syntax highlighting is nice
  2. wouldn’t it be nice to copy /paste it

The main problem is many terminals don’t for no especially good reason allow copy/paste of formatted text. There are three tools which help:

  1. xclip -t text/html this eats standard input and allows pasting of it as HTML, so it can include formatting and color and so on. By default, xclip does plain text, so you have to specify that it’s HTML.
  2. aha this takes ANSI formatted text (i.e. the formatting scheme used by terminals and turns it into HTML).
  3.  unbuffer many tools will only write color output to a terminal, not a pipe. Through the magic of PTYs (pseudoterminals) this program fools other programs into thinking they’re running on a terminal when they’re not.

Now all you have to do is pipe them together. So, take a silly C++ program like this:

void silly;

And compile it like this:

unbuffer g++-7 -std=c++1z  | aha | xclip -t text/html

Then a simple paste into a box which accepts HTML (e.g. email, etc) gives: error: variable or field 'silly' declared void
 void silly;


Server Send Events

Server send events are a way of getting data pushed over an http connection to a browser. The JavaScript interface is very simple. Here’s some useful info:

I thought it would be a neat way of debugging some GIS code in C++: all I had to do was write the data to a socket, and have a page that collected it in JavaScript and plugged it into Google maps. It was that simple, sort of, except that it was actually incredibly awkward getting it up and running since web browsers are both finicky and don’t provide much error information.

All you have to do is open a socket, write the right HTTP headers and send the correct data. I eventually ended up sending it in chunked encoding, which means each message is essentially preceded by a length so the browser knows how much data to accept and put together into a message.  The alternative is to use Content-Length and have fixed length messages (like the w3schools example), but I couldn’t manage to get my browser (Firefox) to accept more than one message before closing the connection due to an error. No idea why, but the chunked encoding is much more flexible anyway.

Probably the biggest hurdle is that my HTLM page was just a file, but the server send events were from a (local) server. Either way that meant it was cross domain and so Firefox would block the request because of CORS. Turns out the fix was a simple headre tweak but that had me going for a while!

Anyway, here’s the code:

#include <iostream>
#include <cstring>
#include <cerrno>
#include <vector>
#include <iomanip>

#include <unistd.h>
#include <sys/socket.h>
#include <netinet/in.h>
#include <netinet/ip.h>

using namespace std;

class VeryBasicServer
		int sockfd=-1; 
		int connection_fd{-1};

			//Incantation to create a socket 
			//Address famile: AF_INET (that is internet, aka IPv4)
			//Type is a reliable 2 way stream, aka TCP
			//Note this just creates a socket handle, it doesn't do anything yet.
			//The 0 is the socket type which is pretty redundant for SOCK_STREAM (only one option)
			sockfd = socket(AF_INET, SOCK_STREAM, 0);
			if (sockfd < 0)
				throw runtime_error("Error opening socket:"s + strerror(errno));

			//Set options at the socket level (a socket has a stack of levels, and 
			//all of them have options). This one allows reuse of the address, so that
			//if the program crashes, we don't have to wait for the OS to reclaim it, before
			//we can use this socket again. Useful for debugging!
			int true_ = 1; 
			//Binding a socket is what sets the port number for a particular
			sockaddr_in serv_addr = {};
			serv_addr.sin_family = AF_INET;          //internet address family
			serv_addr.sin_addr.s_addr = INADDR_ANY;  //allow connections from any address
			serv_addr.sin_port = htons(6502);        //Still fighting the 80s CPU wars. 6502 > 8080
			if (bind(sockfd, (struct sockaddr *) &serv_addr, sizeof(serv_addr)) < 0) 
				throw runtime_error("Error binding socket:"s + strerror(errno));
			//This is necessary for stream sockets (e.g. TCP). When a client_ connects
			//a new socket will be created just for that connection on a high port.

		void read_ignore_and_dump_incoming_data()
			vector<char> buf(65536, 0);
			int n = read(connection_fd,,buf.size());
			if (n < 0) 
				throw runtime_error("Error reading from socket:"s + strerror(errno));
			cout << "Message reads: ";
			cout.write(, n);

		void accept()
			//Listen blocks until a cnnection is made, then hands over the newly created
			//socket for the connection
			connection_fd = ::accept(sockfd, nullptr, nullptr);

			if (connection_fd < 0) 
				throw runtime_error("Error accepting:"s + strerror(errno));
			//We can actually ignore the HTTP header just to get up
			//and running. For ServerSendEvents, there's only one accepted
			//MIME type, and uncompressed is always a valid compression option
			//even if the client_ doesn't request it.

			//Construct a valid working header.
			//Two important things to note. One is the access control. Since
			//this server isn't serving anything except SSEs, the web page
			//which is using them ust come from elsewhere. Unless we allow 
			//connections from places other than the originating server, then
			//the web browser will block the attempt for security.
			//The other point is the chunked encodeing. The browser connecting
			//has to know when we've finished sending an event in the stream
			//of data. Chunked encoding allows us to send data blocks along with a 
			//length so the server knows when a block is finished. The other option
			//is to have a fixed Content-Length instead. I never got it working,
			//but it's much less flexible so I didn't try hard.
			//Note also the \r\n line terminations, not just \n. Part of the HTTP spec.
			write("HTTP/1.1 200 OK\r\n"                          //Standard HTTP header line
			      "Content-Type: text/event-stream\r\n"          //This is the only allowed MIME type for SSE
                  "Transfer-Encoding: chunked\r\n"               //Chunked encoding lets it know when an event is done without knowing packet boundaries.
				  "Access-Control-Allow-Origin: *\r\n"           //Because the HTML comes from a file, not this server, we have to allow access
			      "\r\n");                                       //End of header indicator

		//Write data with less than the minimal amount of error checking.
		void write(const string& str)
			int n = ::write(connection_fd,, str.size());

			if(n < 0)
				throw runtime_error("Error writing:"s + strerror(errno));

		void write_chunk(const string& str)
			//Chunked encoding is the size in hex followed by the data. Note that both
			//the size and data fields are terminated by HTTP line terminations (\r\n)
			ostringstream o;
			o << hex << str.size() << "\r\n" << str << "\r\n";

			cerr << "closing\n";

int main()
	VeryBasicServer hax;

	cout << "Press enter to send an event" << endl;
	for(int i=1;;i++)
		if(cin.get() == EOF)
		//Build up an event in server-send-event format. The message consists of 
		//one or more fields of the form:
		//field: value\n
		//Followed by an empty line.
		//The main tag is "data:" which carries the data payload.
		//for more info (e.g. different message types and dispatch)
		string sse;
		sse += "data: " + to_string(i) + "\n";  //Send the current number as the data value
		sse += "\n";                            //Empty field ends the message.

		//Send the message data using chunked encoding



And the corresponding HTML page is:

<!DOCTYPE html>
<head>  <meta charset="UTF-8"> </head>

<h1>Getting server updates</h1>
<div id="result"></div>

if(typeof(EventSource) !== "undefined") {
    var source = new EventSource("");
    source.onmessage = function(event) {
        document.getElementById("result").innerHTML += "Data: " + + "<br>";
    source.onerror = function(event) {
        document.getElementById("result").innerHTML += "Connection failed<br>";
} else {
    document.getElementById("result").innerHTML = "Sorry, your browser does not support server-sent events...";


And that’s basically it. It’s a handy and rather silly chunk you can put in a C++ file. It’s a lot of very poor practices and does all the wrong things for scalability, security, generality and good sense, but it’s handy for debugging hacks.

Code is on github.

Blog as you go: sigma delta DAC

I have some piezo speakers from another project. That one used a bi-level amp to drive them. I figured it would  be fun to try a tri-level drive, using an H bridge allows you to have +, – and 0V across the device. And for fun, why not make it a direct sigma delta encoder?

It’s going to run on a microcontroller (an arduino). It’ll need very precise timings, so I’ll not be using the arduino environment.

Here’s a first pass in C++ for linux:

#include <iostream>
#include <cmath>
#include <utility>

using namespace std;

float signal(int t)
	//Roughly 0-1
	return (1 + sin(t/20. - 1))/2.1;

float quantize(float f, int levels)
	return min(levels-1.f, max(0.f, floor(f*levels)))/(levels-1);

int main()

	float integral = 0;

	for(int i=0; i < 200; i++)
		float output = quantize(integral, 3);
		float difference = signal(i) - output;
		integral += difference;

		cout << signal(i) << " " << output << endl;


And this produces output like this (mildly prettified):

Which looks about right. The code will pretty much stink on an arduino since it’s all floating point. It’s easy to convert to integer code though:

#include <iostream>
#include <cmath>
#include <utility>

using namespace std;

uint16_t signal(int t)
	//Roughly 0-1
	return 65535 * (1 + sin(t/20. - 1))/2.1;

uint16_t quantize(int32_t i)
	if(i < 21845)
		return 0;
	else if (i < 43690)
		return 32768;
		return 65535;

int main()

	int32_t integral = 0;

	for(int i=0; i < 200; i++)
		uint16_t output = quantize(integral);
		float difference = signal(i) - output;
		integral += difference;

		cout << signal(i) << " " << integral << " " << output << endl;


I’ve used a uint16_t for the value, which will effectively represent both positive and negative levels with 32768 being the zero point. Note that the error integral must be both signed and wider since the errors can grow beyond the signal range:


Now to port to the arduino. So, I’ll get my Makefile from here.

I’m going to pick pins 8 and 9 on the UNO, that is PB0,1 on the chip for my outputs. and here’s how you get 3 way opposed outputs as +/-/0. To demo, I connected a pair of LEDs in parallel but facing the other way:

#include <avr/io.h>
#include <util/delay.h>
int main()
	int i=0;

			DDRB = 3;
			PORTB = 1;
		else if(i==1)

			DDRB = 3;
			PORTB = 2;
			PORTB = 0;

		if(i > 2)


So I started the port and BOOM! 😦

It stopped working. The reason was simple: the simple makefile takes one object file and converts it to HEX. Since I’m using sin(), we actually need to link the .o and .a into a .elf file, then convert THAT to HEX. The snippet is:

%.hex: %.elf
avr-objcopy -j .text -j .data -O ihex $< $@

prog.elf: delta_sigma.o
avr-gcc $(FLAGS) -o prog.elf delta_sigma.o -lm

Obvious, really, in hindsight…

So, OK, now to convert the modulator code to the arduino. Lots of things went wrong. But first, here’s the code:

#include <math.h>
#include <stdint.h>
#include <avr/io.h>

uint16_t signal(int32_t t)
	float u = t / 1024.f;

	//Roughly 0-1
	return 65535 * (1 + sin(2*3.151592f*u))/2.1;


uint16_t quantize(int32_t i)
	if(i < 21845)
		return 0;
	else if (i < 43690)
		return 32768;
		return 65535;

int main()

	int32_t integral = 0;


	for(uint32_t i=0; ; i++)
		uint16_t output = quantize(integral);
		int32_t difference = (int32_t)signal(i) - output;
		integral += difference;

		if(output == 0)
			PORTB=1; //Output 1 0
		else if(output == 65535)
			DDRB =255;
			PORTB=2; //Output 0 1
			PORTB=0; //Output 0 0

What didn’t go wrong? Nothing! I wasn’t nearly careful enough with my ints (only 16 bits on AVR), ints of specific width, overflow and that sort of thing. Also, initially, I decided to output a 0 level by tri-stating the two outputs, so they both float to the middleish. Turns out that didn’t work well since they float extremely slowly (not surprising really!). Forcing them both down to 0 worked much better.

After all that, I then connected a simple RC filter across it so you an see the results:

That’s actually a pretty nice sine wave there! It ought to be: there’s really not much room for nonlinearity and other distortions to creep in. I’ve zoomed in a few levels so you can see how it looks in detail.

It is however really really slow. I’m using full floating point, and a transcendental operation every iteration of the sigma delta encoder. That is really slowing down the cycle time since the AVR isn’t very fast. That accidentally solves the other problem which I’ve made no attempt to make sure every path takes the same number of cycles. But that sin() is dominating so heavily that it doesn’t matter.

And that’s it for part 1: a working sigma delta encoder. For part 2, I’ll make it fast enough to generate audio tones which aren’t simply the sigma-delta encoder transitions (I hope).

Oh also here’s tehe obligatory github link.


Automatic dependencies for C++ in 3 very simple lines of make

I like Make. Make is best ponbuild tool. If you don’t know it already, it expresses dependencies as a DAG (directed acyclic graph), so a node lists it’s dependencies and has a rule about how to generate the file(s) from the dependencies. If a file is older than its dependencies, then it’s stale and needs to be rebuilt.

If the DAG is specified correctly, then if you change one source file, you’ll rebuild the minimum necessary with the maximum possible parallelism…

…I meant GNU Make by the way, of course.

In Make, dependencies are generally specified statically, like this:

target: dependency1 dependency2

The problem with this is that in many cases dependencies are not static. That is, the dependencies are a function of the dependencies. What? OK, let me illustrate. Suppose you have a dependency:


OK, so the object file foo.o depends on So far, so good. But it’s almost certainly incomplete. If #include’s anything then foo.o also relly depends on that too. In other words with the incomplete dependencies, if you modify that header and type “make”, foo.o won’t rebuild because make doesn’t know it ought to. This has the annoying problem that when you’re modifying headers, you keep having to “make clean” and rebuild everything.

Fortunatley, make provides a solution for you in typical make fashion: it provides the mechanism to deal with dynamic dependencies and you provide the language specific tools. You can just do:

include morestuff

Naturally, morestuff is going to be dynamically generated. GCC is nice and since it knows about the dependencies (it has to when it actually does the compile), and will emit them in make format while it does the build, ready to be included next time. So if a source file changes, the .o is rebuilt and new dependencies are generated. Next time we come to build, it checks those fresh new dependencies.

.PHONY: all clean

all: prog

	rm -f *.o *.d prog .deps

prog: a.o b.o
	$(CXX) -o prog $^

#We can't use the built in rule for compiling c++.
#Let's say a header changed to include a new header.
#That change would cause a rebuild of the .cc file
#but not the .d file. 
# Arguments mean:
# -MMD output #include dependencies, but ignore system headers, while
#      the build continues
# -MF output dependencies to a file, so other crud goes to the screen
# -MP make a blank target for dependencies. That way if a dependency is
#     deleted, make won't fail to build stuff 
%.o %d:
	$(CXX) -c $< $(CXXFLAGS) -MMD -MP -MF $*.d

#Search for all .d files and add them
#Note that .d files are built in with .o files. If there's no .o yet,
#then there's no .d. But that's OK; there are no precise dependencies,
#but we know we have to rebuild it anyway because the .o is not there.
#So, it's OK to only pick up .d files that have alredy been generted.
include $(wildcard *.d)

And the full set of files.

Now try creating an empty “c.h” and make b.h include it. Type make (you’ll see a.o rebuild). Now touch “c.h”. Make will rebuild a.o again as it should. Revert b.h and remove c.h. Make again builds it correctly. Touch c.h again and type make and nothing (correctly) gets rebuilt.

Actually, the mechanism is is subtle and interesting. It looks like a simple include directive, but the key is that if morestuff doesn’t exist, make will build it first, include it and then process the rest of the build. You cn do cool stuff there, but fortunately it’s not needed here.

A simple custom iterator example

Someone recently asked me about custom iterators in C++. It turns out most of the examples and tutorials are not ideal, especially if you get the principle but not the mechanics.

Without getting bogged down in boilerplate, here’s a complete standalone example which compiles and runs. It provides an iterator that iterates over a string one chunk at a time, returning those chunks (words and gaps):

#include <string>
#include <cctype>  //for isgraph
#include <iterator>
#include <iostream>
#include <algorithm>
#include <vector>

using namespace std;

//Iterate over a string one chunk at a time.
//Deriving from std::iterator is necessary to make the code work with 
//the standard algorithms (for now). Either that or specialising the
//std::iterator_traits type. We're telling it that we have a forward 
//iterator and that the value returned from dereferencing is a string.
//Note that this derivation is not necessary for simple for-loops or
//using the iterator in any custom context.
class ChunkStringIterator: public iterator<forward_iterator_tag, string>
	//Represent the string remaining to be iterated over as a
	//string iterator pair.
	string::const_iterator b, e;
		ChunkStringIterator(string::const_iterator b_, string::const_iterator e_)
		:b(b_), e(e_)
		//Return the chunk starting at the current position
		string operator*()
			auto n = find_next();
			return string(b, n);
		//Advance to the next chunk
		ChunkStringIterator& operator++()
			b = find_next();
			return *this;

		//We need this operator for the classic iteration loop
		bool operator!=(const ChunkStringIterator& i) const
			return b != i.b;

		//This is the meat: it returns an iterator to one-past-the-end of
		//the current chunk. Chunk is defigned as a contiguous region of 
		//chars of the same category.
		string::const_iterator find_next()
			auto current = b;
			bool printable = isgraph(*current);
			while(current != e && (bool)isgraph(*current) == printable)
			return current;


//This class is a trivial wrapper for string that dishes out
//the custom iterators.
class ChunkString
	const string& data;
		ChunkString(const string& s)
		ChunkStringIterator begin()
			return ChunkStringIterator(data.begin(), data.end());

		ChunkStringIterator end()
			return ChunkStringIterator(data.end(), data.end());

int main()
	string foo = "Hello, world. This is a bunch of chunks.";
	//A for-each loop to prove it works.
	for(auto w: ChunkString(foo))
		cout << "Chunk: -->" << w << "<--\n";

	//Use a standard algorithm to prove it works
	vector<string> chunks;
	copy(ChunkString(foo).begin(), ChunkString(foo).end(), back_inserter(chunks));
	cout << "there are " << chunks.size() << " chunks.\n";


Note: this is not complete. I haven’t implemented anything that’s not necessary to the the for each example working. It’s basically missing the other comparison operator (==) and the other dereference operator (->). Probably some other things too.