~gabe/flashlight-firmware/anduril2

« back to all changes in this revision

Viewing changes to Werner/candleflicker/candleflicker.c

  • Committer: Selene Scriven
  • Date: 2015-02-16 20:01:11 UTC
  • Revision ID: ubuntu@toykeeper.net-20150216200111-my1dk6eqz1mv9cif
Added Werner's candle flicker firmware.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
// Candle-flicker LED program for ATtiny13A based on analysis at[0]
 
2
// and IIR filter proposed and tested at [1].
 
3
// 
 
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.
 
8
// 
 
9
// [0] http://inkofpark.wordpress.com/2013/12/15/candle-flame-flicker/
 
10
// [1] http://inkofpark.wordpress.com/2013/12/23/arduino-flickering-candle/
 
11
#define F_CPU 9600000
 
12
#include <avr/eeprom.h>
 
13
#include <avr/interrupt.h>
 
14
#include <avr/io.h>
 
15
#include <avr/sleep.h>
 
16
#include <util/delay.h>
 
17
 
 
18
#include <stdbool.h>
 
19
#include <stdint.h>
 
20
 
 
21
// simple 32-bit LFSR PRNG with seed stored in EEPROM
 
22
#define POLY 0xA3AC183Cul
 
23
 
 
24
#define PWM OCR0B  //PWM-value 
 
25
 
 
26
 
 
27
 
 
28
 
 
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
 
39
 
 
40
 
 
41
}
 
42
 
 
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
 
46
 
 
47
 
 
48
static uint32_t lfsr = 1;
 
49
static void init_rand() {
 
50
    static uint32_t EEMEM boot_seed = 1;
 
51
    
 
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);
 
57
}
 
58
static inline uint8_t rand(uint8_t bits) {
 
59
    uint8_t x = 0;
 
60
    uint8_t i;
 
61
    
 
62
    for (i = 0; i < bits; i++) {
 
63
        x <<= 1;
 
64
        lfsr = (lfsr >> 1) ^ (-(lfsr & 1ul) & POLY);
 
65
        x |= lfsr & 1;
 
66
    }
 
67
    
 
68
    return x;
 
69
}
 
70
 
 
71
 
 
72
 
 
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.
 
75
// 
 
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.
 
82
// 
 
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
 
87
    int8_t center = -128;
 
88
    uint8_t i;
 
89
    for (i = 0; i < 16; i++) {
 
90
        center += rand(1) << 4;
 
91
    }
 
92
    
 
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));
 
99
    return center + fuzz;
 
100
}
 
101
 
 
102
#ifdef INKOFPARK
 
103
 
 
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).
 
107
//
 
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).
 
111
// 
 
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].
 
116
//
 
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) {
 
122
    const int16_t
 
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))
 
128
    static int16_t
 
129
        /* 15:13 */ d1 = 0,
 
130
        /* 15:14 */ d2 = 0;
 
131
    
 
132
    int16_t /* 15:13 */ y;
 
133
    
 
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));
 
141
    
 
142
    return y;
 
143
}
 
144
 
 
145
#else
 
146
 
 
147
// scipy-designed filter with 4Hz cutoff and 100Hz sample rate.
 
148
// parameters set using the following little snippet:
 
149
// >>> rate=156.25
 
150
// >>> cutoff = 2.2/(rate/2.)
 
151
// >>> b, a = signal.butter(2, cutoff, 'low')
 
152
//
 
153
// which yields:
 
154
//
 
155
//  b = 0.00184036,  0.00368073,  0.00184036
 
156
//  a = 1.        , -1.87503716,  0.88239861
 
157
//
 
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)
 
163
//
 
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)
 
168
// base-2 exponent.
 
169
// 
 
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
 
175
// as follows:
 
176
// 
 
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
 
183
// version.
 
184
// 
 
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.
 
192
//
 
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) {
 
198
    const int16_t
 
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
 
202
    static int16_t
 
203
        /* 15:6 */ d1 = 0,
 
204
        /* 15:6 */ d2 = 0;
 
205
    
 
206
    int16_t /* 15:6 */ y;
 
207
    
 
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);
 
211
    
 
212
    return y;
 
213
}
 
214
 
 
215
#endif
 
216
 
 
217
static uint8_t next_intensity() {
 
218
    const uint8_t m = 171, s = 2;
 
219
    const int16_t scale = (s * FILTER_STDDEV) / (255 - m);
 
220
    
 
221
    int16_t x = m + flicker_filter(normal()) / scale;
 
222
    return x < 0 ? 0 : x > 255 ? 255 : x;
 
223
}
 
224
 
 
225
// use timer so update rate isn't
 
226
// affected by calculation time
 
227
static volatile bool tick = false;
 
228
ISR(TIM0_OVF_vect) {
 
229
    static uint8_t cycles = 0;
 
230
    if (!cycles--) {
 
231
        tick = true;
 
232
        cycles = TIMER0_OVF_RATE / UPDATE_RATE;
 
233
    }
 
234
}
 
235
 
 
236
int main(void)
 
237
{
 
238
    init_rand();
 
239
    //init_pwm(); //old standard init
 
240
        pwminit(); //changed to nanjg standard PB1 leg 6
 
241
    
 
242
    
 
243
 //enable timer overflow interrupt timer
 
244
    TIMSK0 = 1 << TOIE0;
 
245
    sei();
 
246
    
 
247
    while(1)
 
248
    {
 
249
        // sleep till next update, then perform it.
 
250
        while (!tick) sleep_mode();
 
251
        tick = false;
 
252
        
 
253
        PWM = next_intensity();
 
254
                //PWM= 222;
 
255
    }
 
256
}