Preformance of CKMS+

I keep all my scripts out of the way in a file in this directory.

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

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

(.:199411): Gdk-CRITICAL **: 09:55:21.200: gdk_cursor_new_for_display: assertion 'GDK_IS_DISPLAY (display)' failed

(.:199411): Gdk-CRITICAL **: 09:55:21.204: gdk_cursor_new_for_display: assertion 'GDK_IS_DISPLAY (display)' failed
Out[1]:
Main.VariousUtilities

As a working test case I pick the AR(2) process $y_n = 5/4y_{n-1} - 3/8y_{n-2} + e_n$ as signal to be estimated and $x_n = y_n + u_n$

In [28]:
N = 10^4; D = 10^3; p = 2;
e = randn(N+D)

l = [1, -5/4, 3/8]
w = [1]
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^12; Θ = 2pi*(0:Nex-1)/Nex;
Spred_num = at.z_crossspect_dm(pred,pred;Nex, L = 150)
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[28]:
2×10000 Array{Complex{Float64},2}:
 -0.827652+0.0im  1.90188+0.0im  …  -0.409354+0.0im  -1.04874+0.0im
 -0.827652+0.0im  1.37657+0.0im        2.9424+0.0im   1.61986+0.0im
In [2]:
sig, pred = sc.KSEdata(10^5; pred_ind = 1:1, sig_ind= 1:1, geni = 1)

[sig; pred]

Sol load location: ../../../data/KSE_Data/ks_pred_lin1e5_r1.jld
Out[2]:
2×100000 Array{Complex{Float64},2}:
 -0.0622791-0.0665026im  -0.0602835-0.0630604im  …  -0.191963+0.0385498im
 -0.0643694-0.0698753im  -0.0622791-0.0665026im     -0.188049+0.0347455im
In [29]:
# Benchmark
h = mrb.get_wf_bs(sig,pred; M_out = 20)

# plot(f,"k.-",label = "exact")
plot(h[:],"r.:",label = "approx")
xlabel("lag")
Out[29]:
PyObject Text(0.5, 24.0, 'lag')
In [30]:
sig_hat = at.my_filt(h,pred);

vew = 90:200
res = util.TakeLook(sig,sig_hat; vew)

var(res[100:end])
Out[30]:
1.0561665295558965e-30

Now using old CKMS

In [31]:
# Benchmark
h = mr.get_wf(sig,pred; M_out = 20)

# plot(f,"k.-",label = "exact")
plot(h[:],"r.:",label = "approx")
xlabel("lag")
Out[31]:
PyObject Text(0.5, 24.0, 'lag')
In [32]:
sig_hat = at.my_filt(h,pred);

vew = 90:2000
res = util.TakeLook(sig,sig_hat; vew)

var(res[100:end])
Out[32]:
2.8613959268596034e-6

Now we look at new CKMS.

In [33]:
h_num_fft = sc.vector_wiener_filter_fft(sig,pred, maxit = 2)

# plot(f,"k.-",label = "exact")
plot(h_num_fft[:],"r.:",label = "approx")
xlabel("lag")

Starting iterations
steps10000
steps7000
Ending iterations
Time taken for crossspect: 8.860561581
Bytes Allocated: 23151258416

Time taken for spectfact: 0.192028765
Bytes Allocated: 348270912
Out[33]:
PyObject Text(0.5, 24.0, 'lag')
In [34]:
sig_hat = at.my_filt(h_num_fft,pred);

res = util.TakeLook(sig,sig_hat; vew)

var(res[100:end])
Out[34]:
6.136185938451246e-9