Lomb-Scargle periodogram for Processing, ported from Matlab code, so it's very idiosyncratic and non-idiomatic. Original code published in "Ice Ages and Astronomical Causes: Data, Spectral Analysis and Mechanisms" available at http://books.google.com/books?id=P8ideTkMQisC&pg=PA289&dq=spectral+lomb+scargle&hl=en

_screenshot.png

lomb.pde

float[] lomb(float[] t, float[] y, float[] freq) {
  int nfreq = freq.length;
  float fmax = freq[nfreq - 1];
  float fmin = freq[0];
  float[] power = new float[nfreq];
  float[] f4pi = mult(4*PI, freq);
  float var = cov(y);
  float[] yn = add(-mean(y), y);
  for(int fi = 0; fi < nfreq; fi++) {
    float sinsum = sum(sin(mult(f4pi[fi], t)));
    float cossum = sum(cos(mult(f4pi[fi], t)));
    float tau = atan2(sinsum, cossum);
    float[] argu = mult(TWO_PI * freq[fi], add(-tau, t));
    float[] cosarg = cos(argu);
    float cfi = sum(dotasterisk(yn, cosarg));
    float cosnorm = sum(dotasterisk(cosarg, cosarg));
    float[] sinarg = sin(argu);
    float sfi = sum(dotasterisk(yn, sinarg));
    float sinnorm = sum(dotasterisk(sinarg, sinarg));
    power[fi] = (cfi*cfi/cosnorm+sfi*sfi/sinnorm)/(2*var);
  }
  return power;
}

Matlab.pde

float[] hann(float[] x) {
  int n = x.length;
  float[] result = new float[n];
  for(int i = 0; i < n; i++) {
    float weight = .5 * (1 - cos((TWO_PI * i) / (n - 1)));
    result[i] = x[i] * weight;
  }
  return result;
}

float[] range(float start, float end) {
  return range(start, 1, end);
}

float[] range(float start, float step, float end) {
  float[] result = new float[]{};
  for(float t = start; t < end; t += step) {
    result = append(result, t);
  }
  return result;
}

float[] rand(int count) {
  float[] result = new float[count];
  for(int i = 0; i < count; i++) {
    result[i] = random(1);
  }
  return result;
}

float[] mult(float x, float[] y) {
  int n = y.length;
  float[] result = new float[n];
  for(int i = 0; i < n; i++) {
    result[i] = x * y[i];
  }
  return result;
}

float[] dotasterisk(float[] x, float[] y) {
  int n = y.length;
  float[] result = new float[n];
  for(int i = 0; i < n; i++) {
    result[i] = x[i] * y[i];
  }
  return result;
}

float[] add(float x, float[] y) {
  int n = y.length;
  float[] result = new float[n];
  for(int i = 0; i < n; i++) {
    result[i] = x + y[i];
  }
  return result;
}

float[] add(float[] x, float[] y) {
  int n = y.length;
  float[] result = new float[n];
  for(int i = 0; i < n; i++) {
    result[i] = x[i] + y[i];
  }
  return result;
}

float[] interp1(float[] x, float[] v, float[] xq) {
  int n = xq.length;
  int m = x.length;
  float[] result = new float[n];
  for(int i = 0; i < n; i++) {
    float curq = xq[i];
    int prej = 0, curj = 0;
    float prex = x[0], curx = x[0];
    for(int j = 1; j < m; j++) {
      curj = j;
      curx = x[j];
      if(prex <= curq && curx > curq) {
        break;
      } else {
        prej = curj;
        prex = curx;
      }
    }
    if(prex == curx) {
      result[i] = v[prej];
    } else {
      result[i] = map(curq, prex, curx, v[prej], v[curj]);
    }
  }
  return result;
}

float[] cos(float[] x) {
  int n = x.length;
  float[] result = new float[n];
  for(int i = 0; i < n; i++) {
    result[i] = cos(x[i]);
  }
  return result;
}

float[] sin(float[] x) {
  int n = x.length;
  float[] result = new float[n];
  for(int i = 0; i < n; i++) {
    result[i] = cos(x[i]);
  }
  return result;
}

float sum(float[] x) {
  float sum = 0;
  int n = x.length;
  for(int i = 0; i < n; i++) {
    sum += x[i];
  }
  return sum;
}

float mean(float[] x) {
  return sum(x) / x.length;
}

float[] sq(float[] x) {
  int n = x.length;
  float[] result = new float[n];
  for(int i = 0; i < n; i++) {
    result[i] = sq(x[i]);
  }
  return result;
}

float cov(float[] x) {
  float ssd = sum(sq(add(-mean(x), x)));
  return ssd / x.length;
}

float std(float[] x) {
  return sqrt(cov(x));
}

test_lomb.pde

float[] o18, power;

void setup() {
  size(512, 512);
  prep();
}

void prep() {
  int n = 600;
  float[] age = range(0, n);
//  randomSeed(0);
  float[] ager = add(age, add(-.15, mult(0.3, rand(age.length))));
  ager[0] = age[0]; ager[n - 1] = age[n - 1];
  float[] depth = mult(1/10., age);
  float[] bkg = interp1(range(0, 10, 600), rand(60), ager);
  float f1 = 1./mouseX, f2 = 1./125;
  float[] sig = add(
    cos(mult(TWO_PI*f1, ager)),
    cos(add(PI, mult(TWO_PI*f2, ager))));
  o18 = add(sig, bkg);
//  o18 = hann(o18);
  // pick frequencies to evaluate spectral power
  float[] freq = range(0, 0.0001, 0.02);
  power = lomb(age, o18, freq);
  power[0] = 0;
  // normalize to average 1
  power = mult(1./std(power), power);
}

void draw() {
  prep();
  background(255);
  noFill();
  stroke(0);
  
  // plot the results
  beginShape();
  vertex(0, height);
  for(int i = 0; i < power.length; i++) {
    float x = map(i, 0, power.length, 0, width);
    float y = map(power[i], 0, 6, height, 0);
    vertex(x, y);
  }
  vertex(width, height);
  endShape();
  
  beginShape();
  for(int i = 0; i < o18.length; i++) {
    float x = map(i, 0, o18.length, 0, width);
    float y = map(o18[i], -10, 10, 0, height);
    vertex(x, y);
  }
  endShape();
}