spectrogramf
elapsed time: 0.017406119 seconds (6052416 bytes allocated)
Array{Float64,1}
4096
spectrogramf2
elapsed time: 0.002346886 seconds (2663856 bytes allocated)
Array{Float64,1}
4096
spectrogram
elapsed time: 0.084957317 seconds (11150800 bytes allocated)
Array{Float64,2}
4096
welch_pgram
elapsed time: 0.040620592 seconds (7973488 bytes allocated, 55.53% gc time)
Array{Float64,1}
4096
using DSP, Base.Test
using Base.FFTW: fftwNumber
typealias dspNumber fftwNumber
function spectrogramf{T<:dspNumber}(s::Array{T,1}; n=int(length(s)/8), m=int(n/2), r=1, w=(n)->ones(n))
w = w(n)
ind = arrayspliti(s, n, m)
sf = zeros(eltype(s),n)
for i=ind
sf += periodogramf(s[i:i+n-1].*w)
end
t=( (0:length(ind)-1)*(n-m) + n/2) / r
f=(0:n-1)/n*r
sf/(r*length(ind)), t, f
end
function spectrogramf2{T<:dspNumber}(s::Array{T,1}; n=int(length(s)/8), m=int(n/2), r=1, w=(n)->ones(n))
w = w(n)
ind = arrayspliti(s, n, m)
sf = zeros(eltype(s),n)
per = zeros(typeof(1.0im),n)
for i=ind
per[:] = s[i:i+n-1].*w
periodogramtl!(per)
sf += real(per)
end
t=( (0:length(ind)-1)*(n-m) + n/2) / r
f=(0:n-1)/n*r
sf/(r*length(ind)*length(per)), t, f
end
function periodogramf{T<:dspNumber}(s::Array{T,1})
s_fft = fft(s)
abs2(s_fft)/length(s_fft)
end
function periodogramtl!{T<:dspNumber}(s::Array{T,1})
fft!(s)
s[:] = abs2(s)
return nothing
end
function arraysplit_su(sn::Integer, n::Integer, m::Integer)
# n = m is a problem - the algorithm will not terminate.
if !(0 <= m < n)
error("m must be between zero and n.")
end
# the length of the non-overlapping array stride
l = n - m
# total number of strides is the total length of the signal divided
# by the unique number of elements per stride. extra elements at the
# end of of the signal are dropped.
k = ifloor(sn/l - n/l + 1)
l, k
end
function arrayspliti{T<:dspNumber}(s::Array{T,1}, n::Integer, m::Integer)
l, k = arraysplit_su(length(s), n, m)
[a*l + 1 for a=0:(k-1)]
end
# test
N=1024*32
x0=randn(N)
nn=int(N/8)
mm=int(N/16)
# initialize variables before timing
sf, tf, ff = spectrogramf(x0, n=nn, r=10, m=mm)
sf2, tf2, ff2 = spectrogramf2(x0, n=nn, r=10, m=mm)
p, t, f = spectrogram(x0, n=nn, r=10, m=mm)
sw= welch_pgram(x0, nn, mm)
println(spectrogramf)
@time begin sf, tf, ff = spectrogramf(x0, n=nn, r=10, m=mm)
end
println(typeof(sf))
println(length(sf))
println(spectrogramf2)
@time begin sf2, tf2, ff2 = spectrogramf2(x0, n=nn, r=10, m=mm)
end
println(typeof(sf))
println(length(sf))
println(spectrogram)
@time begin p, t, f = spectrogram(x0, n=nn, r=10, m=mm)
s = 1/size(p,2)*sum(p,2)
end
println(typeof(s))
println(length(s))
println(welch_pgram)
@time begin sw= welch_pgram(x0, nn, mm)
end
println(typeof(sw))
println(length(sw))
@test_approx_eq s sf2
@test_approx_eq s sf
@test_approx_eq t tf
@test_approx_eq f ff
@test_approx_eq s sw/10