Hexagonal Blanket Dimensions


Look at the right-most figure in diagram. The small triangle defines the all the major dimensions of the hexagon. Assuming the user measures the edge-to-edge dimension $c$, he or she can calculate the rest of the measurements. A simple right-triangle expression gives the relationship between $s/2$ and $c/2$, and is readily solved for $s$ in terms of $c$,
$$ \frac{s}{2} = \frac{c}{2}\tan\left(30°\right) $$
$$  s = c\tan\left(30°\right) = c\frac{\sqrt{3}}{3}. $$

Similarly, the Pythagorean formula gives the relationship between $c$ and the two other dimensions,
$$  d = \sqrt{ c^2 + s^2} = \frac{2c}{\sqrt{3}}.$$

Simplifying this equation by substituting for $s$ into the previous equation and simplifying produces
$$  \text{long pitch} = s + \frac{d-s}{2} = \frac{s+d}{2}  = \frac{\sqrt{3}}{2}c \approx 0.87c. $$

A blanket $n$ hexagons by $m$ hexagons will be approximately
0.87c\,n \times m\,c, $$
where the $m$ and $n$ dimensions are as shown in the next figure.


If a blanket will be wider at the ends, as shown in the figure, then $n$ will be odd. The total number of hexagons will then be
$$ \text{number of hexagons} = m \frac{n+1}{2} + (m-1)\frac{n-1}{2}. $$

In the example there are $(n+1)/2=5$ tall columns and $(n-1)/2=4$ short columns. So there are 6×5=30 hexagons in tall columns and (6-1)×4=20 hexagons in the short columns, or 50 hexagons overall. Using equation for blanket size and assuming $c=10$ inches, the approximate dimensions of the finished blanket are 9×10×0.87 by 6×10, or 78.3 by 60 inches.

Container Gardening

Container gardening is ideal for the small suburban backyard without much space. It is even better suited to the arid climate because the required water is much less—in my system I water less than twice a week and use less than a gallon. I’m using a container gardening system worked out by others, and it is really impressive. My father put me onto good video documentation as well as EarthTainer’s excellent guide for soil composition. The video discusses using aluminized bubble wrap, sold as insulation. I did not initially wrap my buckets in insulation. An experiment is in order.


My tomatoes are dying. The leaves have started drying, browning, and curling. Older leaves are failing first. My diagnosis is an untreatable fungal infection called Verticillium wilt fungus or Fusarium wilt fungus. The symptoms, according to gardeningknowhow.com, include worsening wilting after watering. My cherry tomato has produced three fruits, none of which are likely to be eatable. Look carefully, you can see a hornworm with more hope than promise.


Apparently soil temperature can influence the spread of Fusarium. Temperatures above 30 C (86 F) contribute to the spread.

I put a thermal probe in each of my two containers for a day. After a day I wrapped one of them in insulation and placed white foam core on top. The buckets are in different positions, so I’m not looking at the difference from one bucket to the next, but rather the change in the difference. I plotted the soil temperature in each bucket over the course of the experiment, along with the air temperature.


The result shows that the insulation helps quite a lot. It lowers the peak container temperature by 3 to 4 degrees F. More importantly, it kept the soil temperature below 80F despite higher air temperature; that should help control fungal growth. In the New Mexico High Desert the insulation should not be considered optional. Partial shade may also be beneficial.

3D Printing the Mountains


Hanging on the wall, under the low ceiling in the basement of my childhood home was a map. It was plastic, and it had a 3D texture of the Colorado Rocky Mountains. Flying over the range in my imagination I could see the high mountain valleys bowled in by the snowcapped peaks.

As an adult I want a map of the Sandia Mountains that I watch from my window. Mountains with all the moods of an impressionist painter. Sunshine illuminating the mountain, standing proud before a black cumulus mountain range. With a 3D printer, and the wonderful availability of open source software, I can. Here is the recipe.

First acquire digital elevation model (DEM) data from the source of your choice. I got my data from viewfinderpanoramas.org, you have to drill down to the coverage map. I like these data, at 3 arcsecond resolution each pixel represents roughly a 30 meter square. Most publically available DEM data is coarser, which is fine if you are printing large areas.

Download the data, unzip it, and add the data to your QGIS session. With a text layer, I loaded a reference location so I could select the printable region. A screenshot of my loaded map appears below.


Certainly the region you want to print is smaller than the DEM file tiles. To extract a region of the DEM file, use the Raster > Extraction > Clipper tool. Choose the rectangular region, draw what you want, and save the new file as a GeoTIFF. Load the new GeoTIFF into QGIS. Scale the vertical range of the new layer so that all values are between 0 and 255. Go to Raster > Raster Calculator and then use the layer’s minimum and maximum values to apply a formula

255*(value – minimum)/(maximum – minimum)

The layer summary shows the minimum and maximum values. In this screencap the values range between 1648 and 3091.


Convert the scaled GeoTIFF to PNG with Raster > Conversion > Translate. OpenSCAD can import PNGs and generate the STL file that most slicing programs want. I needed a recent version of OpenSCAD to load PNGs successfully, the following screenshot was made with 2015-03.1.


The OpenSCAD code is relatively straightforward, the only difficult part is getting the scales correct. You need the pixel extents of the image and the geographic extent (km) of the print area.

// OpenSCAD file to scale and display a raster

vert_min = 1.609; // km
vert_max = 3.255; // km
horiz_x_extent_km = 14.335; // km
horiz_x_extent_px = 154; // 
horiz_y_extent_px = 150; // pixels

h_factor = (1/horiz_x_extent_px)*horiz_x_extent_km;//puts it in km
v_factor = (vert_max - vert_min)/100;

// To keep the object manifold I had to impose a lower layer, this
// also provides strength for the thinnest parts of the plot.
	  horiz_x_extent_km * horiz_y_extent_px/horiz_x_extent_px,
scale([h_factor, h_factor, v_factor]){
    surface(file = "myrasterclipped4.png", 
            center = true);

From OpenSCAD you can render and export an STL. I print with the excellent Repetier Host software. Repetier loads the STL to produce something like the following screenshot.


And finally 3D print it.


This image is the Sandia Mountains crest, my house is located approximately at the lower right of the print. Each day I look out my back window at the this scene.

My Daughter’s Nightstand

This week I finished a cabinet carpentry project started after Thanksgiving in 2014. It was a stretch of my skills, the first time building a project with full panel construction, where all the panels were solid wood. Except for the bottom of the drawer, the entire structure is made of solid wood. To control cost, the choice of wood was pine and fir. The legs and rails are made from fir, bought as a 2×8 and re-milled by me with a tablesaw and a thickness planer. The top, panels, drawer sides, and the rest are made from 1×6 lumber.


I finished the wood with fourteen thousand coats of milk paint, sanding between each one. The milk paint, in cream color, was top-coated with Danish oil finish because milk paint’s rough finish can stain easily. The green glass knobs came from a collection of my grandfather’s.


The early design planned for through tenons, to leave little exposed pegs. I had also intended to carve a sun-and-moon theme, later to be painted with bright colors.


The panels were cut with the tablesaw. I mounted the panels in the tenoning jig. The runout was ridiculous. In spite of using a micrometer to true the tenoning jig, I had a very difficult time keeping the panel thickness constant. Next time I want to make panels I’ll buy a router jig for the purpose.

The drawer slides are also custom. I had roller suspended drawer slides, but neglected to design the drawer with enough clearance. The drawer front attaches to the drawer sides with hand-cut half-blind dovetails. I was pleased that my first attempt at half-blind dovetails worked out well.


The carcase is entirely mortise and tenon joints. I cut the tenons on the tablesaw wit the jig. In the entire cabinet the only metal fasteners are the clips that hold the table tops on, and the L-brackets that hold the draw slides on.


3D-Printer and Your Coffee Grinder


Two weeks ago I acquired a 3D printer; specifically the Printrbot Simple Metal Kit and, after some warranty help, it’s printing well. The first thing I printed was the fan shroud for the 3D printer itself, as recommended by Printrbot, but the next things were for my own design. I printed a set of three sieves, to fulfill my long-time dream of quantifying the performance of coffee mills.

For those of you who read my coffee blog posts years ago, I was frustrated because nobody quantifies the grind performance. Vendors and coffee pundits are happy to talk about the merits of a conical burr grinder or fret about the cheap blade grinder you got from your lost year with Gevalia. Quantifying grinder performance should be pretty straightforward. Take a set of sieves and sort the grinds by size, the more consistent grinders will produce more grounds in a narrow range of sizes. Spoiler: I have only measured my whirling blade grinder so far, not my Capresso burr grinder or any of the commercial grinders.

The printed sieves are just cylinders with a printed mesh bottom. They don’t stack well, since I used too little taper. My finest one had holes so small that it plugged with the finest coffee dust, so there is more to do in the sieve design. Nevertheless, I started with the mesh from Thingiverse, and added 15 mm of wall height. When I get a version of these that stacks I’ll post the design to the Thingiverse too.


I arranged the sieves in a stack and shook. The sieves clogged almost immediately, so I took a small brush and worked the from the top layer down until I had good separation.


The results, in the following picture, show that almost no particles were larger than my largest mesh, less than half were larger than my medium mesh, and the rest were larger than my fine mesh.


So what are the mesh sizes? I took macro photos of the meshes and then measured them optically—you know, counted pixels. The composite photo below shows the basic idea, and below that is a zoomed-in version of the medium sieve.


The medium sieve is made from “threads”, where each “thread” is two passes with the extruder head. It should be possible to do a single extruder path, but I have not yet tuned the OpenSCAD file to get a consistent result.


Mesh Cell Diagonal Coffee Percent
Coarse 2.5 mm <0.1 g 0%
Medium 1.4 mm 1.0 g 28%
Fine 0.48 mm 2.6 g 72%
  fall through <0.1 g 0%

Whether grind consistency can actually be identified by a taster in a blind test is an open, and wonderful, question. Happy brewing!

The Giving Tree and An Arduino Clock

When we bought our house we hired an electrician to install ceiling fans and lights. He worked, worked, and worked. My illusions of electrician’s work dissolved. Their work is a great deal more like laying bricks than solving equations. Oddly, the truth for electricians is also the truth for electronics.

Nearly two years after starting, more than 18~months after finishing the electronics, source code, and cardboard mock-up I finally finished my clock. I have recreated—nay, improved upon—the greatest alarm clock I have ever owned.


In some ways my clock was started about thirty years ago, when some foresighted person planted a lovely little apple tree in what is now my back yard. The lovely little apple tree was a gorgeous mature giant shady apple tree until the spring of 2012. The spring came warm and early but also trisected by two severe cold spells. The frost killed several of my trees major limbs, and may eventually kill the tree.

One of the major limbs became firewood. Then one of the firewood logs became boards, and finally one of the boards became the main faces of the clock. Making a board from firewood is, no doubt, ancient. For me it was a new experience with old techniques. I used a hand powered bow saw, affixed the firewood in the vice, and ripped boards by hand. It is a tiresome process and the boards were not machined to the parallel faces you get in dimensioned lumber. No matter, hand planes helped me get one face smooth and flat. I marked a constant thickness, flipped the board over, and planed the uneven face. The finished board was smooth, clean, and beautifully figured.

The original vision for the box of the clock had hand-cut dovetails visible on the front. However, the apple wood was unsupportive. It was brittle, and prone to fracturing on detail peices. It was also almost too small. The cherry boards I selected for the sides were not long enough to dovetail join the front and back.


The final result is lovely to me. It has a lid so the interface is not visible on the bedside, and the lid showcases the figure of the apple wood. The display is a simple red that does not keep me awake at night. The plain box design weighs enough to keep from migrating around the nightstand, and the interface is actually good.


Design and Use

I set out to imitate the functions of my beloved but deceased clock, see Behold! A number pad. I started the design with a set of use cases, Set Time, Set Alarm, Turn on Alarm, and Turn off Alarm. The notes below are from the original design before I bought the first part, with no edits but formatting.

Set Time

User opens decorative lid. User sets mode switch “set time” mode. The display is turned to current time. User enters numbers starting with most significant digit, in 24 hour time format. After entering the first number, only that digit is displayed.

If user enters a mistake, he hits “clear (#),” and the process restarts at “user enters numbers”. Clear: display is blanked. To keep the setting, user leaves the “set time” mode. On error, the time remains unchanged. Error conditions are:

  • User exits "set time" mode without entering a new time at all.
  • User exits "set time" mode without entering a valid time.

Set Alarm

Same as “Set Time,” but using “set alarm” mode.

Turn on Alarm

User turns on alarm. Display flashes time the alarm is set for about 5 seconds. Display resumes normal operation. The “alarm on” LED lights. When the time reaches the set alarm time, a beeper will sound. The beeper will sound until the alarm is shut off.

Turn off Alarm

Alarm LED is turned off. If the alarm is currently beeping, the beeping ceases.

Electronics Design

I ended up with a Sparkfun keypad, and care of Andrew’s diligent work I was able to configure the Arduino microprocessor to use the keypad almost painlessly. I arranged the pinout as follows:

Pin Use
0 not available
1 not available
2 keypad 0
3 keypad 1
4 keypad 2
5 keypad 3
6 keypad 4
7 keypad 5
8 keypad 6
9 alarm on/off switch
10 alarm tone (PWM)
14 (A0)  
15 (A1) display brightness potentiometer
16 (A2)  
17 (A3) mode switch (3-state) alarm/time
18 (A4) i2c – display and realtime clock
19 (A5) i2c – display and realtime clock

Pin A3 was used for analog read to determine the state of a simple resistor network. Switching the mode switch changes the resistors in a voltage divider. I designed this before I learned about the Arduino’s internal pull-up resistor, so the analog electronics are more difficult than I would make them today.

In case anyone wants it, the code is here also.

#include &lt;Wire.h&gt;
#include &quot;Adafruit_LEDBackpack.h&quot;
#include &quot;Adafruit_GFX.h&quot;
#include &lt;Keypad.h&gt;
#include &quot;RTClib.h&quot;

DateTime now;

#define PIN_ALARMTONE 10

#define PIN_BRIGHTNESSPOT A2 // on the left

#define COLON 2
#define PM 8
#define ALARM 4
#define EXTRA 16

int brightness( void){
  // the analog reading is between 0 and 1024, so divide by 64 to get the
  // values between 0 and 15
  return analogRead( PIN_BRIGHTNESSPOT) &gt;&gt; 6; 

class ClockDisp : public Adafruit_7segment {
  uint8_t alarmOn;
  uint8_t pmOn;
  uint8_t colonOn;
  void writeState( void);
  void alarm( bool);
  void pm( bool);
  void colon( bool);
  void toggleColon( void);
  void printTime( DateTime now);
  void blank( void);

void ClockDisp::blank( void){
  print(10000, DEC);

void ClockDisp::writeState( void){
  writeDigitRaw( 2, alarmOn | pmOn | colonOn );

void ClockDisp::alarm( bool newState ){
  if( newState){
    alarmOn = 4;
    alarmOn = 0;

void ClockDisp::pm( bool newState){
  if( newState){
    pmOn = 8;
    pmOn = 0;

void ClockDisp::colon( bool newState){
  if( newState){
    colonOn = 2;
    colonOn = 0;

void ClockDisp::toggleColon( void){
  colonOn ^= 2;

void ClockDisp::printTime( DateTime now){
  //DateTime now = RTC.now();
  int hour = now.hour();
  if( hour &gt; 12){
    pmOn = 8;
    hour %= 12;
  } else{
   pmOn = 0;
   if( hour == 0) hour = 12;
  int decimalTime = hour * 100 + now.minute();
  print( decimalTime);
const byte ROWS = 4; //four rows
const byte COLS = 3; //four columns
//define the cymbols on the buttons of the keypads
char keys[ROWS][COLS] = {
byte rowPins[ROWS] = {7,2,3,5}; //connect to the row pinouts of the keypad
byte colPins[COLS] = {6,8,4}; //connect to the column pinouts of the keypad

//initialize an instance of class NewKeypad
Keypad customKeypad = Keypad( makeKeymap(keys), rowPins, colPins, ROWS, COLS); 

//Adafruit_7segment matrix = Adafruit_7segment();
ClockDisp clockDisp = ClockDisp();
int alarmNum = -1;
bool isAlarming = 0;
bool hasAlarmed = 0;

void setup(){


  if (! RTC.isrunning()) {
    Serial.println(&quot;RTC is NOT running!&quot;);
    // following line sets the RTC to the date &amp; time this 
    // sketch was compiled, except on Mac.
    RTC.adjust(DateTime(__DATE__, __TIME__));
  clockDisp.setBrightness( 10);

void loop() {
  Serial.println(&quot;Top of loop&quot;);
  now = RTC.now();
  Serial.print(&quot;Brightness pin: &quot;);
  Serial.println( analogRead( PIN_BRIGHTNESSPOT));
  clockDisp.setBrightness( brightness());
  // Print the time, and update the colon
  Serial.print(&quot;Time decimal: &quot;); Serial.println( now.hour() * 100 + now.minute());
  clockDisp.printTime( now);
  if( millis()%1000 &lt; 500){
  delay( 20);
  /* The alarm function is actually a little bit complicated. It is
  easy to test if the current time is equal to the alarm time. We 
  want the alarm to start beeping and not to stop until the switch 
  is thrown. However, if we immediately turned the switch back on,
  the current time would still exceed the alarm time, and it would
  start beebing. Therefore, I created the hasAlarmed variable, which
  records that the alarm has gone off. hasAlarmed will reset to FALSE
  when the time is less than the alarm time, indicated that we've come
  around the clock again.
  if( digitalRead( PIN_ALARMSWITCH)){
    clockDisp.alarm( 1);
    if( isAlarming){
      if( millis()%1000 &lt; 500){
        tone( PIN_ALARMTONE, 1000);
        noTone( PIN_ALARMTONE);
    } else if( hasAlarmed){
      if( (now.hour() * 100 + now.minute()) &lt; alarmNum ){
        hasAlarmed = 0;
    } else {
      // See if we need to turn on the alarm
      if((now.hour() * 100 + now.minute()) == alarmNum){
        isAlarming = 1;
  }else if( isAlarming){
    clockDisp.alarm( 0);
    noTone( PIN_ALARMTONE);
    isAlarming = 0;
    hasAlarmed = 1;
   // Contol of setting state switch
   int mode = analogRead( A3);
   Serial.print(&quot;Mode reading: &quot;); Serial.println( mode);
   if( mode &gt; 900 ){ // set the time
    Serial.println(&quot;Setting the time&quot;);
    // blank the display  
    int numDigits = 0;
    int timeNum = 0;
    while( analogRead(A3) &gt; 900){
      // wait for keys, display them, check state
      // wait for entry of time
      char customKey = customKeypad.getKey();
      if( customKey &gt; 0){
        if( numDigits &gt; 0){
          timeNum = timeNum*10 + (int)(customKey - '0');
          timeNum = (int)(customKey - '0');
        numDigits ++;  
        clockDisp.print( timeNum);
    // Test to see if we got a valid time:
    if( (timeNum &lt; 2500) &amp;&amp; (timeNum &gt; 0)){
      RTC.adjust( DateTime( now.year(), now.month(), now.day(),
              timeNum / 100, timeNum%100, 0));
      // An error, flash bars at the user
      clockDisp.print( 10000);
      delay( 300);
  }else if( mode &gt; 490){ // state 2
    clockDisp.print( alarmNum);
    int numDigits = 0;
    while( 1 ){
      int ar = analogRead( A3);
      if( (ar &lt; 490) || (ar &gt; 900)) break;

      // wait for keys, display them, check state
      // wait for entry of time
      char customKey = customKeypad.getKey();
      if( customKey &gt; 0){
        if( numDigits &gt; 0){
          alarmNum = alarmNum*10 + (int)(customKey - '0');
          alarmNum = (int)(customKey - '0');
        numDigits ++;  
        //clockDisp.print( (int)(customKey - '0'));
        clockDisp.print( alarmNum);
    // Test to see if we got a valid time:
    if( (alarmNum &gt;= 2500) || (alarmNum &lt; 0)){
      // Invalid time entered, flash error bars at user
      clockDisp.print( 10000);
      delay( 300);
    } else{
      hasAlarmed = 0;

int getDecimalTime(){
  DateTime now = RTC.now();
  int decimalTime = now.hour() * 100 + now.minute();
  return decimalTime;

Wiring an Arduino to a Five-state Guitar Switch

I am rebuilding my Arduino-based data collection system. Named the RIMU for some long-forgotten and pointless acronym. My old one used a momentary switch—a button—to change which sensor is displayed on the LCD. It is awful because if the Arduino is busy when you push the switch then nothing happens. I’m fixing it with a 5-position switch sold for changing the pickup combinations on electric guitars.

I was thinking about some really good machine interfaces, and some really bad. I walked about the house with a camera. Below are some photo collections of good and bad interfaces around the house.


Good interfaces control devices unambiguously with touch. You should not have to look at the device to know what mode you have changed it to. Really great interfaces, like a light switch, reveal their state by touch.

My switch is almost as good as a light switch. With five states it is hard to tell which state it is in by feeling its position; however, it is easy to change state. Even better, its position is an absolute indication of its state. There is no issue of the Arduino being busy when the state changes, because the state is fixed by the switch. It is not toggled by some transient state.

I bought my switch from Sparkfun, and also lifted this product image from them.


If you can find old rotary switches with detents you can use the analog read function of the Arduino to get all the states with a single pin. In this case, though, I will use three pins to represent the states. I believe it is possible to configure this switch with a  few resistors to represent all its states with a single analog input, but it is more complex to figure out and I am not short on pins. My definition of pins. The ones with wires are the only ones I need to use.


With just three digital inputs, the Arduino can read all five switch states. The Arduino’s microprocessor includes built-in pull-up resistors. These are easy to use, but a little bit confusing. A pull-up resistor means that if the pin is not attached to anything, the Arduino will read the pin high.

In the following table I show a Y to indicate how the switch is connected internally. However the ABC column shows what you would read from the pins in the Arduino.

Pos 1 2 3 ABC
1 Y 110 6
2 Y Y 100 4
3 Y 101 5
4 Y Y 001 1
5 Y 011 3
// Define the spins to read the switch
#define pinSwitch1 2 
#define pinSwitch2 3 
#define pinSwitch3 4

unsigned char switchState( void){
  unsigned char state=0, in;
  in = digitalRead( pinSwitch1);
  state |= in &lt;&lt; 2;
  in = digitalRead( pinSwitch2);
  state |= in &lt;&lt; 1;
  in = digitalRead( pinSwitch3);
  state |= in;
  switch( state){
    case 6:
    return 1;
    case 4:
    return 2;
    case 5:
    return 3;
    case 1:
    return 4;
    case 3:
    return 5;
  return 0;
void setup( void){
  pinMode( pinSwitch1, INPUT);
  pinMode( pinSwitch2, INPUT);
  pinMode( pinSwitch3, INPUT);
  digitalWrite( pinSwitch1, HIGH);
  digitalWrite( pinSwitch2, HIGH);
  digitalWrite( pinSwitch3, HIGH);
  Serial.begin( 57600);
  Serial.println( 'Started');

void loop( void){
  Serial.print( &quot;Switch: &quot;);
  Serial.println( switchState());
  Serial.print( &quot; &quot;);
  Serial.print( digitalRead(pinSwitch3));
  Serial.print( &quot; &quot;);
  Serial.print( digitalRead( pinSwitch2));
  Serial.print( &quot; &quot;);
  Serial.print( digitalRead( pinSwitch1));
  Serial.print( &quot; &quot;);
  delay( 1000);


In the desert it almost never rains. I know, you know. When I lived in upstate New York, in my foolish youth, I bought a motorcycle in the beginning of April. I was anxious to learn to ride. Too bad, because it rained every day for an entire month. Now I live in New Mexico. Two years ago we had not one drop between the end of December and July.

Where you live the weather turns with a unique step. Here, the summer breaks in July when the monsoon comes. Moist air from the Gulf of California and the Gulf of Mexico flow north over the state and give rise to afternoon thunderstorms. It is the most beautiful weather of the year here. The clouds rise miles—literally—into the sky and continuously billow in fractal glory.

Last week I configured my Raspberry Pi computer with its PiCam to take time-lapse video out my back window. The view is to the east over the Sandia Mountains. I took pictures about every six seconds from 10:30 am until dark, around 8:30 pm. The whole day compressed into four and a half minutes. The best video was from July 14, others are below.

It is hard to understand the desert if you have never lived in one. In the picture below you can see what it would be like if Leeds were in Albuquerque. A nation that ruled the world for a few hundred years fits comfortably in the desert southwest.composite

The rest of the videos I created are here, with the most interesting at the top.

July 11

July 13

July 8

July 7

July 10

July 16

July 15

Arduino Flickering Candle

Update December 24, 2013: Mokus refined his code so that the distribution is now well-behaved (nearly normal) and the PSD no longer turns up at high frequencies). The plots and post have been updated to reflect this change. He will push code to the same link as available.

In my previous post on Candle Flame Flicker I describe the statistics of the optical intensity flicker from a candle in terms of the probability density of the brightness samples and in terms of the power spectral density of the brightness. In this article I discuss how to make an Arduino drive an LED to flicker in a way that is statistically similar to the measurements I took on a real candle.

In the measurements I observed the spectral roll-off of the candle to start between about 5 and 8 Hz, and to decline at a nominal rate of around 44 dB for each decade increase in frequency. The 2nd-order infinite impulse response filter is computationally efficient in terms of using a small amount of memory and requiring few calculations. However, the Arduino is not good at floating point arithmetic. On the Arduino, floating point is done in software, has relatively few bits of precision, and is about 4 to 40 times slower than fixed point (integer) math. It is quite difficult to find useful benchmarks. The basic process is to create white noise and then filter it to the correct spectral shape. Afterward, the signal is scaled to have the appropriate variance and offset with a constant to represent the average brightness of the “flame”.

The approach I used was to design the IIR in Python with scipy’s signal module. I specified a 2nd order lowpass Butterworth filter, specifying a cutoff frequency of 8 Hz, and a sample frequency of 60 Hz. I normalized the coefficients to work in a 16 bit integer space, following Randy Yates’ 2010 Practical Considerations in FIR Fixed Filter Implementations, mainly. From a synthesis perspective, there is some prior art. Philip Ching, a student at Cornell synthesized candle noise quite cleverly, though he neither reported nor replicated the correct statistics. A fellow with the handle Mokus did a very, very tiny implementation for a microcontroller with only 64 bytes of RAM. He helped me modify his code so I could compare his statistics, and after adjustment his spectrum and distribution match fairly well. The high-frequency of his PSD looks a little different from the other methods, but these may not be noticeable to the observer. Finally, there was Eric Evenchick’s coincidental post on hackaday. Mokus reported that Evanchick’s implementation has too slow an update rate; I noticed that Evanchick did not report on the statistics he was targeting nor what he achieved. I did not recreate his work to test.

Then, on to the tests. I really was interested in building and comparing statistics from a 16 bit implementation, a 32 bit implementation in both a Direct Form I and a Direct Form II implementation. Indeed, I had great difficulty getting the filters to perform because I kept misdesigning the integer coefficients and overflowing the registers. As I sought a solution to the 2nd-order filter approach, I also created a 4-stage digital equivalent of an RC filter. The 4-stage RC approach was always stable though it needed a higher sample rate and used much more processor power to create statistically consistent results. On the other hand, it has a more accurate spectrum. A comparison of three different 16-bit methods to Mokus’ and to the actual measurements is shown in the figure below. The legend shows the mean, standard deviation, and their ratio to the right of the label. The All my filters did a credible job of reconstructing the histogram.


The power spectral density (PSD) of the different methods tells a different story. The Direct Form II 16 bit filter is the most visually appealing of the methods I tried. It rolls off more like the physical data than the other methods, except compared to the 4-stage RC filter. The Direct Form II filter is more computationally efficient.


The results for the 32-bit versions show greater variance than the 16-bit versions, but the quality is not markedly better.




I wrote a proof code for the Arduino UNO both to see it flicker and to test the processor speed—separate parts of the code. The results are that compiling with 1.0.3 resulted in a 4,722 byte program that calculated 10,000 new values in 6,292 ms, or 629 microseconds per value. In theory this could produce a sample rate of nearly 1.6 KHz. Or another way of thinking about this is that the final code uses about 629 us/17 ms or about 4% of the processor capability of the Arduino UNO. That leaves a lot of resources available to do other Arduino work or maybe means it can fit in a cheaper processor.

I have released two pieces of code under the GNU Public License 3, you can get the Python I used for filter design and the Arduino test code at the links. If you want the data, please contact me through the comments and I am willing to provide it.

Candle Flame Flicker

For a project I wanted to make an LED flicker like a candle. I searched for the signal statistics of candle flicker, and I found no data. One student web site suggests that candle flame flicker is a 1/f-type random signal with roll-off of 20 dB per decade increase in frequency. Similar processes are typical for turbulence, so this student’s plan seemed reasonable. However, the student did not discuss whether the signal is Gaussian or not, and did not describe the low-frequency characteristics of the signal. 1/f noises may have a spectrum that follows the 1/f curve to very low frequencies, but because they would require infinite power at f=0 they always have some lower frequency change. I had no data.

Candles aren’t hard to get, and I already had a silicon photodiode for an absorptive smoke density measurement system I’m working on. I also had an analog-to-digital converter (ADC) in an ADS1015 on a breakout board from Adafruit. I made the very simple circuit shown below and attached it to a Raspberry Pi to sample the data. There are a lot more details, but those are later.


The voltage measured by the ADC is directly proportional to the current through the resistor. The current through the reverse-biased diode is directly proportional to the incident optical power. Very simple. I recorded a minute of data at a low 250 Hz sample rate, and subjected the data to analysis. The setup was a nominally dark room with the sensor a few inches from the flame. I flapped my hand at the candle to get it to flicker while recording.




The first graph below is a histogram of sample values.


The histogram is about normally distributed, but the right-hand tail is not Gaussian; it is too fat. In other words the candle is occasionally much brighter than normal. The time series (below) shows the same properties. There are clearly visible large excursions.


For human eyes the candle’s flicker will be well represented if the frequency spectrum of the flicker and the distribution of the flicker approximately match a natural source. The spectrum below can be thought of as an average brightness (the peak at the left), a flat spectral region out to about 4 Hz, and then a 1/f-style roll off at 44 dB per decade.


Numerically, we have the recipe for a flickering candle:

  • Samples should be normally distributed with a standard deviation equal to 0.25 of the mean.
  • The power spectral density of the signal should roll off at about 40 dB per decade with a 3 dB cutoff frequency around six cycles per second.

To make this work on the Arduino using the pulse-width modulated outputs, we can further constrain the problem:

  • The maximum value cannot exceed 255 counts—or make the limit that the mean plus two standard deviations is 255.
  • From this, we can derive that 2s+m = 255, use the fact that 0.25m=s to find that the mean m=170 and the standard deviation is about 42.
  • A sample rate of between 30 and 120 Hz should be more than adequate to satisfy the Nyquist criterion for human vision (see Wikipedia).
  • Values may not go below zero or above 255
  • A second-order infinite impulse response (IIR) filter has a roll-off of 40 dB per decade

If we can find a numerically efficient way to generate a time-series of Gaussian random variables inside the Arduino, filter them with a 2nd-order fixed-point IIR, scale them (if needed) then we should be able to make a flickering candle.

Unfortunately, I have failed repeatedly to get a stable fixed point IIR filter. The cookbook solution specification above should be easy to implement, but I have not found it so. A better solution will have to wait for another post.