Seismograph on an Arduino Uno

// Seismograph: an MPU6050 accelerometer, a SW-420 vibration switch and an OLED.
//
//   A4, A5   MPU6050 and OLED on one I2C bus (addresses 0x68 and 0x3C)
//   pin 2    SW-420 DO, the vibration switch: goes LOW when it is shaken
//   pin 13   the Uno's own LED, lit while an earthquake is going on
//
// # What a seismograph does
//
// A seismograph is a very sensitive accelerometer with a pen. It measures how
// hard the ground is being shaken, many times a second, and draws the
// wiggles. The MPU6050 reports acceleration in "g", where 1 g is gravity.
// Standing still on a table it reads 1 g straight down, which is no use for
// finding a shake, so the sketch keeps a slow running average of each axis
// (the part that does not change: gravity and any tilt) and subtracts it. What
// is left is the shaking.
//
// # Intensity is not magnitude
//
// The size of an earthquake at its source, its magnitude (Richter or moment
// magnitude), is the energy it released. It takes several instruments in
// different places to work out. One sensor cannot know it. What one sensor can
// measure is how hard the ground shook right HERE, the peak ground
// acceleration. That is what the screen shows as PEAK, in g, and the
// Modified Mercalli Intensity (MMI) it corresponds to, from I (not felt) to
// X and up (destroyed buildings). The table is the US Geological Survey's
// instrumental intensity: 0.0017 g is about the weakest shake a person
// notices, 0.092 g makes things fall off shelves, over 1 g throws you.
//
// # The traces
//
// Two strips: the top is the sideways (X) shaking, the bottom is forwards and
// back (Y), each with a dotted line at zero. The pen sweeps left to right
// every six seconds, one column at a time, and rubs out the old picture just
// ahead of itself, like the paper drum of a real seismograph. Every column is
// the highest and lowest reading in its 50 milliseconds, drawn as one line, so
// a short sharp jolt is never missed between two columns.
//
// The SW-420 is the burglar alarm of the circuit. It is a spring and a metal
// pin that touch when shaken, which pulls DO low. An interrupt notices that
// even while the sketch is busy drawing, and counts it as a TRIP. A trip, or
// the MPU6050 reading more than 0.02 g, starts an EARTHQUAKE: the sketch
// records the peak and the time until the ground has been quiet for a second
// and a half.
//
// The Serial Plotter shows the same traces live: x, y, z and peak, in g.
//
// Tilt the MPU6050's pitch and roll back and forth to shake the ground.

#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>

const uint8_t MPU = 0x68;
const int TRIP_PIN = 2;
const int LED_PIN = 13;

const float G_PER_COUNT = 1.0 / 16384.0;  // +-2 g range: 16384 counts per g
const float START_G = 0.02;               // shaking harder than this starts an earthquake
const float QUIET_G = 0.008;              // under this is quiet
const unsigned long QUIET_MS = 1500;      // quiet for this long ends it

// The screen: two strips of 22 rows under the text, ten pixels each way from
// the centre line, which is 0.15 g at full scale.
const int CENTRE_X = 31;
const int CENTRE_Y = 52;
const int SPAN = 10;
const float PIXELS_PER_G = SPAN / 0.15;

Adafruit_SSD1306 display(128, 64, &Wire, -1);

volatile uint8_t trips = 0;
volatile unsigned long lastTripAt = 0;
uint8_t seenTrips = 0;

float base[3];          // the slow average: gravity and tilt
bool haveBase = false;
float dyn[3];           // what is left: the shaking
float mag = 0;          // its size, in g

bool quake = false;
float peak = 0;         // the biggest shaking of the last (or current) earthquake, g
unsigned long quakeStart = 0;
unsigned long quakeLength = 0;  // ms, final once it has ended
unsigned long lastLoud = 0;

int column = 0;
int lowest[2];          // the lowest and highest pixel offset seen this column
int highest[2];
int samplesInColumn = 0;
int columnsSinceShown = 0;

unsigned long lastSample = 0;
unsigned long lastPlot = 0;

// The spring in a SW-420 bounces: one knock makes a burst of pulses, so only
// the first pulse in a quarter of a second counts as a trip.
void onTrip() {
  unsigned long now = millis();
  if (now - lastTripAt > 250 || trips == 0) {
    trips++;
    lastTripAt = now;
  }
}

bool readMpu() {
  // The running average follows with a time constant of 0.6 s, however long
  // the loop took since the last reading (drawing the screen takes a while).
  static unsigned long lastRead = 0;
  unsigned long now = millis();
  float follow = (now - lastRead) / 600.0;
  if (follow > 1) follow = 1;
  lastRead = now;
  Wire.beginTransmission(MPU);
  Wire.write(0x3B);  // ACCEL_XOUT_H: six bytes follow, X, Y, Z, high byte first
  if (Wire.endTransmission(false) != 0) return false;
  if (Wire.requestFrom(MPU, (uint8_t)6, (uint8_t)true) != 6) return false;
  for (int i = 0; i < 3; i++) {
    int16_t raw = (Wire.read() << 8) | Wire.read();
    float g = raw * G_PER_COUNT;
    if (!haveBase) base[i] = g;
    // A slow average: it follows a tilt but not a shake, which comes and goes
    // too fast for it.
    base[i] += (g - base[i]) * follow;
    dyn[i] = g - base[i];
  }
  haveBase = true;
  mag = sqrt(dyn[0] * dyn[0] + dyn[1] * dyn[1] + dyn[2] * dyn[2]);
  return true;
}

// The Modified Mercalli Intensity for a peak ground acceleration in g, after
// the USGS instrumental intensity table.
const char *intensity(float g) {
  if (g < 0.0017) return "I";
  if (g < 0.014) return "II-III";
  if (g < 0.039) return "IV";
  if (g < 0.092) return "V";
  if (g < 0.18) return "VI";
  if (g < 0.34) return "VII";
  if (g < 0.65) return "VIII";
  if (g < 1.24) return "IX";
  return "X+";
}

int pixelsFor(float g) {
  int p = (int)(g * PIXELS_PER_G);
  if (p > SPAN) p = SPAN;
  if (p < -SPAN) p = -SPAN;
  return p;
}

void drawText() {
  display.fillRect(0, 0, 128, 19, SSD1306_BLACK);
  display.setTextColor(SSD1306_WHITE);
  display.setTextSize(1);
  // Line one: the biggest shaking and how many times the switch has tripped.
  display.setCursor(0, 0);
  display.print("PEAK ");
  display.print(peak, 3);
  display.print("g");
  display.setCursor(90, 0);
  display.print("TRIP ");
  display.print(seenTrips);
  // Line two: what is going on, and the Mercalli intensity of the peak.
  display.setCursor(0, 10);
  if (quake) {
    display.print("QUAKE ");
    display.print((millis() - quakeStart) / 100 / 10.0, 1);
    display.print("s");
  } else if (quakeLength > 0) {
    display.print("LAST ");
    display.print(quakeLength / 100 / 10.0, 1);
    display.print("s");
  } else {
    display.print("ARMED");
  }
  display.setCursor(68, 10);
  display.print("MMI ");
  display.print(intensity(peak));
}

void startQuake() {
  quake = true;
  peak = 0;
  quakeStart = millis();
  digitalWrite(LED_PIN, HIGH);
  Serial.println("EARTHQUAKE start");
}

void endQuake() {
  quake = false;
  quakeLength = lastLoud - quakeStart;
  digitalWrite(LED_PIN, LOW);
  Serial.print("EARTHQUAKE end peak ");
  Serial.print(peak, 3);
  Serial.print(" g MMI ");
  Serial.print(intensity(peak));
  Serial.print(" lasted ");
  Serial.print(quakeLength / 100 / 10.0, 1);
  Serial.println(" s");
}

// Draw the column just finished, rub out the three ahead of it, and move on.
void finishColumn() {
  for (int a = 0; a < 2; a++) {
    int centre = a == 0 ? CENTRE_X : CENTRE_Y;
    for (int k = 1; k <= 3; k++) {
      int ahead = (column + k) % 128;
      display.drawFastVLine(ahead, centre - SPAN, 2 * SPAN + 1, SSD1306_BLACK);
      if ((ahead & 3) == 0) display.drawPixel(ahead, centre, SSD1306_WHITE);
    }
    // up on screen is the positive side, so higher values have smaller y
    int top = centre - highest[a];
    int height = highest[a] - lowest[a] + 1;
    display.drawFastVLine(column, centre - SPAN, 2 * SPAN + 1, SSD1306_BLACK);
    if ((column & 3) == 0) display.drawPixel(column, centre, SSD1306_WHITE);
    display.drawFastVLine(column, top, height, SSD1306_WHITE);
    lowest[a] = 127;
    highest[a] = -127;
  }
  column = (column + 1) % 128;
  samplesInColumn = 0;
  // Each display() sends only the 8-column blocks that changed, but it takes
  // time, and the accelerometer is not read while it runs. Showing every
  // fourth column keeps the gaps short.
  if (++columnsSinceShown >= 4) {
    columnsSinceShown = 0;
    drawText();
    display.display();
  }
}

void sample() {
  if (!readMpu()) return;
  unsigned long now = millis();

  // A trip from the vibration switch.
  uint8_t t = trips;
  bool tripped = t != seenTrips;
  seenTrips = t;
  if (tripped) {
    Serial.print("TRIP ");
    Serial.println(t);
  }

  if (mag > START_G || tripped) {
    if (!quake) startQuake();
    lastLoud = now;
  } else if (mag > QUIET_G) {
    lastLoud = now;  // still shaking, just below the start threshold
  }
  if (quake) {
    if (mag > peak) peak = mag;
    if (now - lastLoud > QUIET_MS) endQuake();
  }

  // This sample into the column's highest and lowest, for X and Y.
  for (int a = 0; a < 2; a++) {
    int p = pixelsFor(dyn[a]);
    if (p < lowest[a]) lowest[a] = p;
    if (p > highest[a]) highest[a] = p;
  }
  if (++samplesInColumn >= 5) finishColumn();
}

void setup() {
  Serial.begin(115200);
  Serial.println("Mokxi Uno: seismograph");
  pinMode(LED_PIN, OUTPUT);
  pinMode(TRIP_PIN, INPUT_PULLUP);
  attachInterrupt(digitalPinToInterrupt(TRIP_PIN), onTrip, FALLING);

  Wire.begin();
  Wire.beginTransmission(MPU);  // wake the MPU6050: it starts asleep
  Wire.write(0x6B);             // PWR_MGMT_1
  Wire.write(0);
  Wire.endTransmission();
  Wire.beginTransmission(MPU);  // accelerometer range +-2 g
  Wire.write(0x1C);             // ACCEL_CONFIG
  Wire.write(0);
  Wire.endTransmission();

  display.begin(SSD1306_SWITCHCAPVCC, 0x3C);
  display.clearDisplay();
  for (int a = 0; a < 2; a++) {
    lowest[a] = 127;
    highest[a] = -127;
    int centre = a == 0 ? CENTRE_X : CENTRE_Y;
    for (int x = 0; x < 128; x += 4) display.drawPixel(x, centre, SSD1306_WHITE);
  }
  drawText();
  display.display();
}

void loop() {
  unsigned long now = millis();
  if (now - lastSample >= 10) {  // 100 readings a second
    lastSample = now;
    sample();
  }
  if (now - lastPlot >= 100) {
    lastPlot = now;
    Serial.print("x:");
    Serial.print(dyn[0], 3);
    Serial.print(" y:");
    Serial.print(dyn[1], 3);
    Serial.print(" z:");
    Serial.print(dyn[2], 3);
    Serial.print(" peak:");
    Serial.println(peak, 3);
  }
}