DIY extremly Sensitive and cheap Arduino Seismometer

1//
2// nerdaqII  version 0.30  6 September 2012  Martin L. Smith
3//
4// This code replaces bsudaq and includes additional filtering capability
5// documented below. 
6//
7//
8// I Oversampling a/d
9//
10// for a 16 MHz clock, the adc rate is 9615.4 sps and the 512-stack is
11// 18.78 sps.  the final sample rate is our only assurance that we
12// haven't skipped interrupts.
13//
14// we actually use a 1024-long window stepped at 512.  using the
15// double window should (i think) put the rectangular window's first
16// spectral zero at the nyquist frequency for the final sample rate.
17//
18// 24 July 2011  went to 4 x 512 windows.  the aggregate window is 2048 points
19//   and is stepped at 512.  the first zero should be at 18.78/4 or about
20//   4.7 Hz.
21//
22// worth noting that we use the 512 window to set the output sample rate
23// while the number of 512s that are combined scales the spectral response.
24//
25//
26// II Denoising Filter (N)
27//
28// N is a low-pass filter with a transition zone of 0.25 to 2.5Hz.
29// It's function is simply to remove some high-frequency electronic noise
30// without introducing objectionable ringing..
31// We hope it won't substantially alter seismic signals.
32// It's implemented as a 16-tap minimum-phase FIR with a design stopband
33// attenuation of 30 dB and an achieved attenuation of 40+ dB.
34// This is not the filter used in bsudaq.
35//
36//
37// III Detrending Filter (B)
38//
39// B is a 2 pole IIR butterworth high-pass with a 30 second cutoff.  It's
40// function is to remove baseline drift and ameliorate step offsets.  It
41// should not alter seismic signals.
42//
43//
44// IV Long-period Boost Filter (L)
45//
46// L is an second-order butterworth band-pass filter with cutoffs of 0.1 and
47// 0.05 Hz and a stopband attenuation of 60 dB.  It's function is to boost signals
48// at periods of 5-20 seconds by an amount set by the adjustable
49// boost factor, BfdB, below.
50//
51//
52// V Boost Factor, BfdB
53//
54// BfdB is a scale factor, specified in dB, by which the output of the
55// long-period boost filter, L, is multiplied before be combined with the
56// non-boosted output.  It's value roughly determines how much longer-period
57// surface waves are boosted before being added to the body wave data.
58// It's specified in dB so BfdB = 0 corresponds to a multiplier of 1.0 and
59// BfdB = 20 corresponds to a mutiplier of 10.0.
60//
61
62#include <avr/io.h>
63#include <avr/interrupt.h>
64
65
66// debugging support
67
68// add free memory info to output stream
69#include "MemoryFree.h"
70#define FREEMEMCHECK
71#undef FREEMEMCHECK
72
73// allow/suppress detrend filter
74const int do_detrend = 1;
75
76
77// arduino definitions
78
79const int ledPin = 13;
80
81// global variables for oversampling and the running window
82
83volatile unsigned long runningsum;
84volatile unsigned short runningcount;
85volatile unsigned long prev3;
86volatile unsigned long prev2;
87volatile unsigned long prev1;
88volatile unsigned long current_sum;
89
90volatile boolean next_sample_ready;
91volatile unsigned int next_sample;
92
93//
94// coefficients for filter N
95//
96extern "C" {
97#include "denoise_filter.h"
98};
99
100const unsigned int ncoeff = sizeof(coeff) / sizeof(coeff[0]);
101float lagarray[ncoeff];
102unsigned short first_adc = 10; // initial samples to skip
103unsigned short first_loop = 1; // deal with the first composite sample
104
105
106// load defs for the two iir biquad sets
107
108extern "C" {
109#include "biquad.h"
110#include "boost_filter.h"
111#include "detrend_filter.h"
112};
113
114biquad_z_t boost_z[BOOST_BIQUADS_SIZE];
115biquad_z_t detrend_z[DETREND_BIQUADS_SIZE];
116
117
118// the boost factor and filter mode: not const so we can alter them at
119// runtime though at the moment we just keep them fixed.
120
121float BfdB = 20.0;
122int filtermode = 4;
123float Bfmult;
124
125
126void initADC() {
127
128  // don't depend on global initializers
129  runningsum = 0;
130  runningcount = 0;
131  prev1 = prev2 = prev3 = 0;
132  current_sum = 0;
133
134  next_sample_ready = false;
135  next_sample = 0;
136
137  // internal AVcc ref, no left adj, channel 0
138  ADMUX = _BV(REFS0);
139    
140  // ACME off, free-running mode
141  ADCSRB = 0;
142
143  // enable ad, intr, auto
144  // set divisor to 128
145  // start first conversion
146  ADCSRA = _BV(ADEN) | _BV(ADIE) | _BV(ADATE)
147    | _BV(ADPS0) | _BV(ADPS1) | _BV(ADPS2)
148    | _BV(ADSC);
149}
150
151
152ISR(ADC_vect) {
153  byte low = ADCL;
154  byte high = ADCH;
155  if(first_adc) {
156    --first_adc;
157    return;
158  }
159  runningsum += low + (high << 8);
160  if(++runningcount == 512) {
161    if(first_loop)
162      prev3 = prev2 = prev1 = runningsum;
163    else {
164      prev3 = prev2;
165      prev2 = prev1;
166      prev1 = current_sum;
167    }
168    current_sum = runningsum;
169    // if we sum two we have to unshift 4
170    // if we sum four we have to unshift 5
171    next_sample =
172      (unsigned int) ((current_sum + prev1 + prev2 + prev3) >> 5);
173    runningsum = 0;
174    runningcount = 0;
175    next_sample_ready = true;
176  }
177}
178
179
180void setup() {
181  unsigned int i;
182  Serial.begin(9600);
183  Serial.flush();
184  first_loop = 1;
185  first_adc = 10;
186  initADC();
187  pinMode(ledPin, OUTPUT);
188  digitalWrite(ledPin, HIGH);
189}
190
191
192float process_sample(const float y) {
193    float z;
194    float fL;
195    unsigned int i;
196
197    if(filtermode == 1)
198        return y;
199
200    if(first_loop) {
201
202      first_loop = 0;
203      
204        // flood the lag array with the first value
205        for(i = 1; i < ncoeff; i++)
206            lagarray[i] = y;
207
208        // initialize the biquad delay lines
209        biquad_clear(detrend_z, DETREND_BIQUADS_SIZE, (biquad_sample_t) y);
210        biquad_clear(boost_z, BOOST_BIQUADS_SIZE, (biquad_sample_t) y);
211
212        // compute the multiplicative gain factor
213        Bfmult = pow(10.0, BfdB / 20.0);
214
215    } else {
216      
217        // update the bucket brigade
218        for(i = ncoeff - 1; i > 0; i--)
219            lagarray[i] = lagarray[i - 1];
220	lagarray[0] = y;
221
222    }
223
224    // apply N to the raw sample series
225    z = 0.0;
226    for(i = 0; i < ncoeff; i++)
227        z += lagarray[i] * coeff[i];
228
229    if(filtermode == 2)
230        return z;
231
232    // apply B to the output of N (if allowed)
233    if(do_detrend)
234        z = biquad_filter(z, detrend_biquads, detrend_z, DETREND_BIQUADS_SIZE)
235            * detrend_biquads_g;
236
237    if(filtermode == 3)
238        return z;
239
240    // apply L to the output of BN
241    fL = biquad_filter(z, boost_biquads, boost_z, BOOST_BIQUADS_SIZE)
242        * boost_biquads_g;
243    // compute and return the weighted trace
244    return z + Bfmult * fL;
245}
246
247
248// we drive the led pin (13) high during inter-sample idle times.  i hope that
249// we'll be able to tell if we have enough cpu overhead by looking to see if
250// the led is flashing at 18.78 Hz.
251
252// the halfscale correction adjusts the B filter output to be 32768
253// instead of 0.
254
255// for version 0.1 of the code, freeMemory() returns 1319 (out of a
256// maximum of 2000).
257
258void loop() {
259  unsigned int filtered_signal;
260  const float halfscale = 32768.0;
261
262  if(next_sample_ready == false) {
263      digitalWrite(ledPin, HIGH);
264      return;
265  }
266  digitalWrite(ledPin, LOW);
267  next_sample_ready = false;
268
269  filtered_signal = (unsigned int) (process_sample(float(next_sample)
270                                                   - halfscale) + halfscale);
271
272  #ifdef FREEMEMCHECK
273  Serial.print(freeMemory(), DEC);
274  Serial.print("  ");
275  #endif
276
277  Serial.println(filtered_signal, DEC);
278  Serial.flush();
279}