1
// Candle-flicker LED program for ATtiny13A based on analysis at[0]
2
// and IIR filter proposed and tested at [1].
4
// Hardware/peripheral usage notes:
5
// LED control output on pin 5.//changed to fit nanjg layout pin 6 PB1
6
// Uses hardware PWM and the PWM timer's overflow to
7
// count out the frames. Keeps the RNG seed in EEPROM.
9
// [0] http://inkofpark.wordpress.com/2013/12/15/candle-flame-flicker/
10
// [1] http://inkofpark.wordpress.com/2013/12/23/arduino-flickering-candle/
12
#include <avr/eeprom.h>
13
#include <avr/interrupt.h>
15
#include <avr/sleep.h>
16
#include <util/delay.h>
21
// simple 32-bit LFSR PRNG with seed stored in EEPROM
22
#define POLY 0xA3AC183Cul
24
#define PWM OCR0B //PWM-value
29
// set up PWM on pin 5 (PORTB bit 0) using TIMER0 (OC0A)
30
#define TIMER0_OVF_RATE 2343.75 // Hz
31
static void init_pwm() {
32
// COM0A = 10 ("normal view")
33
// COM0B = 00 (disabled)
34
// WGM0 = 001 (phase-correct PWM, TOP = 0xff, OCR update at TOP)
35
// CS = 010 (1:8 with IO clock: 2342 Hz PWM if clock is 9.6Mhz )
36
TCCR0A = 0x81; // 1000 0001
37
TCCR0B = 0x02; // 0000 0010
38
//DDRB |= 1 << DDB0; // pin direction = OUT
43
// set up PWM on pin 6 (PORTB bit 1) using TIMER1 (OC0A)
44
#define TIMER0_OVF_RATE 2343.75 // Hz
45
#define pwminit() do{ TCCR0A=0b00100001; TCCR0B=0b00000010; DDRB |= 1 << DDB1;}while(0) //changed to nanjg pinlayout..hopefully correct
48
static uint32_t lfsr = 1;
49
static void init_rand() {
50
static uint32_t EEMEM boot_seed = 1;
52
lfsr = eeprom_read_dword(&boot_seed);
53
// increment at least once, skip 0 if we hit it
54
// (0 is the only truly unacceptable seed)
55
do {lfsr++;} while (!lfsr);
56
eeprom_write_dword(&boot_seed, lfsr);
58
static inline uint8_t rand(uint8_t bits) {
62
for (i = 0; i < bits; i++) {
64
lfsr = (lfsr >> 1) ^ (-(lfsr & 1ul) & POLY);
73
// approximate a normal distribution with mean 0 and std 32.
74
// does so by drawing from a binomial distribution and 'fuzzing' it a bit.
76
// There's some code at [0] calculating the actual distribution of these
77
// values. They fit the intended distribution quite well. There is a
78
// grand total of 3.11% of the probability misallocated, and no more than
79
// 1.6% of the total misallocation affects any single bin. In absolute
80
// terms, no bin is more than 0.05% off the intended priority, and most
81
// are considerably closer.
83
// [0] https://github.com/mokus0/junkbox/blob/master/Haskell/Math/ApproxNormal.hs
84
static int8_t normal() {
85
// n = binomial(16, 0.5): range = 0..15, mean = 8, sd = 2
86
// center = (n - 8) * 16; // shift and expand to range = -128 .. 112, mean = 0, sd = 32
89
for (i = 0; i < 16; i++) {
90
center += rand(1) << 4;
93
// 'fuzz' follows a symmetric triangular distribution with
94
// center 0 and halfwidth 16, so the result (center + fuzz)
95
// is a linear interpolation of the binomial PDF, mod 256.
96
// (integer overflow corresponds to wrapping around, blending
97
// both tails together).
98
int8_t fuzz = (int8_t)(rand(4)) - (int8_t)(rand(4));
104
// fixed-point 2nd-order Butterworth low-pass filter.
105
// The output has mean 0 and std deviation approximately 0.53
106
// (4342 or so in mantissa).
108
// integer types here are annotated with comments of the form
109
// "/* x:y */", where x is the number of significant bits in
110
// the value and y is the base-2 exponent (negated, actually).
112
// The inspiration for the filter (and identification of basic
113
// parameters and comparison with some other filters) was done
114
// by Park Hays and described at [0]. Some variations and
115
// prototypes of the fixed-point version are implemented at [1].
117
// [0] http://inkofpark.wordpress.com/2013/12/23/arduino-flickering-candle/
118
// [1] https://github.com/mokus0/junkbox/blob/master/Haskell/Math/BiQuad.hs
119
#define UPDATE_RATE 60 // Hz
120
#define FILTER_STDDEV 4.3e3
121
static int16_t /* 15:13 */ flicker_filter(int16_t /* 7:5 */ x) {
123
/* 3:3 */ a1 = -7, // round ((-0.87727063) * (1L << 3 ))
124
/* 4:5 */ a2 = 10, // round ( 0.31106039 * (1L << 5 ))
125
/* 7:10 */ b0 = 111, // round ( 0.10844744 * (1L << 10))
126
/* 7:9 */ b1 = 111, // round ( 0.21689488 * (1L << 9 ))
127
/* 7:10 */ b2 = 111; // round ( 0.10844744 * (1L << 10))
132
int16_t /* 15:13 */ y;
134
// take advantage of fact that b's are all equal with the
135
// chosen exponents (this is a general feature of 2nd-order
136
// discrete-time butterworth filters, actually)
137
int16_t /* 15:14 */ bx = b0 * x >> 1;
138
y = (bx >> 1) + (d1 >> 0);
139
d1 = (bx >> 0) - (a1 * (y >> 3)) + (d2 >> 1);
140
d2 = (bx >> 0) - (a2 * (y >> 4));
147
// scipy-designed filter with 4Hz cutoff and 100Hz sample rate.
148
// parameters set using the following little snippet:
150
// >>> cutoff = 2.2/(rate/2.)
151
// >>> b, a = signal.butter(2, cutoff, 'low')
155
// b = 0.00184036, 0.00368073, 0.00184036
156
// a = 1. , -1.87503716, 0.88239861
158
// B was then rescaled to 1, 2, 1 in order to improve the
159
// numerical accuracy by both reducing opportunities for
160
// round-off error and increasing the amount of precision
161
// that can be allocated in other places. (it also
162
// conveniently replaces 3 multiplies with shifts)
164
// This function implements the transposed direct-form 2, using
165
// 16-bit fixed-point math. Integer types here are annotated
166
// with comments of the form "/* x:y */", where x is the number
167
// of significant bits in the value and y is the (negated)
170
// The inspiration for the filter (and identification of basic
171
// parameters and comparison with some other filters) was done
172
// by Park Hays and described at [0]. Some variations and
173
// prototypes of the fixed-point version are implemented at [1].
174
// The specific parameters of this filter, though, are changed
177
// The higher sampling rate puts the nyquist frequency at 50 Hz
178
// (digital butterworth filters appear to roll off to -inf at
179
// the nyquist rate, though the fixed-point implementation has
180
// a hard noise floor due to the fixed quantum). This probably
181
// doesn't make that much difference visually but it seems to
182
// greatly improve the numerical properties of the fixed-point
185
// SciPy's butterworth parameter designer also seems to put the
186
// cutoff frequency higher than advertised (or maybe that's an
187
// effect of the limited sampling rate, but PSD plots seem to
188
// show the cutoff much higher than expected, and an 8Hz cutoff
189
// was just far too fast-moving for my taste), so I pulled it
190
// down to something that looked pleasant when running, yielding
191
// the following filter.
193
// [0] http://inkofpark.wordpress.com/2013/12/23/arduino-flickering-candle/
194
// [1] https://github.com/mokus0/junkbox/blob/master/Haskell/Math/BiQuad.hs
195
#define UPDATE_RATE 156.25 // Hz
196
#define FILTER_STDDEV 6.2e3
197
static int16_t /* 15:13 */ flicker_filter(int8_t /* 7:5 */ x) {
199
/* 9:8 */ a1 = -480, // round ((-1.87503716) * (1L << 8 ))
200
/* 8:8 */ a2 = 226; // round ( 0.88239861 * (1L << 8 ))
201
// b0 = 1, b1 = 2, b2 = 1; multiplies as shifts below
206
int16_t /* 15:6 */ y;
208
y = ((int16_t) x << 1) + (d1 >> 0);
209
d1 = ((int16_t) x << 2) - (a1 * (int32_t) y >> 8) + (d2 >> 0);
210
d2 = ((int16_t) x << 1) - (a2 * (int32_t) y >> 8);
217
static uint8_t next_intensity() {
218
const uint8_t m = 171, s = 2;
219
const int16_t scale = (s * FILTER_STDDEV) / (255 - m);
221
int16_t x = m + flicker_filter(normal()) / scale;
222
return x < 0 ? 0 : x > 255 ? 255 : x;
225
// use timer so update rate isn't
226
// affected by calculation time
227
static volatile bool tick = false;
229
static uint8_t cycles = 0;
232
cycles = TIMER0_OVF_RATE / UPDATE_RATE;
239
//init_pwm(); //old standard init
240
pwminit(); //changed to nanjg standard PB1 leg 6
243
//enable timer overflow interrupt timer
249
// sleep till next update, then perform it.
250
while (!tick) sleep_mode();
253
PWM = next_intensity();