Skip to content

Instantly share code, notes, and snippets.

@kylemcdonald

kylemcdonald/Matlab.pde

Last active Aug 29, 2015
Embed
What would you like to do?
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
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;
}
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));
}
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();
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment