전산물리학:압축_센싱

개요

섀넌-나이키스트 표본화 정리(Shannon-Nyquist sampling theorem)에 의하면, 신호의 최대 주파수가 B일 때 이를 온전히 재구성하려면 2B 이상의 주파수로 표본을 추출해야 한다.그런데 우리가 보는 대부분의 신호는 큰 패턴이 중요한 역할을 하므로 주파수 밀도 측면에서 희박(sparse)하고, 이를 이용하면 섀넌-나이키스트 표본화 정리가 요구하는 것보다 적은 수의 표본으로도 신호를 재구성할 수 있다. 이러한 아이디어는 Candès, Romberg, and Tao (2006)에 의해 정식화되었고 압축 센싱(compressed sensing)이라 불리고 있다.

기본적인 설명

길이 $N$인 '희박한' 신호 $\mathbf{x}$가 미지수로서 있다. $k\times N$의 '관찰' 행렬 $A$를 곱하자. 이 행렬의 원소는 난수로 주어져 있고 $k \ll N$이라고 가정한다. 그 곱의 결과로 $k$ 차원 벡터 $\mathbf{y} = A \mathbf{x}$를 관찰 결과로 가지고 있다.

$A$와 $\mathbf{y}$만을 가지고 미지수 $\mathbf{x}$를 복원하는 것은 불가능해 보인다. $\mathbf{y}=A\mathbf{x}$를 만족하는 $\mathbf{x}$가 무수히 많기 때문이다. 하지만 그러한 $\mathbf{x}$ 중에서 가장 '희박한' 것, 구체적으로는 $\sum_i |x_i|$를 최소화하는 $\mathbf{x}$를 택하면 원래의 $\mathbf{x}$를 상당히 잘 복원해낼 수 있다.

아래는 참고문헌 항목 중 codeproject.com에 있는 Matlab 코드를 python으로 옮긴 것이다.

from __future__ import print_function,division
from numpy import zeros,dot,transpose
from numpy.random import normal,choice
from pylab import plot,show
from l1eq_pd import l1eq_pd

# Initialize constants and variables
N = 256             # length of signal
P = 5               # number of non-zero peaks
K = 64              # number of measurements to take
x = zeros(N,float)  # original signal (P-sparse)

# Generate signal with P randomly spread values
peaks = choice(N, P, replace=False)
x[peaks] = normal(0.0, 1.0, P)

# Obtain K measurements
A = normal(0.0, 1.0, (K,N))
y = dot(A,x)

# Perform Compressed Sensing recovery
x0 = dot(transpose(A), y)
xp = l1eq_pd(x0, A, [], y)
plot(x)
show()
plot(xp-x)
show()

신호의 길이는 $N=256$인데 이 중 $P=5$개 지점에서만 0이 아닌 값을 가진다. 이 신호는 그런 의미에서 '희박'하다. $K=64$는 관찰의 횟수로서 $N$의 1/4에 불과하다. choice 명령을 써서 5개의 점을 임의로 골라내고 그 곳들의 값을 정규(normal) 분포로부터 결정한다.

관찰 행렬 $A$도 정규 분포로부터 만든 다음, 관찰 결과인 $\mathbf{y}=A\mathbf{x}$를 만들어낸다. 이제 $\mathbf{x}$를 모른다고 가정한 채, $A$와 $\mathbf{y}$로부터 신호를 복원해낸다. 여기에서 절대값 합의 최소화를 위해 쓰인 l1eq_pd라는 함수는 아래 페이지에서 내려받을 수 있다: https://github.com/nikcleju/pyCSalgos

실행 결과는 다음과 같다. 첫 번째 그림은 원래 신호를, 두 번째 그림은 복원된 신호와 원래 신호의 차이를 보여준다. 그 차이는 매우 작다.

이 예는 참고문헌 항목 중 pyrunner.com의 내용으로부터 일부를 수정하여 가져온 것이다. 실공간 상에서는 조밀하지만 실제로는 몇 개의 주파수만이 중요한 신호를 생각해보자. 예컨대 다음과 같은 신호는 단 2개의 주파수 성분만을 지닌다: $$x = \sin(10 \pi t) + \sin (54 \pi t).$$ 이산 사인 변환(discrete sine transform, DCT)를 통해 이 신호를 푸리에 공간으로 옮기고 압축 센싱을 실행한 다음 역변환을 할 것이다. 그런데 앞의 코드를 재활용하기 위해 푸리에 변환을 행렬 곱 형태로 표현하면 편리하다 (계산 시간을 생각하면 비효율적인 선택일 것이다): $$\tilde{x} \equiv \mathcal{F}[x] = \Psi \mathbf{x},$$ $$x \equiv \mathcal{F}^{-1}[\tilde{x}] = \Psi^{-1} \tilde{\mathbf{x}}.$$ 이 떄 $x$는 원래 신호, $\tilde{x}$는 그것의 푸리에 변환($\mathcal{F}$) 결과이고 굵은 글씨로 적은 것은 신호를 벡터로서 적은 것을 의미한다. $\Psi$가 푸리에 변환의 역할을 해주는 행렬이다.

코드는 다음처럼 작성할 수 있다.

from __future__ import print_function,division
from math import pi
from numpy import sin,dot,transpose,identity,linspace
from numpy.random import choice
from pylab import plot,show
from l1eq_pd import l1eq_pd
from scipy.fftpack import dst,idst

# Initialize constants and variables
n = 600            # length of signal
m = 50             # number of measurements to take

# Generate signal with P randomly spread sinusoids
t = linspace(0,1,n)
x = sin(10*pi*t) + sin(54*pi*t)

# Orthonormal basis matrix
Psi = dst(identity(n), norm='ortho', axis=0)
Psi_inv = idst(identity(n), norm='ortho', axis=0)
xt = dot(Psi,x)

# Extract small sample of signal
ri = choice(n, m, replace=False)
ri.sort()
t2, x2 = t[ri], x[ri]
plot(t,x,'b')
plot(t2,x2, 'r.')
show()

A = Psi_inv[ri]      # sensing matrix
y = dot(A,xt)        # measured values

# Perform Compressed Sensing recovery
x0 = dot(transpose(A), y)
X_hat = l1eq_pd(x0, A, [], y)

x_rec = dot(Psi_inv,X_hat)     # IFFT of X_hat(f)
plot(x)
plot(x_rec, ".")
show()

코드의 실행 결과는 아래 그림들처럼 나온다.

첫 번째 그림은 원래 신호와 그 신호를 추출하기 위해서 무작위로 고른 위치들을 보여준다. 두 번째 그림은 원래 신호(실선)와 복원된 신호(점)를 보여주고 있다. 새넌-나이키스트 정리에서 요구하는 것보다 적은 수의 점을 추출했음에도 불구하고 복원된 신호는 원래 신호와 상당히 가깝다.

2012년, L1 norm에 기반한 방법보다 더 적은 데이터를 가지고도 신호를 복원하는 방법이 개발되었다 [Krzakala et al., PRX 2, 021005 (2012)].

참고문헌

  • 전산물리학/압축_센싱.txt
  • Last modified: 2021/12/17 19:52
  • by admin