using DynamicAxisWarping, DSP, LPVSpectral, Distances, Plots, Suppressor

N = 48_000
g(x,N) = exp(-10*(x-N/2)^2/N^2)
t = 1:N
f = range(0.01, stop=1, length=N)
y = sin.(t .* f) .* g.(t, N)
y1 = [y; 0y] .+ 0.1 .* randn.()
y2 = [0y; y] .+ 0.1 .* randn.()

plot([y1 y2][1:100:end,:])

M1,M2 = melspectrogram.((y1,y2), 2048)
@suppress plot(plot(M1), plot(M2))

dtw(.√(M1.power), .√(M2.power))[1]

@suppress dtwplot(.√(M1.power), .√(M2.power), linecolor=:green)

using SpectralDistances
n,m = size(M1.power)
dist = DiscreteGridTransportDistance(SqEuclidean(), n, n)
dtw(.√(M1.power), .√(M2.power), dist)[1]

@suppress dtwplot(.√(M1.power), .√(M2.power), dist, linecolor=:green)

f1 = range(0.01, stop=1, length=N)
f2 = range(0.1, stop=1, length=N).^2
y10 = sin.(t .* f1) .* g.(t, N)
y20 = sin.(t .* f2) .* g.(t, N)
y1 = [y10; 0y10] .+ 0.1 .* randn.()
y2 = [0y20; y20] .+ 0.1 .* randn.()
M1,M2 = melspectrogram.((y1,y2), 2048)
@suppress plot(plot(M1), plot(M2))

dtw(.√(M1.power), .√(M2.power))[1]

@suppress dtwplot(.√(M1.power), .√(M2.power), linecolor=:green)

dtw(.√(M1.power), .√(M2.power), dist)[1]

@suppress dtwplot(.√(M1.power), .√(M2.power), dist, linecolor=:green)

@btime dtw($(M1.power), $(M2.power), $dist)

size(M1.power)