Preformance of CKMS+

VAR(2) with 10 tap filter

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

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")

Unable to init server: Could not connect: Connection refused
Unable to init server: Could not connect: Connection refused

(.:207289): Gdk-CRITICAL **: 14:25:15.825: gdk_cursor_new_for_display: assertion 'GDK_IS_DISPLAY (display)' failed

(.:207289): Gdk-CRITICAL **: 14:25:15.828: gdk_cursor_new_for_display: assertion 'GDK_IS_DISPLAY (display)' failed
Out[1]:
Main.VariousUtilities
In [2]:
N = 10^4; D = 10^3; p = 2; q = 0
e1 = randn(N+D)
e2 = randn(N+D)

Zeros1 = 2*rand(p) .- 1
Zeros2 = 2*rand(p) .- 1
Poles = 2*rand(q) .- 1;
In [7]:
l1 = coeffs(Polynomial([1])*prod(Polynomial([1, -z]) for z in Zeros1))
l2 = coeffs(Polynomial([1])*prod(Polynomial([1, -z]) for z in Zeros2))
w = [1] #coeffs(Polynomial([1])*prod(Polynomial([1, -z]) for z in Poles))
r = 1.0

pred1 = at.ARMA_gen(; l = l1, w, r, e = e1, steps = N, discard = D);
pred2 = at.ARMA_gen(; l = l2, w, r, e = e2, steps = N, discard = D);

pred1 = reshape(pred1, 1,:)
pred2 = reshape(pred2, 1,:)

pred = [pred1; pred2]

A1 = 2*rand(10) .- 1
A2 = 2*rand(10) .- 1

f1 = coeffs(Polynomial([1])*prod(Polynomial([1, -a]) for a in A1))
f2 = coeffs(Polynomial([1])*prod(Polynomial([1, -a]) for a in A2))

f = zeros(1,2,11); f[1,1,:] = f1; f[1,2,:] = f2; 
sig = at.my_filt(f,pred)
Out[7]:
1×10000 Array{Complex{Float64},2}:
 2.06109+0.0im  3.49735+0.0im  …  0.000212496+0.0im  -5.54882+0.0im
In [9]:
Nex = 2^13; Θ = 2pi*(0:Nex-1)/Nex;

Spred1_num = at.z_crossspect_dm(pred1,pred1;Nex, L = 1500)
Hpred1(z) = Polynomial(w)(z)/Polynomial(l1)(z)
Spred1_ana = map(z -> abs2(Hpred1(z)),exp.(im*Θ))

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

figure()
Spred2_num = at.z_crossspect_dm(pred2,pred2;Nex, L = 1500)
Hpred2(z) = Polynomial(w)(z)/Polynomial(l2)(z)
Spred2_ana = map(z -> abs2(Hpred2(z)),exp.(im*Θ))

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


Ssig_num = at.z_crossspect_dm(sig,sig;Nex)
Hsig(z) = Polynomial(f1)(z)*Hpred1(z) + Polynomial(f2)(z)*Hpred2(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; pred]
Out[9]:
3×10000 Array{Complex{Float64},2}:
  2.06109+0.0im   3.49735+0.0im  …  0.000212496+0.0im  -5.54882+0.0im
  1.37699+0.0im  0.801088+0.0im         1.35416+0.0im  0.408652+0.0im
 0.684102+0.0im  0.126085+0.0im       -0.203989+0.0im   1.67233+0.0im
In [11]:
function testdisp(h,sig,pred; vew = 90:200, f)
    figure(figsize=(12,4))
    title("The filters")
    if f!= 0
        plot(f[1,1,:],"k.-",label = "exact1")
        plot(h[1,1,:],"r.:",label = "approx1")
        plot(f[1,2,:],"k.-",label = "exact2")
        plot(h[1,2,:],"r.:",label = "approx2")
        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[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: 2.0093432440332154e-29
Time: 1.836755451 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: 0.002625931516744872
Time: 5.377012438 sec

CKMS+ (1 iteration)

In [14]:
# 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.36277273165177e-5
Time: 7.427863306 sec

CKMS+ (2 iterations)

In [15]:
# 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: 2.307980013582431e-6
Time: 15.662041842 sec

CKMS+ (3 iterations)

In [16]:
# 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: 29.796485332
Bytes Allocated: 80850964743
Time taken for spectfact: 30.56945486
Bytes Allocated: 86645927248

MSE: 36.37196324033341
Time: 60.499290498 sec

CKMS+ (4 iterations)

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