Preformance of CKMS+

AR(2) with 10 tap filter

In [17]:
using DSP, PyPlot, Polynomials, FFTW, Statistics, Suppressor

@suppress begin
sc = include("WFbyCKMSplus_scripts.jl")
at = include("../AnalysisToolbox.jl")
mrb = include("../WFMR_bs.jl")
mr  = include("../WFMR.jl")

whf = include("../WhiteningFilters.jl")
util = include("../Utils.jl")
end
Out[17]:
Main.VariousUtilities
In [8]:
N = 10^4; D = 10^3; p = 2; q = 0
e = randn(N+D)

Zeros = 2*rand(p) .- 1
Poles = 2*rand(q) .- 1

Zeros
Out[8]:
2-element Array{Float64,1}:
 0.844929538844911
 0.5974535445716462
In [10]:
l = coeffs(Polynomial([1])*prod(Polynomial([1, -z]) for z in Zeros))
w = [1] #coeffs(Polynomial([1])*prod(Polynomial([1, -z]) for z in Poles))
r = 1.0

pred = at.ARMA_gen(;l, w, r, e, steps = N, discard = D);

A = 2*rand(10) .- 1
f = coeffs(Polynomial([1])*prod(Polynomial([1, -a]) for a in A))
sig = filt(f,pred)

Nex = 2^13; Θ = 2pi*(0:Nex-1)/Nex;
Spred_num = at.z_crossspect_dm(pred,pred;Nex, L = 1500)
Hpred(z) = Polynomial(w)(z)/Polynomial(l)(z)
Spred_ana = map(z -> abs2(Hpred(z)),exp.(im*Θ))

semilogy(Θ,Spred_num, label = "Spred_num")
semilogy(Θ,Spred_ana, label = "Spred_ana")
legend()

Ssig_num = at.z_crossspect_dm(sig,sig;Nex)
Hsig(z) = Polynomial(f)(z)*Hpred(z)
Ssig_ana = map(z -> abs2(Hsig(z)),exp.(im*Θ))

figure()
semilogy(Θ,Ssig_num, label = "Ssig_num")
semilogy(Θ,Ssig_ana, label = "Ssig_ana")
legend()

sig  = reshape(sig, 1,:) 
pred = reshape(pred, 1,:)

[sig; pred]
Out[10]:
2×10000 Array{Complex{Float64},2}:
 -1.47492+0.0im  -2.06159+0.0im  …   -2.2798+0.0im   -1.22804+0.0im
 -1.47492+0.0im  -2.03207+0.0im     0.475345+0.0im  -0.415156+0.0im
In [11]:
function testdisp(h,sig,pred; vew = 90:200, f)
    figure(figsize=(12,4))
    title("The filters")
    f == 0 || plot(f,"k.-",label = "exact")
    plot(h[:],"r.:",label = "approx")
    xlabel("lag"); legend()

    sig_hat = at.my_filt(h,pred);

    vew = 90:200
    res = util.TakeLook(sig,sig_hat; vew)
    suptitle("A trajectory and Error over a Window")
    println("MSE: ",var(res[100:end]))
end
Out[11]:
testdisp (generic function with 1 method)

Backslash

In [12]:
# Benchmark
h = @timed mrb.get_wf_bs(sig,pred; M_out = 20)
tim = h.time; h = h.value

testdisp(h,sig,pred; f)

println("Time: ",tim," sec")

MSE: 1.0928494165248141e-30
Time: 1.814714545 sec

Old CKMS

In [13]:
# Benchmark
h = @timed mr.get_wf(sig,pred; M_out = 20)
tim = h.time; h = h.value

testdisp(h,sig,pred; f)

println("Time: ",tim," sec")

MSE: 2.050450707184427e-5
Time: 4.911224731 sec

CKMS+ (1 iteration)

In [18]:
# Benchmark
h = @timed sc.vector_wiener_filter_fft(sig,pred, maxit = 1)
tim = h.time; h = h.value

testdisp(h,sig,pred; f)

println("Time: ",tim," sec")

par = 1500
MSE: 
2.6486774783910133e-6
Time: 3.791689984 sec

CKMS+ (2 iterations)

In [19]:
# Benchmark
h = @timed sc.vector_wiener_filter_fft(sig,pred, maxit = 2)
tim = h.time; h = h.value

testdisp(h,sig,pred; f)

println("Time: ",tim," sec")

par = 1500
par = 1500

MSE: 5.419098218083698e-8
Time: 8.711294077 sec

CKMS+ (3 iterations)

In [22]:
# Benchmark
h = @timed sc.vector_wiener_filter_fft(sig,pred, maxit = 3, verb = true)
tim = h.time; h = h.value

testdisp(h,sig,pred; f)

println("Time: ",tim," sec")

Starting iterations
steps: 10000
par = 1500
steps: 8500
par = 1500
steps: 7000
par = 1500
Ending iterations
Time taken for crossspect: 17.401634227
Bytes Allocated: 42750091184
Time taken for spectfact: 1.996582589
Bytes Allocated: 3386962800

MSE: 0.00011543047100693456
Time: 19.4894837 sec

CKMS+ (4 iterations)

In [21]:
# Benchmark
h = @timed sc.vector_wiener_filter_fft(sig,pred, maxit = 4, verb = true)
tim = h.time; h = h.value

testdisp(h,sig,pred; f)

println("Time: ",tim," sec")

Starting iterations
steps: 10000
par = 1500
steps: 8500
par = 1500
steps: 7000
par = 1500
steps: 5500
par = 1500
Ending iterations
Time taken for crossspect: 30.077784239
Bytes Allocated: 67027797200
Time taken for spectfact: 12.723491174
Bytes Allocated: 21800027568

MSE: 3.2391207700903157
Time: 42.882749775 sec
In [ ]:
In [ ]: