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.
forked from efekarakus‘s block: Find Peaks
<!DOCTYPE html>
<meta charset="utf-8">
<title>Find Peaks</title>
<style>
.axis path,
.axis line {
fill: none;
stroke: black;
shape-rendering: crispEdges;
}
text {
font-family: sans-serif;
font-size: 11px;
}
path {
fill: none;
stroke: black;
stroke-width: 1.5px;
}
path.cardinal {
stroke: crimson;
stroke-width: 3px;
}
path.area {
stroke: black;
stroke-width: 0.5px;
fill: crimson;
fill-opacity: 0.2;
}
circle.peak {
fill: crimson;
}
</style>
<body>
<script src="//d3js.org/d3.v3.min.js" charset="utf-8"></script>
<script src="d3-peaks.js" charset="utf-8"></script>
<script>
var width = 920,
height = 464,
margin = {top: 10, left: 10, bottom: 20, right: 10};
var ricker = d3_peaks.ricker;
var findPeaks = d3_peaks.findPeaks()
.kernel(ricker)
.gapThreshold(1)
.minLineLength(2)
.minSNR(1.0)
.widths([1,2,3]);
var signal = d3.range(0,200).map(function(d) {return Math.random()*20;});
var svg = d3.select("body").append("svg")
.attr("width", width + margin.left + margin.right)
.attr("height", height + margin.top + margin.bottom);
var g = svg.append("g")
.attr("transform", "translate(" + margin.left + "," + margin.top + ")");
var x = d3.scale.linear()
.domain ([0, signal.length])
.range([20, width]);
var y = d3.scale.linear()
.domain([0, d3.max(signal)])
.range([height, 0]);
var xAxis = d3.svg.axis()
.scale(x)
.ticks(20);
var line = d3.svg.line()
.x(function(d, i) { return x(i); })
.y(function(d) { return y(d); });
var area = d3.svg.area()
.x(function(d) { return x(d.x); })
.y0(height)
.y1(function(d) { return y(d.y); });
g.append("g")
.datum(signal)
.append("path")
.attr("d", line);
g.append("g")
.attr("class", "axis")
.attr("transform", "translate(0," + y(0) + ")")
.call(xAxis);
function draw(peaks,i){
g.selectAll(".peak"+i)
.data(peaks)
.enter().append("circle")
.attr("cx", function(peak) { return x(peak.index); })
.attr("cy", function(peak) { return y(signal[peak.index]); })
.attr("r", 4)
.attr("class", "peak"+i);
line.interpolate("cardinal")
.x(function(d) { return x(d.index); })
.y(function(d) { return y(signal[d.index]); });
peaks.sort(function(a, b) { return a.index - b.index; });
g.append("g")
.datum(peaks)
.append("path")
.attr("d", line)
.attr("class", "cardinal");
}
draw(findPeaks(signal),1)
draw(findPeaks(signal.map(function(d){return -d})),2)
var areas = [];
findPeaks(signal).forEach(function(peak) {
var width = peak.width;
var lowerBound = Math.max(0, peak.index - width);
var upperBound = Math.min(signal.length, peak.index + width + 1);
var a = [];
for (var i = lowerBound; i < upperBound; i++) {
a.push({
x: i,
y: signal[i]
});
}
areas.push(a);
});
areas.forEach(function(data) {
g.append("g")
.datum(data)
.append("path")
.attr("d", area)
.attr("class", "area");
})
</script>
</body>(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;
}));{"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, 3, 6, 6, 4, 4, 8, 11, 8, 9, 3, 3, 3, 5, 6, 4, 7, 9, 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]}