압축센싱의 정식화는 단극음원을 기반으로 한 등가음원법(equivalent source method)으로부터 출발한다. 즉, 음장(sound field)을 여러 개의 단극음원에서 방사되는 파(wave)의 조합으로 근사한다. Fig. 1과 같이, 어떤 관심공간 내 음장이 N개의 등가음원으로 구성되어 있다고 가정하자. 이때 m(=1, 2, ..., M)번째 센서에서 계측되는 음압값 y_m은 그린 함수(green function)⁽¹⁶⁾를 통해 주파수 영역에서 식 (1)과 같이 나타낼 수 있다.

Fig. 1 
Distribution of N equivalent monopole sources in the interested domain. n-th equivalent source with strength x_n is located at each grid point. m-th sensor in the array receives the signal y_m

여기서 x_n [Pa⋅m]은 n번째 등가 음원의 강도를 나타내고, r_mn [m]은 n번째 등가 음원과 m번째 센서 사이의 거리이며, k(=2πf/c) [1/m]는 주파수 f [Hz]에 해당하는 파수를 나타낸다. 단, c=343 [m/s]는 음속을 의미한다. 식 (1)을 M개의 배열 센서에 대하여 확장하면, 관측벡터(measurement vector) y∈ℂ^M은 식 (2)와 같이 나타낼 수 있다.

여기서 A∈ℂ^M×N은 관측행렬(measurement matrix)이라고 한다.

그러나 실제 환경에서의 계측은 다양한 원인으로 유발된 잡음의 영향을 포함하기 때문에, 식 (2)를 일반화하면 식 (3)과 같다.

여기서 e=[e₁, ..., e_n]^T∈ℂ^M은 서로 상관관계가 없는 요소로 구성된 가우시안(gaussian) 잡음벡터이며, 위 첨자 T는 전치(transpose)를 의미한다.

일반적으로 관심 공간 내 등가 음원의 개수 N은 배열 센서의 개수 M보다 크기(M < N) 때문에. 식 (3)의 관측행렬 A는 미지수의 개수(N)가 식의 개수(M)보다 많은 부족결정 시스템(underdetermined system)을 구성하게 된다. 따라서 무한가지 경우의 해를 가지게 되나, 만약 실 음원의 개수 s가 등가 음원의 개수 N보다 매우 작다면(s ≪ N), 해 x∈ℂ^N는 대부분의 요소가 0 또는 0에 가까운 값을 가지는 희소벡터(sparse vector)로 볼 수 있으며, 식 (4)와 같은 l₀-norm 최소화 문제를 통해 해를 구할 수 있다.

여기서 l₀-norm ‖x‖₀은 0이 아닌 요소의 개수를 나타내며, l₂-norm은 ‖x‖₂=(∑x_n²)^1/2이고, ε은 잡음벡터 e의 l₂-norm을 나타낸다.

그러나 식 (4)의 최소화 문제는 NP-hard(Non-deterministic Polynomial-time hard) 문제에 해당하기에 사실상 계산이 불가능함이 알려져 있다⁽¹⁷⁾. 따라서 서론에서 언급한 바와 같이 볼록완화(convex relaxation)를 이용한 l₁-norm 최소화 문제(또는 기저추적이라고도 한다.)로 식 (4)를 근사하거나, OMP나 IHT 등을 통해 x의 지지집합을 반복 갱신하여 해를 추정하게 된다.

음원의 형태는 크게 3가지로 구분할 수 있으며, 모든 방향으로 균일한 음을 방사하는 점음원(point source), 무한개의 점음원이 직선의 형태로 분포하는 선음원(line source), 그리고 2차원 평면의 형태를 가지는 면음원(surface source) 등이 있다. 이 연구에서는 하나의 음원을 관심공간 내 하나의 격자점으로 표현 가능한 점음원에 대해서만 고려하기로 한다.

Fig. 4에는 제안 알고리듬의 검증을 위한 모사실험의 조건을 나타내었다. 센서 배열의 중심을 좌표계의 원점으로 설정하였고, 원점에서 z-축 방향으로 1.0 m 떨어진 1.0×1.0 m²의 영역을 5 mm 간격의 441개 격자로 분할하여 관심공간으로 지정하였다. 센서 배열에는 Fig. 5(a)와 같이 24개의 센서를 사용하여 랜덤한 배열을 구성하였고, 이로부터 M=24, N=441인 관측행렬 A를 생성하였다. 관측행렬 A의 RIP 조건 부합⁽⁷⁾에 대한 분석을 0 kHz ~ 12 kHz 주파수 대역에 사전 수행하였고, 분석 결과로부터 candes criteria⁽⁷⁾를 만족하는 주파수 대역은 5 kHz ~ 12 kHz로 나타났다. 이는 생성된 관측행렬 A, 다시 말해 현 실험조건에서 추정 가능한 음원의 주파수 대역이 5 kHz ~ 12 kHz임을 의미하며, 음원의 하나 또는 여러 주파수 성분이 해당 대역에 속하면 제안 알고리듬을 주파수 별로 반복 적용함으로써 음원위치추정이 가능함을 나타낸다. 여기서는 8 kHz의 순음(pure tone)을 음원의 주파수로 선정하였다.

Fig. 4 
Schematic drawing for the arrangement of interested domain and the sensor array

Fig. 5 
The random sensor array used in anti-aliasing low pass filter

Table 1에는 기존 IHT와 제안 알고리듬의 성능 비교를 위한 실험조건을 나타내었다. 음원 개수 s를 1에서 8까지 증가시키며 위치추정을 1000회씩 반복 수행하였으며, 이때 음원 위치는 각 시행마다 무작위로 지정하였고, 음원의 강도는 모든 시행에 걸쳐 동일하게 유지하였다.

Fig. 6에는 모사실험 결과를 나타내었다. 각 그림의 x-축은 음원 개수 s를 센서 개수 M으로 정규화화한 상대 음원수를 뜻하며, 제안 알고리듬(실선)과 기존 IHT⁽¹³⁾(점선)의 결과를 동시에 비교하였다. 신호 대 잡음비(signal-to-noise ratio, SNR)가 10 dB로 비교적 높은 환경(Fig. 6(a))과 0 dB로 낮은 환경(Fig. 6(b))에 대해, 그림의 좌측에는 반복계산 횟수의 평균을, 우측에는 복구 오차율 ψ의 평균을 도시하였다. 단, 복구 오차율 ψ는 식 (8)과 같다.

Fig. 6 
Comparison of the number of iterations (left) and recovery error (right) between the proposed- and previous algorithm

여기서 아래 첨자 t와 e는 각각 참값과 추정치를 의미한다.

제안 및 기존 알고리듬의 반복계산 횟수와 복구 오차율 ψ는 음원 개수 s가 증가할수록 비례하는 경향을 띠었다. 이는 음원이 많아질수록 음원위치추정에 필요한 반복계산 횟수, 즉 계산시간(혹은 계산비용)이 증가하며 해가 수렴하였어도 실제와는 다른 위치를 추정했을 가능성이 큼을 의미한다. 그러나 제안 알고리듬의 결과(실선)는 높은 SNR 환경(Fig. 6(a))에서 반복계산 횟수가 기존 IHT(점선) 대비 약 50 % 감소한 것을 알 수 있다. SNR과 무관하게 복구 오차율 ψ가 유사한 점을 미루어 동일 환경에서 제안 알고리듬은 기존보다 2배 빠르게 수렴한다고 볼 수 있으며, IHT 알고리듬의 계산시간 단축이 가능함을 다시 확인하였다. SNR이 상대적으로 열악한 환경(Fig. 6(b))에서는 알고리듬의 개선 효과가 반감되기는 하나, 여전히 기존 대비 성능 우위에 있음을 알 수 있다.

알고리듬 개선 효과는 상대 음원수(s/M)가 0.2 이상에서 급격하게 확대됨을 모사실험을 통해 볼 수 있었다. 이를 실제실험을 통해 검증하기 위해, 3(=s)개의 스피커 (SONY SRS-XB10)를 음원으로 사용하였고, 9(=M)개의 마이크로폰 유닛으로 구성된 센서 배열(Fig. 5(b))을 구성하였다. 각 유닛은 MEMS 마이크로폰(Knowles SPU0414HR5H-SB)과 고주파 잡음의 영향 억제를 위한 아날로그 안티 에일리어싱 필터(analogue anti-aliasing filter)의 조합으로 구성되어 있다(Fig. 5(c)).

모사실험의 경우와 마찬가지로 관측행렬 A의 RIP 조건에 대한 만족 여부를 사전 검토하였으며, candes criteria를 만족하는 주파수 대역 7 kHz ~ 12 kHz 중 8 kHz의 순음을 음원의 주파수로 선정하였다. Fig. 7에는 스피커 3개를 배치했을 때 각 배열 센서에서 계측된 신호의 평균 파워 스펙트럼을 도시하였고, 주파수 분해능 25 Hz, 해닝 창문함수(hanning window function), 75 %의 오버랩(overlapping)을 적용하여 계산하였다. 안티 에일리어싱 필터를 통해 고조파(harmonic) 성분이 어느 정도 억제된 것을 알 수 있으며, 8 kHz 부근에서 약 63 dB의 피크가 나타났다.

Fig. 7 
Averaged power spectrum of measurement signal

Table 2에는 실제실험의 조건을 나타내었다. 스피커의 음량을 조절하여 음원의 강도를 가급적 동일하게 유지하였으며, 모든 실험에 대해 샘플링 주파수 25.6 kHz로 신호를 수집하였다. 음원의 개수가 2개, 그리고 3개인 경우로 나누어 각각 30회씩 위치 추정을 반복하였으며, 이때 음원의 위치는 각 시행마다 임의로 배치하였다.

종료조건과 관련한 ε은 정규화 요소(regularization parameter)의 역할을 하며 그 값에 따라 추정결과가 크게 달라질 수 있다. 잡음의 l₂-norm으로부터 적절한 값을 선정함이 바람직하나, 모사실험과 달리 실제 환경에서는 잡음의 특성을 알 수 없다. 여기서는 시행착오적으로 여러 ε값을 시도한 결과 10^-6으로 결정하였다. 참고로, 제시된 값에서 상당 정도의 변화가 있더라도 추정결과에는 큰 변화가 없었다.

Table 3에는 음원 개수별로 두 성능지표에 대한 결과를 정리하였다. 제안 알고리듬의 평균 반복계산 횟수와 평균 계산시간은 기존 IHT 대비 음원의 개수가 2개일 때 각각 35.8 %, 24.7 %, 음원의 개수가 3개일 때 각각 27.6 %, 23.3 % 감소하였고, 모사실험의 결과로부터 예상할 수 있었듯 호모토피법의 스텝크기 결정방식을 이용한 제안법을 통해 IHT 알고리듬의 성능이 향상되었음을 알 수 있었다. 이때 평균 반복계산 횟수와 평균 계산시간의 감소율이 동일하지 않은 이유는 알고리듬이 한 번 반복될 때 소요되는 시간이 다르기 때문으로 추측된다.

Fig. 8에는 스피커 3개를 배치한 경우의 계측 신호를 바탕으로 계산된 반복계산 횟수의 비교 예시를 나타내었다. Fig. 2와 같이, 제안 알고리듬은 스텝크기를 큰 폭으로 변화시켜 해의 수렴을 촉진하는 반면, 기존 알고리듬에서는 도중부터 발생한 채터링 현상으로 해가 더디게 수렴함을 알 수 있다. 이를 통해 채터링 현상의 제거가 계산시간 저감에 직접적인 영향을 미쳤음을 확인할 수 있었다.

Fig. 8 
Comparison of the step size between the proposed- and previous algorithm in experiment

Fig. 9에는 Fig. 8에서의 위치추정결과를 나타내었다. 제안(Fig. 9(a)) 및 기존(Fig. 9(b)) 알고리듬 모두 음원의 위치, 다시 말해서 스피커의 위치를 잘 추정한 것을 알 수 있다. 음원의 추정 강도는 제안 알고리듬이 기존 대비 최대 2 dB 정도의 차이를 보였는데, 스피커의 음량을 가급적 동일하게 설정하였음을 고려하면 제안 알고리듬의 강도 추정 정확도가 상대적으로 뛰어나다고 할 수 있다.

Fig. 9 
Acoustic source map for algorithm