Preformance of CKMS+

ARMA(10,10) with 10 tap filter

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

@suppress_out 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[109]:
Main.VariousUtilities
In [43]:
Zeros = 2*rand(10) .- 1
Poles = 2*rand(10) .- 1

[Zeros Poles]
Out[43]:
10-element Array{Float64,1}:
 -0.9533171783748062
  0.6006760516695353
 -0.9471600026794325
  0.027280322740684593
 -0.05904056784962419
 -0.626624700050439
  0.07397089978187665
  0.00973589339259906
  0.8189080074717072
  0.7819765781164727
In [83]:
l = coeffs(Polynomial([1])*prod(Polynomial([1, -z]) for z in Zeros))
w = coeffs(Polynomial([1])*prod(Polynomial([1, -z]) for z in Poles))
r = 1.0

N = 10^4; D = 10^3; p = 2;
e = randn(N+D)

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[83]:
2×10000 Array{Complex{Float64},2}:
 65.6576+0.0im  -123.915+0.0im  …  -3.30248+0.0im  -1.65636+0.0im
 65.6576+0.0im   62.5172+0.0im      232.188+0.0im   199.595+0.0im
In [98]:
function testdisp(h,sig,pred; vew = 90:200, f)
    @suppress begin
    figure(figsize=(12,4))
    title("The filters")
    f == 0 || plot(f,"k.-",label = "exact")
    plot(real(h[:]),"r.:",label = "approx")
    xlabel("lag"); legend()
    end
    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[98]:
testdisp (generic function with 1 method)

Backslash

In [99]:
# 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.1581138833115986e-26
Time: 0.006065444 sec

Old CKMS

In [100]:
# 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: 1.492989912142938
Time: 0.313833642 sec

CKMS+ (1 iteration)

In [101]:
# 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")

Starting iterations
steps: 10000
Ending iterations
Time taken for crossspect: 2.484085618
Bytes Allocated: 8750431200
Time taken for spectfact: 0.168156353
Bytes Allocated: 262391520

MSE: 1.4988421297939916
Time: 2.726050193 sec

CKMS+ (2 iterations)

In [105]:
# 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")

Starting iterations
steps: 10000
steps: 7000
Ending iterations
Time taken for crossspect: 8.83260085
Bytes Allocated: 23151259152
Time taken for spectfact: 0.327119022
Bytes Allocated: 529402320

MSE: 0.01197143634248616
Time: 9.234497356 sec

CKMS+ (3 iterations)

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

testdisp(h,sig,pred; f)

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

Starting iterations
steps: 10000
par = 1000
steps: 7000
par = 1000
steps: 4000
par = 1000
Ending iterations
Time taken for crossspect: 12.961039386
Bytes Allocated: 27216246806
Time taken for spectfact: 9.059043232
Bytes Allocated: 14667153109

MSE: 10.782586111812499
Time: 22.16997829 sec