In [1]:
# Packages
using WAV
using LinearAlgebra
using FFTW
using SpecialMatrices
using LinearMaps
using Plots
using Arpack
In [2]:
s, Fs = wavread("files/piano_A41.wav",100001);
In [3]:
Fs
Out[3]:
In [3]:
# wavplay(s,Fs)
In [4]:
function monocomponents(x::Array{Float64}, xx::Array{Float64}, Fs::Float32, k::Int64, d::Int64)
# w is output from wavread(), package WAV, m is the number of samples we actually take
# x=vec(s) where s is the output of wavread()
# Fs is the output of wavread()
# k (even) is the number of eigenpairs computed by eigs()
# d is the level of recursion
# x=vec(s) where s is the output of wavread()
# There is no grouping of repeating components!
println("d = ",d)
if d<=2
H=Hankel(x)
n=size(H,1)
xsimple=Array{Float64}(undef,2*n-1)
# Fast Hankel EVD solver
f(x)=H*x
A=LinearMap(f,n,issymmetric=true)
λ,U=eigs(A, nev=k, which=:LM)
# @show λ
# Significant threshold percentage test
τ=0.1
L=map(Int,round((sum(abs.(λ).>(τ*maximum(abs.(λ))))/2)))
if L==1
@show L
# Return the signal as monocomponent.
# The signal is reconstructed using only ONE element of the skew-diagonals
xsimple=[(λ[1]*U[1,1])*U[:,1]; (λ[1]*U[n,1])*U[2:n,1]]
xsimple+=[(λ[2]*U[1,2])*U[:,2]; (λ[2]*U[n,2])*U[2:n,2]]
xx=[xx xsimple]
else
d+=1
# Analyze the i-th component
for i=1:L
# how L,i
l1=2*i-1
l2=2*i
xsimple=[(λ[l1]*U[1,l1])*U[:,l1]; (λ[l1]*U[n,l1])*U[2:n,l1]]
xsimple+=[(λ[l2]*U[1,l2])*U[:,l2]; (λ[l2]*U[n,l2])*U[2:n,l2]]
# Recursion - we do not need more than 2*L eigenpairs
xx=monocomponents(xsimple,xx, Fs,2*L,d)
end
end
end
xx
end
Out[4]:
In [5]:
x=vec(s)
xx=x
k=40
xmono=monocomponents(x,xx,Fs,k,0)
Out[5]:
We computed 8 mono-components (xmono[:,1] is the original signal).
Let us plot and listen to the mono-components.
In [6]:
norm(x-xmono[:,1])
Out[6]:
In [9]:
k=6
plot(xmono[:,k][1:400])
Out[9]:
In [14]:
t=range(0,stop=length(s)/Fs,length=length(s))
Out[14]:
In [30]:
l=10000
k=9
fs=abs.(fft(xmono[:,k]))
m,ind=findmax(fs[1:l])
println("Frequency = ", ind*Fs/length(fs) ," Hz, Amplitude = ", m)
plot(t[1:l],fs[1:l])
Out[30]:
In [9]:
wavplay(xmono[:,k],Fs)
In [10]:
xx=sum(xmono[:,2:9],dims=2)
wavplay(xx,Fs)
In [ ]: