Differences

This shows you the differences between two versions of the page.

--- 전산물리학:압축_센싱 [2019/01/22 11:21] – [코드 구현] admin
+++ 전산물리학:압축_센싱 [2023/09/05 15:46] (current) – external edit 127.0.0.1
@@ Line 19: / Line 19: @@
 N = 256             # length of signal
 P = 5               # number of non-zero peaks
-K = 64              # number of measurements to take (N < L)
+K = 64              # number of measurements to take
 x = zeros(N,float)  # original signal (P-sparse)
@@ Line 44: / Line 44: @@
 실행 결과는 다음과 같다. 첫 번째 그림은 원래 신호를, 두 번째 그림은 복원된 신호와 원래 신호의 차이를 보여준다. 그 차이는 매우 작다.
-{{ ::전산물리학::sensing_x.png?400 |}}{{ ::전산물리학::sensing_d.png?400 |}}
+{{::전산물리학::sensing_x.png?400|}}{{::전산물리학::sensing_d.png?400|}}
 =====푸리에 공간에서 희박한 신호=====
+이 예는 참고문헌 항목 중 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$가 푸리에 변환의 역할을 해주는 행렬이다.
+코드는 다음처럼 작성할 수 있다.
+<Code:python>
+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()
+</Code>
+코드의 실행 결과는 아래 그림들처럼 나온다.
+{{::전산물리학::sensing_fourier1.png?400|}}
+{{::전산물리학::sensing_fourier2.png?400|}}
+첫 번째 그림은 원래 신호와 그 신호를 추출하기 위해서 무작위로 고른 위치들을 보여준다. 두 번째 그림은 원래 신호(실선)와 복원된 신호(점)를 보여주고 있다. 새넌-나이키스트 정리에서 요구하는 것보다 적은 수의 점을 추출했음에도 불구하고 복원된 신호는 원래 신호와 상당히 가깝다.
+=====통계물리적 접근=====
+년, L1 norm에 기반한 방법보다 더 적은 데이터를 가지고도 신호를 복원하는 방법이 개발되었다 [Krzakala et al., PRX 2, 021005 (2012)].
 ======참고문헌======
   * https://www.codeproject.com/Articles/852910/Compressed-Sensing-Intro-Tutorial-w-Matlab
   * http://www.pyrunner.com/weblog/2016/05/26/compressed-sensing-python/
   * https://dilawarnotes.wordpress.com/2017/09/06/compressed-sensing-a-python-demo/
+  * https://doi.org/10.1103/PhysRevX.2.021005