전산물리학:압축_센싱

Differences

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

Link to this comparison view

Both sides previous revision Previous revision
Next revision		 Previous revision
전산물리학:압축_센싱 [2019/01/22 11:36] – [푸리에 공간에서 희박한 신호] admin전산물리학:압축_센싱 [2024/07/24 11:17] (current) – admin
Line 9: Line 9:
 =====코드 구현===== =====코드 구현=====
 아래는 참고문헌 항목 중 codeproject.com에 있는 Matlab 코드를 python으로 옮긴 것이다. 아래는 참고문헌 항목 중 codeproject.com에 있는 Matlab 코드를 python으로 옮긴 것이다.
-<Code:python>+<code:python>
 from __future__ import print_function,division from __future__ import print_function,division
 from numpy import zeros,dot,transpose from numpy import zeros,dot,transpose
Line 19: Line 19:
 N = 256             # length of signal N = 256             # length of signal
 P = 5               # number of non-zero peaks 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) x = zeros(N,float)  # original signal (P-sparse)
  
Line 37: Line 37:
 plot(xp-x) plot(xp-x)
 show() show()
-</Code>+</code>
 신호의 길이는 $N=256$인데 이 중 $P=5$개 지점에서만 0이 아닌 값을 가진다. 이 신호는 그런 의미에서 '희박'하다. $K=64$는 관찰의 횟수로서 $N$의 1/4에 불과하다. choice 명령을 써서 5개의 점을 임의로 골라내고 그 곳들의 값을 정규(normal) 분포로부터 결정한다. 신호의 길이는 $N=256$인데 이 중 $P=5$개 지점에서만 0이 아닌 값을 가진다. 이 신호는 그런 의미에서 '희박'하다. $K=64$는 관찰의 횟수로서 $N$의 1/4에 불과하다. choice 명령을 써서 5개의 점을 임의로 골라내고 그 곳들의 값을 정규(normal) 분포로부터 결정한다.
  
Line 44: Line 44:
  
 실행 결과는 다음과 같다. 첫 번째 그림은 원래 신호를, 두 번째 그림은 복원된 신호와 원래 신호의 차이를 보여준다. 그 차이는 매우 작다. 실행 결과는 다음과 같다. 첫 번째 그림은 원래 신호를, 두 번째 그림은 복원된 신호와 원래 신호의 차이를 보여준다. 그 차이는 매우 작다.
-{{ ::전산물리학::sensing_x.png?400 |}}{{ ::전산물리학::sensing_d.png?400 |}}+{{::전산물리학::sensing_x.png?400|}}{{::전산물리학::sensing_d.png?400|}}
  
 =====푸리에 공간에서 희박한 신호===== =====푸리에 공간에서 희박한 신호=====
 이 예는 참고문헌 항목 중 pyrunner.com의 내용으로부터 일부를 수정하여 가져온 것이다. 이 예는 참고문헌 항목 중 pyrunner.com의 내용으로부터 일부를 수정하여 가져온 것이다.
 실공간 상에서는 조밀하지만 실제로는 몇 개의 주파수만이 중요한 신호를 생각해보자. 예컨대 다음과 같은 신호는 단 2개의 주파수 성분만을 지닌다: 실공간 상에서는 조밀하지만 실제로는 몇 개의 주파수만이 중요한 신호를 생각해보자. 예컨대 다음과 같은 신호는 단 2개의 주파수 성분만을 지닌다:
-$$x = \sin(240 \pi t) + \sin (400 \pi t).$$ +$$x = \sin(10 \pi t) + \sin (54 \pi t).$$ 
-이산 코사인 변환(discrete cosine transform, DCT)를 통해 이 신호를 푸리에 공간으로 옮기고 압축 센싱을 실행한 다음 역변환을 할 것이다. 그런데 앞의 코드를 재활용하기 위해 푸리에 변환을 행렬 곱 형태로 표현하면 편리하다 (단, 계산 시간을 생각하면 비효율적인 선택일 것이다):+이산 사인 변환(discrete sine transform, DCT)를 통해 이 신호를 푸리에 공간으로 옮기고 압축 센싱을 실행한 다음 역변환을 할 것이다. 그런데 앞의 코드를 재활용하기 위해 푸리에 변환을 행렬 곱 형태로 표현하면 편리하다 (계산 시간을 생각하면 비효율적인 선택일 것이다):
 $$\tilde{x} \equiv \mathcal{F}[x] = \Psi \mathbf{x},$$ $$\tilde{x} \equiv \mathcal{F}[x] = \Psi \mathbf{x},$$
 $$x \equiv \mathcal{F}^{-1}[\tilde{x}] = \Psi^{-1} \tilde{\mathbf{x}}.$$ $$x \equiv \mathcal{F}^{-1}[\tilde{x}] = \Psi^{-1} \tilde{\mathbf{x}}.$$
Line 56: Line 56:
  
 코드는 다음처럼 작성할 수 있다. 코드는 다음처럼 작성할 수 있다.
-<Code:python>+<code:python>
 from __future__ import print_function,division from __future__ import print_function,division
 from math import pi from math import pi
Line 63: Line 63:
 from pylab import plot,show from pylab import plot,show
 from l1eq_pd import l1eq_pd from l1eq_pd import l1eq_pd
-from scipy.fftpack import dct,idct+from scipy.fftpack import dst,idst
  
 # Initialize constants and variables # Initialize constants and variables
 n = 600            # length of signal n = 600            # length of signal
-m = 150            # number of measurements to take (N < L)+m = 50             # number of measurements to take
  
 # Generate signal with P randomly spread sinusoids # Generate signal with P randomly spread sinusoids
-t = linspace(0,1/16,n) +t = linspace(0,1,n) 
-x = sin(240*pi*t) + sin(400*pi*t)+x = sin(10*pi*t) + sin(54*pi*t)
  
 # Orthonormal basis matrix # Orthonormal basis matrix
-Psi = dct(identity(n), norm='ortho', axis=0) +Psi = dst(identity(n), norm='ortho', axis=0) 
-Psi_inv = idct(identity(n), norm='ortho', axis=0)+Psi_inv = idst(identity(n), norm='ortho', axis=0)
 xt = dot(Psi,x) xt = dot(Psi,x)
  
Line 95: Line 95:
 x_rec = dot(Psi_inv,X_hat)     # IFFT of X_hat(f) x_rec = dot(Psi_inv,X_hat)     # IFFT of X_hat(f)
 plot(x) plot(x)
-plot(x_rec) +plot(x_rec, ".")
-plot(x-x_rec)+
 show() show()
-</Code> +</code> 
-실행 결과는 아래와 같다.+코드의 실행 결과는 아래 그림들처럼 나온다. 
 + 
 +{{::전산물리학::sensing_fourier1.png?400|}} 
 +{{::전산물리학::sensing_fourier2.png?400|}} 
 + 
 +첫 번째 그림은 원래 신호와 그 신호를 추출하기 위해서 무작위로 고른 위치들을 보여준다. 두 번째 그림은 원래 신호(실선)와 복원된 신호(점)를 보여주고 있다. 새넌-나이키스트 정리에서 요구하는 것보다 적은 수의 점을 추출했음에도 불구하고 복원된 신호는 원래 신호와 상당히 가깝다. 
 + 
 +=====통계물리적 접근===== 
 + 
 +2012년, L1 norm에 기반한 방법보다 더 적은 데이터를 가지고도 신호를 복원하는 방법이 개발되었다 [Krzakala et al., PRX 2, 021005 (2012)]. 
  
 ======참고문헌====== ======참고문헌======
Line 105: Line 114:
   * http://www.pyrunner.com/weblog/2016/05/26/compressed-sensing-python/   * http://www.pyrunner.com/weblog/2016/05/26/compressed-sensing-python/
   * https://dilawarnotes.wordpress.com/2017/09/06/compressed-sensing-a-python-demo/   * https://dilawarnotes.wordpress.com/2017/09/06/compressed-sensing-a-python-demo/
 +  * https://doi.org/10.1103/PhysRevX.2.021005
  • 전산물리학/압축_센싱.1548124563.txt.gz
  • Last modified: 2023/09/05 15:46
  • (external edit)