Instantly share code, notes, and snippets.

Find Peaks

This example illustrates how to use the findPeaks API. You can download the library from https://github.com/efekarakus/d3-peaks.

The algorithm is based on "Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching" by Pan Du, Warren A. Kibbe and Simon M. Lin. The paper can be found here.

 (function (global, factory) { typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports) : typeof define === 'function' && define.amd ? define(['exports'], factory) : (factory((global.d3_peaks = global.d3_peaks || {}))); }(this, function (exports) { 'use strict'; /** * See https://en.wikipedia.org/wiki/Mexican_hat_wavelet */ function ricker() { var σ = 1; var ricker = function(t) { var t2 = t*t, variance = σ*σ; var C = 2.0 / ( Math.sqrt(3 * σ) * (Math.pow(Math.PI, 0.25)) ); var norm = (1.0 - (t2)/(variance)); var gauss = Math.exp( -(t2) / (2*variance) ); return C*norm*gauss; } ricker.std = function(_) { return arguments.length ? (σ = _, ricker) : σ; } /** * Range of points to sample from the wavelet. [-reach, reach] */ ricker.reach = function() { return 5 * σ; } return ricker; }; function convolve() { var kernel = ricker(); /** * y[n] = Sum_k{x[k] * h[n-k]} * y: output * x: input * h: smoother */ var convolve = function(signal) { var size = signal.length, n = -1, convolution = new Array(size); while (++n < size) { var y = 0; var box = boundingBox(n, kernel.reach(), 0, size - 1); box.forEach(function(δ) { var k = n + δ; y += signal[k] * kernel(δ); }); convolution[n] = y; } return convolution; }; convolve.kernel = function(_) { return arguments.length ? (kernel = _, convolve) : kernel; } function range(reach) { reach = +reach; var i = -1, n = 2*reach + 1, range = new Array(n); while(++i < n) { range[i] = (-reach) + i; } return range; } function boundingBox(n, reach, lo, hi) { for (var i = 1; i <= reach; i++) { var left = n - i, right = n + i; if (left >= lo && right <= hi) continue; return range(i - 1); } return range(reach); } return convolve; }; function isLocalMaxima(arr, index) { var current = arr[index], left = arr[index - 1], right = arr[index + 1]; if (left !== undefined && right !== undefined) { if (current > left && current > right) { return true; } else if (current >= left && current > right) { return true; } else if (current > left && current >= right) { return true; } } else if (left !== undefined && current > left) { return true; } else if (right !== undefined && current > right) { return true; } return false; } /** * @param {arr} row in the CWT matrix. * @return Array of indices with relative maximas. */ function maximas(arr) { var maximas = []; var length = arr.length; arr.forEach(function(value, index) { if (isLocalMaxima(arr, index)) maximas.push({x: index, y: value}); }); return maximas; }; function nearestNeighbor(line, maximas, window) { var cache = {}; maximas.forEach(function(d) { cache[d.x] = d.y; }); var point = line.top(); for (var i = 0; i <= window; i++) { var left = point.x + i; var right = point.x - i; if ( (left in cache) && (right in cache) ) { if (cache[left] > cache[right]) { return left; } return right; } else if (left in cache) { return left; } else if (right in cache) { return right; } } return null; } function percentile(arr, perc) { var length = arr.length; var index = Math.min(length - 1, Math.ceil(perc * length)); arr.sort(function(a, b) { return a - b; }); return arr[index]; } function Point(x, y, width) { this.x = x; this.y = y; this.width = width; this.snr = undefined; } Point.prototype.SNR = function(conv) { var smoothingFactor = 0.00001; var signal = this.y; var lowerBound = Math.max(0, this.x - this.width); var upperBound = Math.min(conv.length, this.x + this.width + 1); var neighbors = conv.slice(lowerBound, upperBound); var noise = percentile(neighbors, 0.95); signal += smoothingFactor; noise += smoothingFactor; this.snr = signal/noise; return this.snr; } Point.prototype.serialize = function() { return {index: this.x, width: this.width, snr: this.snr}; } function RidgeLine() { this.points = []; this.gap = 0; } /** * If the point is valid append it to the ridgeline, and reset the gap. * Otherwise, increment the gap and do nothing. * * @param {point} Point object. */ RidgeLine.prototype.add = function(point) { if (point === null || point === undefined) { this.gap += 1; return; } else { this.points.push(point); this.gap = 0; } } /** * @return {Point} Last point added into the ridgeline. */ RidgeLine.prototype.top = function() { return this.points[this.points.length - 1]; } /** * @return {number} Length of points on the ridgeline. */ RidgeLine.prototype.length = function() { return this.points.length; } /** * @return {boolean} True if the gap in the line is above a threshold. False otherwise. */ RidgeLine.prototype.isDisconnected = function (threshold) { return this.gap > threshold; } /** * @param {Array} Smallest scale in the convolution matrix */ RidgeLine.prototype.SNR = function(conv) { var maxSnr = Number.NEGATIVE_INFINITY; this.points.forEach(function(point) { var snr = point.SNR(conv); if (snr > maxSnr) maxSnr = snr; }); return maxSnr; } function findPeaks() { var kernel = ricker, gapThreshold = 1, minLineLength = 1, minSNR = 1.0, widths = [1]; var findPeaks = function(signal) { var M = CWT(signal); var ridgeLines = initializeRidgeLines(M); ridgeLines = connectRidgeLines(M, ridgeLines); ridgeLines = filterRidgeLines(signal, ridgeLines); return peaks(signal, ridgeLines); }; /** * Smoothing function. */ findPeaks.kernel = function(_) { return arguments.length ? (kernel = _, findPeaks) : kernel; } /** * Expected widths of the peaks. */ findPeaks.widths = function(_) { _.sort(function(a, b) { return a - b; }); return arguments.length ? (widths = _, findPeaks) : widths; } /** * Number of gaps that we allow in the ridge lines. */ findPeaks.gapThreshold = function(_) { return arguments.length ? (gapThreshold = _, findPeaks) : gapThreshold; } /** * Minimum ridge line length. */ findPeaks.minLineLength = function(_) { return arguments.length ? (minLineLength = _, findPeaks) : minLineLength; } /** * Minimum signal to noise ratio for the peaks. */ findPeaks.minSNR = function(_) { return arguments.length ? (minSNR = _, findPeaks) : minSNR; } var CWT = function(signal) { var M = new Array(widths.length); widths.forEach(function(width, i) { var smoother = kernel() .std(width); var transform = convolve() .kernel(smoother); var convolution = transform(signal); M[i] = convolution; }); return M; } var initializeRidgeLines = function(M) { var n = widths.length; var locals = maximas(M[n - 1], widths[n - 1]); var ridgeLines = []; locals.forEach(function(d) { var point = new Point(d.x, d.y, widths[n - 1]); var line = new RidgeLine(); line.add(point); ridgeLines.push(line); }); return ridgeLines; } var connectRidgeLines = function(M, ridgeLines) { var n = widths.length; for (var row = n - 2; row >= 0; row--) { var locals = maximas(M[row], widths[row]); var addedLocals = []; // Find nearest neighbor at next scale and add to the line ridgeLines.forEach(function(line, i) { var x = nearestNeighbor(line, locals, widths[row]); line.add(x === null ? null : new Point(x, M[row][x], widths[row])); if (x !== null) { addedLocals.push(x); } }); // Remove lines that has exceeded the gap threshold ridgeLines = ridgeLines.filter(function(line) { return !line.isDisconnected(gapThreshold); }); // Add all the unitialized ridge lines locals.forEach(function(d) { if (addedLocals.indexOf(d.x) !== -1) return; var point = new Point(d.x, d.y, widths[row]); var ridgeLine = new RidgeLine(); ridgeLine.add(point); ridgeLines.push(ridgeLine); }); } return ridgeLines; } var filterRidgeLines = function(signal, ridgeLines) { var smoother = kernel() .std(1.0); var transform = convolve() .kernel(smoother); var convolution = transform(signal); ridgeLines = ridgeLines.filter(function(line) { var snr = line.SNR(convolution); return (snr >= minSNR) && (line.length() >= minLineLength); }); return ridgeLines } /** * Pick the point on the ridgeline with highest SNR. */ var peaks = function(signal, ridgeLines) { var peaks = ridgeLines.map(function(line) { var points = line.points; var maxValue = Number.NEGATIVE_INFINITY, maxPoint = undefined; points.forEach(function(point) { var y = signal[point.x]; if (y > maxValue) { maxPoint = point; maxValue = y; } }); return maxPoint.serialize(); }); return peaks; } return findPeaks; }; var version = "0.0.1"; exports.version = version; exports.ricker = ricker; exports.convolve = convolve; exports.findPeaks = findPeaks; }));
 Find Peaks
 {"signal": [ 1, 10, 12, 12, 8, 15, 10, 4, 13, 11, 8, 7, 8, 9, 3, 10, 8, 5, 9, 6, 9, 7, 17, 5, 8, 6, 4, 5, 3, 5, 9, 6, 3, 8, 1, 7, 4, 8, 7, 5, 3, 6, 3, 13, 6, 6, 9, 6, 8, 9, 5, 6, 6, 5, 9, 6, 5, 7, 7, 7, 7, 6, 6, 7, 6, 4, 7, 8, 6, 18, 6, 11, 12, 7, 19, 13, 2, 14, 16, 4, 11, 10, 11, 10, 8, 6, 5, 6, 5, 6, 6, 6, 4, 5, 3, 7, 6, 7, 5, 5, 6, 3, 7, 2, 3, 1, 3, 4, 3, 2, 4, 4, 1, 2, 1, 6, 1, 1, 1, 4, 3, 1, 5, 1, 2, 3, 2, 6, 2, 2, 3, 2, 2, 2, 1, 4, 7, 6, 5, 7, 7, 6, 4, 12, 9, 7, 5, 2, 4, 7, 5, 4, 4, 4, 2, 5, 3, 3, 3, 5, 7, 4, 4, 4, 2, 4, 1, 3, 2, 3, 4, 13,13, 13, 13, 12,11, 11, 8, 9, 3, 3, 3,6, 6,6, 6, 6, 7, 13, 9, 7, 12, 7, 14, 11, 8, 5, 2, 6, 6, 5, 5, 2, 4, 2, 4, 6, 6, 14, 16, 16, 18, 14, 3, 18, 26, 2, 23, 26, 8, 18, 16, 3, 14, 27, 4, 10, 14, 20, 2, 15, 22, 8, 9, 11, 15, 8, 8, 15]}