[CGSP] Chap 12.4: Node Subsampling for PSD Estimation

using LinearAlgebra, Plots, FFTW, Statistics

columnwise_kron = 
(C,D) -> hcat([kron(C[:,i],D[:,i]) for i in 1:size(C)[2]]...)

#49 (generic function with 1 method)

12.4.1 The Sampling Problem

아래와 같이 길이가 \(N=10\) 인 신호 \({\bf x}\)를 고려하자.

x = rand(10)

10-element Vector{Float64}:
 0.03235208758206609
 0.5069925854414447
 0.5795228508497553
 0.682832351742401
 0.64422613488741
 0.24116013388795854
 0.8439116925218157
 0.6362602319916778
 0.386069828675059
 0.5313655894235898

여기에서 1,3,4,5 번째 원소만 추출하여길이가 \(K=4\) 인 신호 \({\bf y}\)를 만들고 싶다.

y = x[[1,3,4,5]]

4-element Vector{Float64}:
 0.03235208758206609
 0.5795228508497553
 0.682832351742401
 0.64422613488741

이 과정은 아래와 같이 수행할 수도 있다.

Φ= [1 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0]

4×10 Matrix{Int64}:
 1  0  0  0  0  0  0  0  0  0
 0  0  1  0  0  0  0  0  0  0
 0  0  0  1  0  0  0  0  0  0
 0  0  0  0  1  0  0  0  0  0

Φ*x

4-element Vector{Float64}:
 0.03235208758206609
 0.5795228508497553
 0.682832351742401
 0.64422613488741

즉 적당한 \(K\times N\) selection matrix를 선언하여 subsampling을 수행할 수 있다. 이때 매트릭스 \({\bf \Phi}\)를 subsampling matrix 혹은 sparse sampling matrix 라고 부른다.

12.4.2 Compressed LS Estimator

N = 10
V = [i*j for i in 0:(N-1) for j in 0:(N-1)] |> 
    x -> reshape(x,(N,N)) .|> 
    x -> exp(im * (2π/N) * x) 

10×10 Matrix{ComplexF64}:
 1.0+0.0im        1.0+0.0im          …        1.0+0.0im
 1.0+0.0im   0.809017+0.587785im         0.809017-0.587785im
 1.0+0.0im   0.309017+0.951057im         0.309017-0.951057im
 1.0+0.0im  -0.309017+0.951057im        -0.309017-0.951057im
 1.0+0.0im  -0.809017+0.587785im        -0.809017-0.587785im
 1.0+0.0im       -1.0+1.22465e-16im  …       -1.0+1.10218e-15im
 1.0+0.0im  -0.809017-0.587785im        -0.809017+0.587785im
 1.0+0.0im  -0.309017-0.951057im        -0.309017+0.951057im
 1.0+0.0im   0.309017-0.951057im         0.309017+0.951057im
 1.0+0.0im   0.809017-0.587785im         0.809017+0.587785im

G = columnwise_kron(conj(V),V)

100×10 Matrix{ComplexF64}:
 1.0+0.0im        1.0+0.0im          …        1.0+0.0im
 1.0+0.0im   0.809017+0.587785im         0.809017-0.587785im
 1.0+0.0im   0.309017+0.951057im         0.309017-0.951057im
 1.0+0.0im  -0.309017+0.951057im        -0.309017-0.951057im
 1.0+0.0im  -0.809017+0.587785im        -0.809017-0.587785im
 1.0+0.0im       -1.0+1.22465e-16im  …       -1.0+1.10218e-15im
 1.0+0.0im  -0.809017-0.587785im        -0.809017+0.587785im
 1.0+0.0im  -0.309017-0.951057im        -0.309017+0.951057im
 1.0+0.0im   0.309017-0.951057im         0.309017+0.951057im
 1.0+0.0im   0.809017-0.587785im         0.809017+0.587785im
 1.0+0.0im   0.809017-0.587785im     …   0.809017+0.587785im
 1.0+0.0im        1.0+0.0im                   1.0+0.0im
 1.0+0.0im   0.809017+0.587785im         0.809017-0.587785im
    ⋮                                ⋱  
 1.0+0.0im        1.0+0.0im                   1.0+0.0im
 1.0+0.0im   0.809017+0.587785im         0.809017-0.587785im
 1.0+0.0im   0.809017+0.587785im     …   0.809017-0.587785im
 1.0+0.0im   0.309017+0.951057im         0.309017-0.951057im
 1.0+0.0im  -0.309017+0.951057im        -0.309017-0.951057im
 1.0+0.0im  -0.809017+0.587785im        -0.809017-0.587785im
 1.0+0.0im       -1.0-1.11022e-16im          -1.0+2.27596e-15im
 1.0+0.0im  -0.809017-0.587785im     …  -0.809017+0.587785im
 1.0+0.0im  -0.309017-0.951057im        -0.309017+0.951057im
 1.0+0.0im   0.309017-0.951057im         0.309017+0.951057im
 1.0+0.0im   0.809017-0.587785im         0.809017+0.587785im
 1.0+0.0im        1.0+0.0im                   1.0+0.0im

- 방법1

ĉx = vec(x*x')
p̂ = inv(G' * G) * G' * ĉx

10-element Vector{ComplexF64}:
    0.25854107856772546 + 2.245922875954761e-20im
   0.004743491121735806 - 1.3138893409553828e-18im
   0.006946482731189413 - 9.791191432641327e-19im
   0.001721693617954179 - 1.9827974128203887e-18im
   0.011344167525098774 + 2.6827005818057562e-19im
 0.00012662617844242917 - 3.748573865136995e-20im
   0.011344167525098762 + 2.7448152053954017e-18im
  0.0017216936179541913 - 9.35534609073096e-19im
   0.006946482731189404 + 1.954408900185458e-18im
   0.004743491121735756 - 2.561030398375897e-18im

- 방법2

ĉy = vec(y*y')
p̂ = (kron(Φ,Φ)*G)' * ĉy

10-element Vector{ComplexF64}:
   3.759462826821233 + 0.0im
   2.765185174577697 - 2.0816681711721685e-17im
   1.077337414764992 + 2.7755575615628914e-17im
 0.11594812606807317 + 2.0816681711721685e-17im
 0.08838298603932843 + 3.903127820947816e-17im
 0.32863702713833354 + 4.622231866529366e-33im
 0.08838298603932859 + 9.540979117872439e-18im
  0.1159481260680729 - 2.0816681711721685e-17im
  1.0773374147649915 + 0.0im
  2.7651851745776965 - 2.0816681711721685e-17im