diff --git a/docs/src/assets/dist_peaks.png b/docs/src/assets/dist_peaks.png new file mode 100644 index 000000000..ac7be1098 Binary files /dev/null and b/docs/src/assets/dist_peaks.png differ diff --git a/docs/src/assets/prom_peaks.png b/docs/src/assets/prom_peaks.png new file mode 100644 index 000000000..377eacdbb Binary files /dev/null and b/docs/src/assets/prom_peaks.png differ diff --git a/docs/src/contents.md b/docs/src/contents.md index 3b6d3c818..b4b496560 100644 --- a/docs/src/contents.md +++ b/docs/src/contents.md @@ -13,5 +13,6 @@ Pages = ["periodograms.md", "convolutions.md", "lpc.md", "index.md", + "findpeaks.md", ] ``` diff --git a/docs/src/findpeaks.md b/docs/src/findpeaks.md new file mode 100644 index 000000000..b4dc52d33 --- /dev/null +++ b/docs/src/findpeaks.md @@ -0,0 +1,51 @@ +# `Findpeaks` + +`findpeaks` function is inspired by MATLAB's findpeaks function. + +The function allows you to get positions of local maxima of one-dimensional data. +Since real-world data are often noisy, there is usually a large number of local maxima that are of no interest. + +`findpeaks` offers filtering based on 4 parameters (can be combined): +* peak height +* distance between peaks +* peak [prominence](https://en.wikipedia.org/wiki/Topographic_prominence) +* threshold - minimal difference between neighboring data points + +```@docs +findpeaks +``` + +# Example + +```julia +using Findpeaks + +using DelimitedFiles + +using Plots +default(legend=false) + + +data = readdlm("test/data/findpeaks_example_spectrum.txt") +x = data[:, 1] +y = data[:, 2] + +peaks = findpeaks(y, x, min_prom=1000.) + +plot(x, y, title="Prominent peaks") +scatter!(x[peaks], y[peaks]) +``` + +![](assets/prom_peaks.png) + +```julia +sep = 0.2 + +peaks = findpeaks(y, x, min_dist=sep) + +plot(x, y, title="Peaks at least $sep units apart") +scatter!(x[peaks], y[peaks]) +``` + +![](assets/dist_peaks.png) + diff --git a/src/DSP.jl b/src/DSP.jl index da0eb5075..bdb2a9377 100644 --- a/src/DSP.jl +++ b/src/DSP.jl @@ -15,9 +15,10 @@ include("periodograms.jl") include("Filters/Filters.jl") include("lpc.jl") include("estimation.jl") +include("findpeaks.jl") using Reexport -@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Unwrap, .Estimation +@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Unwrap, .Estimation, .Findpeaks include("deprecated.jl") end diff --git a/src/findpeaks.jl b/src/findpeaks.jl new file mode 100644 index 000000000..eb110dc72 --- /dev/null +++ b/src/findpeaks.jl @@ -0,0 +1,320 @@ +module Findpeaks + +export findpeaks, PeakInfo + +struct PeakInfo{T} + height :: T + left_diff :: T + right_diff :: T + prominence :: T + plateau_range :: UnitRange{Int} +end + +""" +Overwrites arrays specified in `args` by `args[inds]`. +""" +macro apply_indices!(inds, args...) + Expr( + :block, + [ :($(esc(args[i])) = $(esc(args[i]))[$(esc(inds))]) for i in 1:length(args) ]...) +end + +""" +Macro specific to findpeaks function." +""" +macro on_empty_return(peaks, x) + quote + if isempty($(esc(peaks))) + return (empty($(esc(x))), empty($(esc(x)), PeakInfo)) + end + end +end + +""" +`findpeaks(y::Array{T}, +x::Array{S}=collect(1:length(y)) +;min_height::T=minimum(y), min_prom::T=minimum(y), +min_dist::S=0, threshold::T=0 ) where {T<:Real,S}`\n + +Finds peaks in 1D array. Similar to MATLAB's findpeaks().\n +Returns:\n +* a vector of peak positions found +* a vector of `PeakInfo` structures with additional information about peaks. + +*Arguments*:\n +`y` -- data\n +*Optional*:\n +`x` -- x-data\n +*Keyword*:\n +`min_height` -- minimal peak height\n +`min_prom` -- minimal peak prominence\n +`min_dist` -- minimal peak distance (keeping highest peaks)\n +`threshold` -- minimal difference (absolute value) between + peak and neighboring points\n +`min_plateau_points` -- minimal number of points a plateau needs to have to be considered a peak +`max_plateau_points` -- maximal number of points a plateau can have to be considered a peak + +The order of filtering is: +**height -> threshold -> prominence -> distance -> plateaus** + +*Warning* `NaN` values are currently not supported and any `NaN`s should be filtered out by the user before applying this function. +""" +function findpeaks( + y :: AbstractVector{T}, + x :: AbstractVector{S} = collect(1:length(y)) + ; + kwargs... + ) where {T, S} + + @on_empty_return(y, x) + + if length(x) != length(y) + lx, ly = length(x), length(y) + throw(ArgumentError("`x` and `y` need to have the same length: x($lx) y($ly)")) + end + + + new_peaks, new_peak_info + min_height = get(kwargs, :min_height, minimum(y)) + min_prom = get(kwargs, :min_prom, zero(y[1])) + min_dist = get(kwargs, :min_dist, zero(x[1])) + threshold = get(kwargs, :threshold, zero(y[1])) + min_plateau_points = get(kwargs, :min_plateau_points, zero(length(y))) + max_plateau_points = get(kwargs, :max_plateau_points, typemax(Int)) + + peaks = 1:length(y) |> collect + + heights = y[peaks] + inds2keep = heights .> min_height + @apply_indices!(inds2keep, peaks, heights) + @on_empty_return(peaks, x) + + inds2keep, ldiffs, rdiffs = in_threshold(peaks, y, threshold) + @apply_indices!(inds2keep, peaks, ldiffs, rdiffs, heights) + @on_empty_return(peaks, x) + + inds2keep, proms = with_prominence(peaks, y, min_prom) + @apply_indices!(inds2keep, peaks, ldiffs, rdiffs, heights, proms) + @on_empty_return(peaks, x) + + inds2keep = with_distance(peaks, x, y, min_dist) + @apply_indices!(inds2keep, peaks, ldiffs, rdiffs, proms, heights) + @on_empty_return(peaks, x) + + peak_info = [PeakInfo{T}(heights[i], ldiffs[i], rdiffs[i], proms[i], 1:1) for i in 1:length(inds2keep)] + + peaks, peak_info = group_plateaus(peaks, peak_info, min_plateau_points, max_plateau_points) + (peaks = peaks, peak_info = peak_info) +end + +function in_threshold( + peaks :: AbstractVector{Int}, + y :: AbstractVector{T}, + threshold :: T, + ) where {T <: Real} + + peak_len = length(peaks) + + inds2keep = 1:peak_len |> collect + rdiffs = zeros(T, peak_len) + ldiffs = zeros(T, peak_len) + + dy = diff(y) + + n_peaks2keep = 0 + for (j, peak_i) in enumerate(peaks) + # resolve edges of the data + if peak_i == 1 + # only check right difference + rd = -dy[peak_i] + if rd >= threshold + n_peaks2keep += 1 + inds2keep[n_peaks2keep] = j + ldiffs[j] = zero(threshold) + rdiffs[j] = rd + end + elseif peak_i == length(y) + # only check left difference + ld = dy[peak_i - 1] + if ld >= threshold + n_peaks2keep += 1 + inds2keep[n_peaks2keep] = j + ldiffs[j] = ld + rdiffs[j] = zero(threshold) + end + else + # check both sides + ld = dy[peak_i - 1] + rd = -dy[peak_i] + if rd >= threshold && ld >= threshold + n_peaks2keep += 1 + inds2keep[n_peaks2keep] = j + ldiffs[j] = ld + rdiffs[j] = rd + end + end + end + + (inds2keep[1:n_peaks2keep], ldiffs, rdiffs) +end + +function with_prominence( + peaks :: AbstractVector{Int}, + y :: AbstractVector{T}, + min_prom::T, + ) where {T <: Real} + + #minimal prominence refinement + proms = prominence(y, peaks) + proms .> min_prom, proms +end + + +function prominence(y :: AbstractVector{T}, + peaks :: AbstractVector{Int} + ) where {T <: Real} + yP = y[peaks] + proms = zero(yP) + + for (i, p) in enumerate(peaks) + lP, rP = 1, length(y) + for j = (i-1):-1:1 + if yP[j] > yP[i] + lP = peaks[j] + break + end + end + ml = minimum(y[lP:p]) + for j = (i+1):length(yP) + if yP[j] > yP[i] + rP = peaks[j] + break + end + end + mr = minimum(y[p:rP]) + ref = max(mr,ml) + proms[i] = yP[i] - ref + end + + proms +end + +""" +Select only peaks that are further apart than `min_dist` +""" +function with_distance( + peaks :: AbstractVector{Int}, + x :: AbstractVector{S}, + y :: AbstractVector{T}, + min_dist::S, + ) where {T <: Real, S} + + peaks2del = zeros(Bool, length(peaks)) + inds = sortperm(y[peaks], rev=true) + sorted_peaks = copy(peaks) + permute!(sorted_peaks, inds) + for i = 1:length(sorted_peaks) + for j = 1:(i-1) + if abs(x[sorted_peaks[i]] - x[sorted_peaks[j]]) <= min_dist + if !peaks2del[j] + peaks2del[i] = true + end + end + end + end + + inds[.!peaks2del] +end + +function find_plateaus(peaks :: AbstractVector{Int}) + n_ranges = 0 + ranges = [1:1 for i in 1:length(peaks)] + is_plateau = diff(peaks) .== 1 + constructed_range = -1:-1 + for (i, t) in enumerate(is_plateau) + if t + if constructed_range == -1:-1 + # start plateau range + constructed_range = i:i + end + else + if constructed_range != -1:-1 + # finish plateau range + constructed_range = constructed_range.start : i + n_ranges += 1 + ranges[n_ranges] = constructed_range + + constructed_range = -1:-1 + else + # add single peak + n_ranges += 1 + ranges[n_ranges] = i:i + end + end + end + # finish the last plateau + if constructed_range != -1:-1 + constructed_range = constructed_range.start : length(peaks) + n_ranges += 1 + ranges[n_ranges] = constructed_range + else + n_ranges += 1 + ranges[n_ranges] = length(peaks) : length(peaks) + end + ranges[1:n_ranges] +end + + +function merge_l_r_diff(l_edge :: PeakInfo, r_edge :: PeakInfo) + (l_edge.left_diff, r_edge.right_diff) +end + +function global_plateau_range(peaks :: AbstractVector{Int}, + peak_info :: Vector{PeakInfo{T}}, + r :: UnitRange{Int}) where T <: Real + l_edge = peaks[r.start] + r_edge = peaks[r.stop] + plateau_range = l_edge : r_edge + + l_edge_info = peak_info[r.start] + r_edge_info = peak_info[r.stop] + left_diff, right_diff = merge_l_r_diff(l_edge_info, r_edge_info) + merged_peak_info = PeakInfo( + l_edge_info.height, + left_diff, + right_diff, + l_edge_info.prominence, + plateau_range, + ) + + (plateau_range, merged_peak_info) +end + +function center_peak(plateau :: UnitRange{Int}) + center_peak = floor(Int, (plateau.start + plateau.stop)/2) +end + +function group_plateaus(peaks :: AbstractVector{Int}, + peak_info :: Vector{PeakInfo{T}}, + min_plateau_points :: Int, + max_plateau_points :: Int) where T <: Real + + inds = sortperm(peaks) + @apply_indices!(inds, peaks, peak_info) + + plateaus = find_plateaus(peaks) + ranges = filter(x -> min_plateau_points <= length(x) <= max_plateau_points, plateaus) + + l = length(ranges) + new_peaks = Vector{Int}(undef, l) + new_peak_info = Vector{PeakInfo}(undef, l) + for (i, r) in enumerate(ranges) + plateau_range, merged_peak_info = global_plateau_range(peaks, peak_info, r) + new_peaks[i] = center_peak(plateau_range) + new_peak_info[i] = merged_peak_info + end + + new_peaks, new_peak_info +end + +end # module diff --git a/test/data/findpeaks_example_spectrum.txt b/test/data/findpeaks_example_spectrum.txt new file mode 100644 index 000000000..04249eeb4 --- /dev/null +++ b/test/data/findpeaks_example_spectrum.txt @@ -0,0 +1,300 @@ +301.12555 372.666656494 +301.13156 379.666656494 +301.13754 360.333312988 +301.14352 382.0 +301.14954 369.0 +301.15552 368.0 +301.16153 373.333312988 +301.16751 369.333312988 +301.17349 382.0 +301.1795 377.0 +301.18549 357.333312988 +301.19147 348.333312988 +301.19748 355.333312988 +301.20346 374.333312988 +301.20947 369.666656494 +301.21545 364.333312988 +301.22144 351.0 +301.22745 339.0 +301.23343 352.666656494 +301.23941 354.666656494 +301.24542 360.0 +301.2514 368.0 +301.25739 344.666656494 +301.2634 394.333312988 +301.26938 408.333312988 +301.27536 399.333312988 +301.28134 421.666656494 +301.28735 490.333312988 +301.29333 619.333312988 +301.29932 607.666656494 +301.30533 498.333312988 +301.31131 386.666656494 +301.31729 388.333312988 +301.32327 385.666656494 +301.32928 374.666656494 +301.33527 330.666656494 +301.34125 368.333312988 +301.34726 387.333312988 +301.35324 375.333312988 +301.35922 377.0 +301.3652 354.0 +301.37122 379.0 +301.3772 368.0 +301.38318 338.666656494 +301.38916 363.0 +301.39517 348.666656494 +301.40115 382.0 +301.40714 358.333312988 +301.41312 365.0 +301.4191 396.666656494 +301.42511 365.666656494 +301.43109 369.666656494 +301.43707 369.0 +301.44305 367.666656494 +301.44907 384.666656494 +301.45505 350.666656494 +301.46103 359.0 +301.46701 381.0 +301.47299 378.333312988 +301.47897 357.0 +301.48499 357.333312988 +301.49097 369.333312988 +301.49695 353.666656494 +301.50293 351.666656494 +301.50891 380.333312988 +301.51489 366.666656494 +301.5209 338.333312988 +301.52689 371.666656494 +301.53287 376.333312988 +301.53885 362.666656494 +301.54483 338.666656494 +301.55081 371.0 +301.55682 345.666656494 +301.56281 389.0 +301.56879 396.0 +301.57477 352.0 +301.58075 373.333312988 +301.58673 382.666656494 +301.59271 367.0 +301.59869 343.333312988 +301.60468 370.666656494 +301.61069 361.666656494 +301.61667 372.666656494 +301.62265 386.666656494 +301.62863 338.666656494 +301.63461 367.0 +301.64059 371.666656494 +301.64658 363.333312988 +301.65256 362.0 +301.65854 357.333312988 +301.66452 379.666656494 +301.6705 369.333312988 +301.67648 394.0 +301.6825 382.333312988 +301.68848 370.666656494 +301.69446 360.0 +301.70044 377.333312988 +301.70642 391.333312988 +301.7124 382.333312988 +301.71838 352.333312988 +301.72437 354.0 +301.73035 363.0 +301.73633 407.333312988 +301.74231 410.333312988 +301.74829 413.0 +301.75427 447.666656494 +301.76025 525.0 +301.76624 757.333312988 +301.77222 913.333312988 +301.7782 679.333312988 +301.78418 474.333312988 +301.79016 399.0 +301.79614 369.666656494 +301.80212 375.666656494 +301.80811 352.666656494 +301.81409 348.333312988 +301.82007 375.333312988 +301.82605 383.666656494 +301.83203 390.333312988 +301.83801 371.666656494 +301.84399 360.666656494 +301.84998 375.333312988 +301.85596 372.0 +301.86194 363.0 +301.86789 374.666656494 +301.87387 350.0 +301.87985 362.666656494 +301.88583 365.333312988 +301.89182 370.333312988 +301.8978 373.666656494 +301.90378 389.0 +301.90976 394.666656494 +301.91574 375.0 +301.92172 373.666656494 +301.9277 348.0 +301.93369 375.333312988 +301.93964 376.666656494 +301.94562 375.666656494 +301.9516 364.333312988 +301.95758 384.333312988 +301.96356 369.0 +301.96954 366.666656494 +301.97552 370.666656494 +301.98151 365.666656494 +301.98746 367.333312988 +301.99344 380.0 +301.99942 372.333312988 +302.0054 387.333312988 +302.01138 381.0 +302.01736 367.666656494 +302.02335 373.666656494 +302.0293 386.0 +302.03528 369.333312988 +302.04126 364.666656494 +302.04724 378.0 +302.05322 379.666656494 +302.0592 391.333312988 +302.06516 377.0 +302.07114 373.666656494 +302.07712 362.0 +302.0831 395.666656494 +302.08908 393.333312988 +302.09503 392.333312988 +302.10101 374.0 +302.10699 379.333312988 +302.11298 352.0 +302.11896 402.333312988 +302.12491 403.0 +302.13089 383.0 +302.13687 381.666656494 +302.14285 389.0 +302.1488 400.0 +302.15479 373.0 +302.16077 375.333312988 +302.16675 368.666656494 +302.17273 392.666656494 +302.17868 378.666656494 +302.18466 369.333312988 +302.19064 370.333312988 +302.19659 363.666656494 +302.20258 387.0 +302.20856 402.333312988 +302.21454 361.333312988 +302.22049 397.666656494 +302.22647 407.333312988 +302.23245 435.333312988 +302.23843 452.666656494 +302.24438 469.0 +302.25037 607.0 +302.25635 918.0 +302.2623 930.333312988 +302.26828 658.666656494 +302.27426 471.0 +302.28021 389.0 +302.28619 391.666656494 +302.29218 375.0 +302.29813 381.333312988 +302.30411 392.666656494 +302.31009 376.666656494 +302.31604 362.666656494 +302.32202 391.666656494 +302.328 380.333312988 +302.33395 422.0 +302.33994 381.333312988 +302.34592 420.333312988 +302.35187 488.666656494 +302.35785 577.333312988 +302.36383 710.333312988 +302.36978 721.666656494 +302.37576 814.666656494 +302.38171 1005.333312988 +302.3877 844.0 +302.39368 594.333312988 +302.39963 511.0 +302.40561 513.0 +302.41156 529.0 +302.41754 551.0 +302.42352 684.666656494 +302.42947 966.666656494 +302.43546 955.666656494 +302.44141 709.333312988 +302.44739 529.666656494 +302.45337 459.0 +302.45932 441.333312988 +302.4653 417.333312988 +302.47125 422.0 +302.47723 440.666656494 +302.48318 471.0 +302.48917 502.0 +302.49512 560.333312988 +302.5011 874.0 +302.50708 1121.666656494 +302.51303 958.333312988 +302.51901 942.666656494 +302.52496 839.666656494 +302.53094 600.0 +302.5369 470.0 +302.54288 440.333312988 +302.54883 419.0 +302.55481 411.666656494 +302.56076 420.0 +302.56674 409.0 +302.57269 413.0 +302.57867 504.0 +302.58463 511.0 +302.59061 585.666656494 +302.59656 694.666656494 +302.60254 1023.333312988 +302.60849 1824.0 +302.61447 1838.0 +302.62042 1314.333312988 +302.6264 1046.333312988 +302.63235 691.666656494 +302.63831 641.666656494 +302.64429 655.0 +302.65024 640.333312988 +302.65622 488.333312988 +302.66217 419.666656494 +302.66815 395.333312988 +302.6741 404.0 +302.68008 396.333312988 +302.68604 420.666656494 +302.69199 437.0 +302.69797 442.0 +302.70392 523.333312988 +302.7099 615.0 +302.71585 737.333312988 +302.7218 867.666656494 +302.72778 1214.666656494 +302.73373 2404.0 +302.73972 3620.333312988 +302.74567 2985.666656494 +302.75162 2194.666656494 +302.7576 1377.0 +302.76355 795.0 +302.76953 599.666656494 +302.77548 670.333312988 +302.78143 835.666656494 +302.78741 861.0 +302.79337 866.0 +302.79932 792.333312988 +302.8053 895.666656494 +302.81125 959.333312988 +302.8172 747.0 +302.82318 563.0 +302.82913 520.0 +302.83508 504.333312988 +302.84106 605.0 +302.84702 696.666656494 +302.85297 761.666656494 +302.85895 862.0 +302.8649 1013.666656494 +302.87085 1316.0 +302.87683 2119.0 +302.88278 4045.666656494 +302.88873 6054.666656494 +302.89468 5485.0 +302.90067 3806.333312988 +302.90662 1878.666656494 +302.91257 918.0 diff --git a/test/findpeaks.jl b/test/findpeaks.jl new file mode 100644 index 000000000..dcaad44f5 --- /dev/null +++ b/test/findpeaks.jl @@ -0,0 +1,153 @@ +using DSP +using Test + +const NAME = "Findpeaks: " + +@testset "$NAME Single Gaussian peak" begin + gaussian(x,μ,σ) = 1/sqrt(2*π)/σ*exp(-((x-μ)^2)/2/σ^2) + x = collect(range(1,stop=1000,length=1000)) + μ = rand(x) + data = gaussian.(x,μ,rand(x)/3) + + p, _ = findpeaks(data) + @test length(p) == 1 + @test p[1] == μ + + p, _ = findpeaks(data,x) + @test length(p) == 1 + @test p[1] == μ +end + + +@testset "$NAME Prominence" begin + gaussian(x,μ,σ) = exp(-((x-μ)^2)/2/σ^2) + + x = range(0., stop = 1., length = 1000) |> collect + + positions = [1, 2, 3, 5, 7, 8]/10; + amps = [3, 7, 5, 5, 4, 5]; + widths = [1, 3, 3, 4, 2, 3]/100; + + base = 4 * cos.(2π * x); + signal = ( a * gaussian.(x, μ, σ) for (a, μ, σ) in zip(amps, positions, widths) ) + data = sum(signal) + base + + peaks, _ = findpeaks(data, x, min_prom=4.) + + @test length(peaks) == 2 + + expected_peak_1 = argmin(abs.(x .- positions[2])) # global highest peak + expected_peak_2 = argmin(abs.(x .- positions[4])) # lowest peak + + # check peaks are around expected -> they may be shifted because of background + @test abs(peaks[1] - expected_peak_1) < 20 + @test abs(peaks[2] - expected_peak_2) < 20 +end + +@testset "$NAME Threshold" begin + y = [0., 1., 3., 1., 4., 5., 3., 0., 2., 5., 4., 0.] + + @test Set(findpeaks(y, threshold=3.1)[1]) == Set([]) + @test Set(findpeaks(y, threshold=2.1)[1]) == Set([]) + @test Set(findpeaks(y, threshold=1.1)[1]) == Set([3]) + @test Set(findpeaks(y, threshold=0.1)[1]) == Set([3, 6, 10]) + @test Set(findpeaks(y, threshold=0.0)[1]) == Set([3, 6, 10]) +end + +@testset "$NAME Min. Distance" begin + y = [0., 1., 3., 1., 4., 5., 3., 0., 2., 5., 4., 0.] + + @test Set(findpeaks(y, min_dist=4)[1]) == Set([6]) + @test Set(findpeaks(y, min_dist=3)[1]) == Set([6, 10]) + @test Set(findpeaks(y, min_dist=2)[1]) == Set([3, 6, 10]) + @test Set(findpeaks(y, min_dist=0)[1]) == Set([3, 6, 10]) + @test Set(findpeaks(y, min_dist=0)[1]) == Set([3, 6, 10]) +end + +@testset "$NAME Min. Height" begin + y = [0., 1., 3., 1., 4., 5., 3., 0., 2., 5., 4., 0.] + + @test Set(findpeaks(y, min_height=6.0)[1]) == Set([]) + @test Set(findpeaks(y, min_height=4.9)[1]) == Set([6, 10]) + @test Set(findpeaks(y, min_height=2.9)[1]) == Set([3, 6, 10]) +end + +@testset "$NAME Plateaus" begin + y = [0., 1., 3., 3., 3., 1., 4., 5., 3., 0., 0., 1., 5., 5., 4., 0.] + + @test Set(findpeaks(y)[1]) == Set([4, 8, 13]) + @test Set(findpeaks(y, min_plateau_points=1)[1]) == Set([4, 8, 13]) + @test Set(findpeaks(y, min_plateau_points=2)[1]) == Set([4, 13]) + @test Set(findpeaks(y, min_plateau_points=3)[1]) == Set([4]) + @test Set(findpeaks(y, min_plateau_points=4)[1]) == Set([]) + + @test Set(findpeaks(y, max_plateau_points=2)[1]) == Set([8, 13]) + @test Set(findpeaks(y, max_plateau_points=1)[1]) == Set([8]) + @test Set(findpeaks(y, max_plateau_points=0)[1]) == Set([]) + + @test Set(findpeaks(y, min_plateau_points=2, max_plateau_points=2)[1]) == Set([13]) +end + +@testset "$NAME Two Plateaus next to each other" begin + y = [0., 3., 3., 5., 5., 1.] + @test Set(findpeaks(y)[1]) == Set([4]) +end + +@testset "$NAME Threshold with height" begin + y = [0., 1., 3., 2., 4., 2., 1., 0.] + @test Set(findpeaks(y, min_height=2.9)[1]) == Set([3, 5]) + @test Set(findpeaks(y, min_height=2.9, threshold=2.)[1]) == Set([5]) +end + +@testset "$NAME PeakInfo" begin + y1 = [0., 1., 3., 1., 4., 5., 3., 0., 2., 5., 4., 0.] + + @test Set([ + PeakInfo(3., 2., 2., 2., 3:3), + PeakInfo(5., 1., 2., 5., 6:6), + PeakInfo(5., 3., 1., 5., 10:10)] + ) == Set(findpeaks(y1)[2]) + + y2 = [0., 1., 3., 3., 3., 1., 4., 5., 3., 0., 0., 1., 5., 5., 4., 0.] + @test Set([ + PeakInfo(3., 2., 2., 2., 3:5), + PeakInfo(5., 1., 2., 5., 8:8), + PeakInfo(5., 4., 1., 5., 13:14)] + ) == Set(findpeaks(y2)[2]) + + y3 = [0., 1., 3., 2., 4., 2., 1., 0.] + @test Set([ + PeakInfo(4., 2., 2., 4., 5:5) + ]) == Set(findpeaks(y3, min_height=2.9, threshold=2.)[2]) + + + # test prominence giving the correct value when it is limited by a boundary value + y4 = [1., 5., 0.] + @test Set([ + PeakInfo(5., 4., 5., 4., 2:2) + ]) == Set(findpeaks(y4)[2]) + +end + + +@testset "$NAME Empty inputs" begin + y1 = Float64[] + x1 = Int64[] + + @test isempty(findpeaks(y1)[1]) + @test findpeaks(y1, x1)[1] == empty(x1) + + y2 = Integer[] + x2 = String[] + + @test isempty(findpeaks(y2)[1]) + @test findpeaks(y2, x2)[1] == empty(x2) +end + +@testset "$NAME Non-equal data lengths" begin + y = [1.0, 2.0, 1.0] + x = [2.0] + + @test_throws ArgumentError findpeaks(y, x) +end + diff --git a/test/runtests.jl b/test/runtests.jl index 74af5b3fb..0bd014e07 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,7 @@ using Random: seed! testfiles = [ "dsp.jl", "util.jl", "windows.jl", "filter_conversion.jl", "filter_design.jl", "filter_response.jl", "filt.jl", "filt_stream.jl", "periodograms.jl", "resample.jl", "lpc.jl", "estimation.jl", "unwrap.jl", - "remez_fir.jl" ] + "remez_fir.jl", "findpeaks.jl" ] seed!(1776)