Transactions of the Korean Society for Noise and Vibration Engineering
[ Article ]
Transactions of the Korean Society for Noise and Vibration Engineering - Vol. 35, No. 4, pp.395-408
ISSN: 1598-2785 (Print) 2287-5476 (Online)
Print publication date 20 Aug 2025
Received 03 Jul 2025 Revised 24 Jul 2025 Accepted 24 Jul 2025
DOI: https://doi.org/10.5050/KSNVE.2025.35.4.395

일반화된 주파수 합 빔형성을 이용한 저주파 대역 도래각 추정

이정훈 ; 박용성* ; 피터 거스토프**
Low-frequency Direction-of-arrival Estimation using Generalized Frequency-sum Beamforming
Jeung-Hoon Lee ; Yongsung Park* ; Peter Gerstoft**
*Applied Ocean Physics and Engineering, Woods Hole Oceanographic Institution, Research Scientist
**Scripps Institution of Oceanography, University of California, Professor

Correspondence to: Member, Dept. of Mechanical Engineering, Changwon National University, Professor E-mail : jhoonlee@changwon.ac.krRecommended by Editor Hanshin Seol

Ⓒ The Korean Society for Noise and Vibration Engineering

Abstract

Estimating the direction-of-arrival (DOA) at low frequencies is challenging owing to the inadequate spatial resolution of conventional delay-and-sum beamforming. This study presents an alternative approach based on generalized frequency-sum (gFS) processing, which constructs higher-order spectral autoproducts up to a selected order Q. By combining frequency components, the method enhances spatial resolution when the acoustic wavelength exceeds the array aperture. The upper limit of Q is governed by the condition that the resulting frequency sum remains below the spatial Nyquist limit; in practice, this bound is further limited by noise. In addition, through a framework based on multinomial expansion, we demonstrate that gFS is inherently unsuitable for resolving multiple DOAs. Numerical simulations and experiments confirm that this technique offers a viable solution for single-source DOA estimation in the low-frequency regime.

Keywords:

Frequency-sum, Beamforming, Autoproduct

키워드:

합주파수, 빔형성, 자기곱

1. 서 론

센서배열을 이용한 도래각(direction of arrival, DOA) 추정은 레이더, 소나, 무선통신 등 다양한 신호처리 분야에서 활발히 연구되고 있는 주제이다(1). 도래각 추정의 대표적인 방법은 지연-합(delay-and-sum) 기반의 빔포밍(conventional beamforming, CBF) 방법이 보편적으로 사용된다(2,3). 구현이 간단하고 직관적이지만, 저주파 대역에서는 해상도(공간 분해능)이 부족하여 적용이 어렵다. 저주파 신호는 아파트 층간 충격음(4), 고래 울음소리(5), 공동현상(cavitation)에 의한 프로펠러 소음 등 다양한 상황에서 발생한다(6~8).

적응형 기법인 MVDR(minimum variance distortionless response)과 MUSIC(MUltiple SIgnal Classification)은 CBF의 해상도 한계를 개선한다는 장점이 있으나, 배경잡음이 낮고, 신호원이 비상관적이며, 스냅샷(snapshot)이 충분한 경우에 한해 효과적이다(9,10). 압축센싱(compressive sensing, CS) 기법은 고해상도 구현이 가능하지만, 저주파 영역에서는 관측행렬의 기저벡터들이 서로 상관성을 가지게 되어 도래각 추정 정확도를 저해한다(11,12).

Abadi et al.은 광대역 신호를 가정하여 주파수합 빔포밍(frequency-sum beamforming, FS)을 제안한 바 있다(13). 이 방법은 동일한 공간 위치에서 대역 내 두 주파수 성분을 곱하는 2차 주파수합 자기곱 (frequency-sum autoproduct)을 활용한다. 자기곱 개념은 주파수 ‘차’ 빔포밍(frequency-difference beamforming, FD)과 함께 제안되었으며(14), FD는 낮은 주파수 성분의 켤레복소를 높은 주파수 성분과 곱하여 주파수차 자기곱을 생성한다(15~17). 주파수차 자기곱은 그것을 구성하는데 사용된 두 성분의 차주파수, 즉, 낮은 주파수에서의 음장을 모사한다는 것이 밝혀졌다(18). 따라서, FD는 대역 내 주파수보다 훨씬 낮은 주파수에서 빔포밍을 수행하며, 이는 센서 간격의 두 배보다 짧은 파장을 가지는 고주파 신호의 도래각을 추정하는데 적합하다(18~20).

FS는 FD와 반대 개념으로 이해할 수 있다. 스펙트럼 성분 두 개의 곱으로 정의되는 주파수합 자기곱은 합주파수에서의 음장을 모사할 수 있다. FS는 합주파수에서 빔포밍을 수행하므로 빔폭(beamwidth)을 좁힐수 있으며, 이는 CBF보다 높은 공간 해상도를 제공한다. 따라서, 배열 크기보다 더 큰 파장을 가지는 저주파 신호에 대해 효과적이다.

이 논문에서 제안하는 Q차 주파수합 자기곱을 이용한 일반화된 주파수합 빔포밍(generalized frequency-sum beamforming, gFS)은 기존 FS의 해상도 한계를 극복하기 위함이다. 이 고차 자기곱은 대역 내 Q개 주파수 성분을 곱하여 구성되며, 자기곱 구성에 사용된 주파수를 모두 더하여 합주파수를 얻는다. 높은 차수의 정보를 높은 합주파수에서 처리함으로써 해상도 향상을 기대할 수 있다. Q차 자기곱은 단순히 Q개의 성분을 곱하는 방식이므로, 계산 복잡도는 기존 FS와 동일하게 유지된다. 따라서, 이 방법은 간단한 구현 대비 해상도 향상의 효용이 크다.

또한, 고차 자기곱은 해상도 개선을 위한 기존의 알고리즘과도 쉽게 결합될 수 있다는 측면에서도 장점을 갖는다. 예를 들어, 저자들은 주파수차 자기곱을 압축센싱 및 MUSIC 알고리즘에 적용함으로써, 근접한 두 음원에 대한 분리가 가능함을 보여준 바 있다(21,22). Abadi et al.은 주파수합 자기곱을 MVDR 알고리즘과 결합하여 공동현상 영상화(passive cavitation imaging)에 적용한 바 있다(13). 유사한 방식으로, 고차 자기곱 역시 CS, MUSIC, MVDR 등의 고해상도 알고리즘과도 결합이 가능하며, 이를 통해 추가적인 성능 개선이 기대된다. 이러한 확장 방향은 이 논문의 범위를 벗어나므로, 여기에서는 gFS 본연의 특징과 장-단점에 대해 집중한다.

gFS는 다중 도래각(multi-DOA) 신호뿐만 아니라, 잡음이 큰 환경에서 성능 저하를 보인다. 자기곱이 유발하는 교차항(crossterm)에 의해 허위 도래각(ghost DOA)이 발생한다. FD 기법에서는 주파수 평균을 통해 허위 성분을 억제할 수 있었지만, 주파수 합 기반의 gFS에서는 그렇지 않다. 이 논문에서는 다항전개(multinomial expansion) 이론을 통해 그 원인을 해석적으로 규명하고자 한다(23).

논문 구성은 다음과 같다. 2절에서는 gFS를 제안하고, 다항 전개를 활용하여 장-단점에 대해 분석한다. 3절에서는 수치 시뮬레이션과 실험을 통해 gFS 성능을 검증하고, 4절에서 결론을 맺는다.


2. 이 론

gFS는 다중 도래각 상황에서 적용이 어려우나(2.2절), 일반화를 위해 K개 평면파(plane waves) 음원을 가정한다. 각 음파는 간격 d, M개의 센서로 구성된 균일 선배열(ULA: uniform linear array)에 입사하고 있다. k(=1, …, K)번째 음원은 도래각 θk∈[-90°, 90°](0°는 배열의 수직 방향)와 복소진폭 xk로 표현된다. 자기곱을 형성하기 위해 음원의 광대역성도 가정에 포함된다. 주파수 범위는 [fL,fU](f의 단위는 Hz임)이나, 저주파 음원을 고려하는 상황이므로, 상한 주파수 fU는 배열의 공간 나이퀴스트 주파수 fN=c/2d (c는 음속, [m/s])보다 작다는 가정도 포함한다. 주파수 구성은 반드시 연속적일 필요는 없으며, 고조파(harmonic)와 같은 이산 대역도 허용된다.

k번째 음원과 센서 사이의 위상 지연은 식 (1)의 조향벡터로 표현된다.

aθk,f=1M1e-j2πf/c1d sinθke-j2πf/cM-1 d sinθkTM(1) 

상첨자 T는 전치(transpose)를, j=-1,MM차원 복소 열벡터를 의미하고, 단순한 해석을 위해 M으로 정규화하였다. 배열출력 y∈ℂM식 (2)와 같은데, n∈ℂM은 가우시안 잡음(Gaussian noise)벡터이다.

y=k=1KaΘk,fxk+n(2) 

2.1 일반화된 주파수 합 빔형성(gFS)

Q차 주파수합 자기곱(이하, 일반화된 자기곱이라 칭함) gQ∈ℂM식 (3)과 같이 정의하자.

gQ=yf1yfQ=q=1Qyfq(3) 

여기서 ⊙기호는 성분별 곱(Hadamard 곱)을 의미한다. gQ 표시에 있어 누적곱기호 ∏를 이용하였는데, 이는 q = 1, …, Q에 해당하는 벡터 y(fq)들에 대한 Hadamard 곱을 의미한다. 일반적으로 ∏ 기호는 스칼라에 한하지만, 간결한 표현을 위해 벡터와 함께 사용하기로 한다. 자기곱을 구성하는 데 사용되는 주파수 fq는 신호 대역폭 내에 포함되어 있어야 하며, 또한 fLf1 ≤ … ≤ fQfU를 만족해야 함을 유의한다.

K개 도래각에 대해 gQ식 (4)와 같이 쓸 수 있다.

gQ=q=1Qk=1Kxkqakqθk,fq+gε(4) 

여기서 xkq = |xkq|exp(kq), ϕkq∈[0, 2π), gε는 잡음 항이다. 이중 하첨자 kqk번째 도래각과 q번째 주파수를 지칭한다. 식 (4)의 우변 첫 항은 신호성분으로서 총 KQ개 항의 합으로 구성된다. 그 중 K개 항은 자기항(self-terms)으로서 하나의 도래각만이 나타나는 항이며, 나머지 (KQ-K)개는 서로 다른 도래각으로 구성되는 교차항(cross-terms)이다.

신호성분에 대한 일반항 gsQ개 항의 Hadamard 곱으로 구성되는데, 다항전개에서 등장하는 선택지수(selection index) pkg를 도입하여 식 (5)와 같이 표현할 수 있다(23).

gs=x11a11p11x1Qa1Qp1Qx21a21p21x2Qa2Qp2QxK1aK1pK1xKQaKQpKQ=q=1Qk=1Kxkqakqθk,fqpkq(5) 

선택지수 pkg는 1 또는 0의 값을 가지며, 각각 xkqakq를 ‘선택’ 또는 ‘배제’함을 의미한다. pkq식 (6)과 같은 규칙을 따른다.

k=1Kq=1Qpkq=Q,k=1Kpkq=1(6) 

앞의 조건은 gsQ개 항으로 구성되도록 하며, 뒤의 조건은 주파수 fq에 대해 K개 도래각 가운데 하나의 도래각이 선택되도록 제한한다. 선택지수는 자기항과 교차항을 구분하고 추출하는 데 유용하다. 식 (7)과 같이 특정 k에 대해 합산한 값이 Q가 되는 경우, gsk번째 도래각에 대응하는 자기항이 되고, 그 합이 Q보다 작은 식 (8)과 같은 경우, gs는 교차항 중의 하나가 된다.

q=1Qpkq=Q(7) 
q=1Qpkq<Q(8) 

그에 대한 예시를 Table 1Table 2에 보였다. 각 열은 값이 1인 단 하나의 비영(non-zero) 요소만 포함하고 있는데, 식 (6)의 두 번째 조건을 만족하는 구성이다. Table 1은 비영 요소가 두 번째 행에만 존재하고 나머지는 모두 0으로 채워져 있다. 즉, 식 (7) q=1Qp2q=Q를 만족한다.

Selection index example for self-term extraction

Selection index example for cross-term extraction

따라서, Table 1에 해당하는 gs는 두 번째 (k = 2) 도래각 θ2에 대한 자기항이 된다. 반면, Table 2에서는 비영 요소의 위치가 열에 따라 다른데, 식 (8) 조건에 부합한다. 결과적으로, Table 2에 대한 gs는 교차항(cross-term) 중 하나가 된다.

주파수차 자기곱은 차주파수에서의 평면파 음장을 모사한다(14,16,19,20). 유사하게, 주파수합 자기곱 gQ는 합주파수 fΣ=q=1Qfq에서 평면파 음장을 모사한다. 이는 평면파 음장 yf,r=x expj2πfcesr을이용하여 증명할 수 있는데, x는 진폭, es는 음파진행 방향과 나란한 단위벡터, r은 위치벡터를 의미한다. 이때, 자기곱 음장은 식 (4)에서 K = 1, gε = 0으로 두고 얻을 수 있다.

식 (9)는 진폭 측면에서 차이가 있지만, 주파수가 f인 평면파와 그 형태가 동일하므로 자기곱 음장이 평면파 음장을 모사한다는 것을 알 수 있다.

gQK=1,gϵ=0=q=1Qxkqexpj2πfΣcesr(9) 

그러므로, 고차 자기곱 gQ와 합주파수 조향벡터 a(θ, f)간 빔형성인 gFS를 제안할 수 있다. 구체적으로, Q차 일반화 주파수 합 빔파워 BQ식 (10)과 같다.

BQ=|aHθ,fΣgQ2fΣ,L(10) 

상첨자 H는 허미시안 전치(Hermitian transpose)를, 기호 <ㆍ>f∑,LfL개 스냅샷에 대한 비상관 평균(incoherent averaging)을 의미한다.

Q = 1인 경우의 자기곱 g1은 단일 주파수 f1에서의 배열출력 y가 되므로, B1은 CBF와 등가하다. Q = 2인 경우의 gFS는 FS와 동일하며, 합주파수 f1+f2는 단일 주파수 (f1 또는 f2) 보다 더 높은 처리 주파수에 해당하므로 CBF 보다 좁은 빔폭을 얻을 수 있다. Q>2 조건에서는 빔폭은 더욱 좁아지며, gFS의 특징이 부각되기 시작한다. 높은 해상도를 위해 높은 Q가 바람직하지만, 지나치게 큰 Q는 공간 에일리어싱(spatial aliasing)에 의한 허위 도래각 발생을 초래한다. 최대 Q는 합주파수 f가 배열의 나이퀴스트 주파수 이하가 되도록 제한해야 한다(ffN). 파장 λ 관점에서 이 조건은 d/λ ≤ 0.5와 같으며, λ = c/f는 합주파수에 대응하는 파장이다. 한편, 잡음의 존재는 Q를 최대치까지 증가시키는 데 있어 추가 제약을 가하는데, 관련 논의는 3.1절에서 다룬다.

2.2 다중 도래각에 대한 제한

다중 도래각 상황 (K ≥ 2)에서는 (KQ-K)개의 허위 도래각이 존재하게 되어, gFS의 적용이 어렵다. 허위 도래각 θg는 교차항에 의한 것이며, 이는 공간 나이퀴스트 기준 위반에 의해 발생하는 허위 도래각과는 개념적으로 구분된다. 조향벡터 a(θ, f)와 신호 일반항 gs 간 내적(빔포밍)을 취한 결과는, 실 도래각 뿐만이 아니라 허위 도래각에서 국소 최대값을 가진다. 따라서, θg식 (11)과 같이 조향벡터 a(θ, f)와 벡터 ag(θk, fq)간 내적 최대화를 통해 계산할 수 있다. 여기서, ag(θk, fq)는 식 (5)에서 xkq = 1로 설정하고 얻은 교차항과 관련한 신호성분을 의미한다.

θg=argmaxθaHθ,fΣagθk,fq=argmaxθaHθ,fΣq=1Qk=1Kakqpkqθk,fq(11) 

식 (1)에서와 같이 조향벡터를 M으로 정규화함에 따라, a(θ, f)와 ag(θk, fq)의 l2-norm은 각각 1과 M1-Q가 되고, 두 벡터간 내적 최대값은 M1-Q가 된다. 따라서, 식 (11)식 (12)θg에 대해 푸는 것과 등가하다.

aHθg,fΣq=1Qk=1Kakqpkqθk,fq=M1-Q(12) 
ag(m)=1MQ/2q=1Qk=1Kexpj2πfq/cm-1d sinθkpkq=1MQ/2expq=1Qk=1Kj2πfq/cm-1d sinθkpkq(13) 

식 (12)식 (1)agm번째 성분에 대한 표현을 식 (13)에 적용하고, exp 함수의 위상항이 동일해야하는 조건을 반영하면, θg식 (14)와 같은 조건을 만족해야함을 알 수 있다.

fΣsinθg=k=1Kq=1Qfqpkqsinθk(14) 

또한, 정규화 주파수(normalized frequency) αk식 (15)와 같이 정의하면, 식 (14)식 (16)과 같이 간단히 쓸 수 있다.

ak=q=1QfqpkqfΣ(15) 
sinθg=k=1Kαksinθk(16) 

pkq의 특성상 0 ≤ αk ≤ 1이며, k=1Kαk=1인 점에 주목하면 식 (16)식 (17)과 같이 귀결되는데, 식 (17)과 같은 부등식은 허위 도래각 θg는 항상 실 도래각 범위 내에 존재함을 시사하고 있다.

sinθ1sinθgsinθKθ1θgθKθ1θ2θK(17) 

간단한 예시로, 두 개의 실 도래각 (θ1, θ2)에 대해 두 주파수 (f1, f2)로 구성된 2차 자기곱을 고려하자. 식 (16)에 따르면 두 개의 허위 도래각 θg1θg2식 (18a)식 (18b)와 같은 조건을 충족한다.

sinθg1=1-α1sinθ1+α1sinθ2sinθg2=α1sinθ1+1-α1 sinθ2(18a) 

또는

sinθg1=α2sinθ1+1-α2sinθ2sinθg2=1-α2sinθ1+α2sinθ2(18b) 

여기서 α1 = f1/(f1+f2), α2 = f2/(f1+f2)이고, 0 ≤ α1α2 ≤ 1을 만족한다. 만약 α1 = 0, (α2 = 1)이면, sinθg1(sinθg2)는 sinθ1(sinθ2)이 되고, 반대로 α1 = 1, (α2 = 0)일때, sinθg1(sinθg2)는 sinθ2(sinθ1)가 된다. 따라서 sinθ1 ≤ (sinθg1, sinθg2) ≤ sinθ2, 혹은 θ1 ≤ (θg1, θg2) ≤ θ2가 성립한다.

주파수 ‘차’ 기법에서는 허위 도래각이 실 도래각의 범위 내부 뿐 아니라 외부도 발생하므로, 식 (10)과 같은 주파수 평균을 통한 억제가 가능하다(15,21,22). 그러나, 주파수 ‘합’ 기법에서는 허위 도래각이 실 도래각 범위 [θ1, θK]에만 국한되므로, 평균화를 통한 허위 성분의 억제가 어렵게 되는 것을 증명할 수 있다.


3. 검 증

3.1 단일 도래각 시뮬레이션

가장 단순한 상황, 경계가 없는 공기 중 환경에서 단일 도래각 (K = 1)이 존재하는 경우를 고려하면, 센서 간격 d = 0.01 m, 센서 개수 M = 20개이다. 배열의 나이퀴스트 주파수는 fN = 17.0 kHz (c = 340 m/s) 이며, 음원의 대역폭은 [fl, fu] = [170 (= 0.01fN), 680 (= 0.04fN)] Hz로서, fN 대비 매우 낮은 주파수 대역을 가정한다.

이 시뮬레이션의 목적은 Q 증가에 따른 gFS의 분해능 향상을 확인하는 데 있다. 신호모델 식 (2)를 이용하여, L = 20개의 스냅샷을 생성하였고, 신호대 잡음비(signal-to-noise ratio, SNR)는 25 dB로 식 (19)와 같이 설정하였다.

SNR=10log10Ek=1Kaθkxk22En22(19) 

빔파워 산출시 빔폭을 일정하게 유지하기 위해 (해상도를 유지하기 위해), 합주파수 f 변화에 대한 평균은 취하지 않는다. 대신, 주어진 Q에 대해 신호대역 |fL, fU|를 Q개의 하위밴드(sub-band)로 균등분할하고, q(= 1, …, Q)번째 하위밴드의 중심 주파수를 fq로 설정하였다. 즉, 식 (20)의 분할조건에서 엘리어싱을 방지하지 하기 위한 최대 Q, 즉, 반파장 센서간격 조건 (d/λ = 0.5)을 만족하는 Q는 40이 된다.

fq=fL+fU-fLqQ+1(20) 

Fig. 1은 실 도래각 θ1 = 0°에 대해, Q 증가에 따른 gFS의 특성 변화를 관찰하고 있다. Q = 1일 때, 즉, Fig. 1(a)의 CBF는 실 도래각 주변으로 평탄하고 넓은 빔을 가지는데, Rayleigh 분해능 한계 ∆θ = sin-1(λ/Md)가 정의되지 않는다는 점 역시 뚜렷하지 않는 주엽의 특성을 반영한다(3). Fig. 2(b)Q = 2(FS)인 경우에도, 그림에 표시된 d/λ값이 작으므로 해상도 향상은 제한적이다.

Fig. 1

Single DOA simulation result: BQ. Dotted lines represent the true DOA (θ1 = 0°)

Q = 10일때, 해석 주파수 증가에 의해 주엽이 수축되고, 부엽이 나타나기 시작한다. Q가 최대차수 40에 도달하면, Fig. 1(d)에서 보듯이 예리한 빔이 형성되고 만족스러운 해상도를 얻을 수 있다. Q가 최대값을 초과하면, 주엽은 더욱 수축하지만 Fig. 1(e)과 같이 엘리어싱에 의한 허위 도래각이 발생한다.

0° 이외 도래각에 대한 검증은 식 (21)의 Gram 행렬 G를 이용한다(11).

G=AHA(21) 

행렬 A=|a(θ1, f) … a(θN, f)|∈ℂM×N는 합주파수 f에서 산출된 모든 조향벡터를 수집하여 구성되므로, G의 각 열(또는 행)은 도래각 θ1∈[-90°, 90°]에 대한 빔패턴(beampattern)으로 간주할 수 있다.

Fig. 1에 사용된 Q에 대해 산출한 Gram 행렬을 Fig. 2에 나타냈다. 배열에 대한 수직방향(broadside, 그림에서 0°에 가까운 영역) 입사조건에서는 빔폭이 좁게 유지되지만, 측면입사(endfire, 그림에서 ±90°에 가까운 영역)에 가까울수록 빔폭이 넓어진다. 중요한 점은, Q 증가에 따라 주엽이 더 예리하고 집중되는 형태로 변하며, 이는 곧 해상도 향상을 의미한다는 것이다. Q가 최대차수에 도달하면 ±90° 부근에서만 엘리어싱이 발생한다(Fig. 2(d)). 최대차수를 크게 웃돌 경우, 실 도래각 외 허위 도래각 성분이 지배적으로 발생한다(Fig. 2(e)).

Fig. 2

Gram matrix G

특정 도래각 θ1 = -45°, 0°, +45°에 대해 Q를 점진적으로 증가시키며 얻은 빔파워를 Fig. 3에 보였다. 에일리어싱이 유발하는 허위 도래각 발생을 방지하기 위해서는 d/λ < 0.5 조건을 만족하도록 Q가 지정되어야 함을 다시 한번 강조하고 있다.

Fig. 3

2D map of beampower BQ as a function of the order Q and look angle θ for the true DOAs. The value in parentheses on the horizontal-axis is the d/λ∑ corresponding to each Q. Dotted line shows the maximum Q, below which is aliasing free

Fig. 4는 SNR을 5 dB 수준으로 낮추어 얻은 결과이며, Fig. 1과 대비하여 볼 때 Q가 증가할수록 잡음에 민감하게 반응하고 있음을 알 수 있다. 해석적 접근을 위해 K = 1(단일 도래각)인 경우의 gQ식 (2)식 (3)을 이용하여 전개하면, 식 (22)이며, 여기서 xqaqx1qa1q(θ1, fq)의 약식 표기이다.

Fig. 4

Single DOA simulation result for a low SNR of 5 dB (L = 20): BQ. Dotted lines represent the true DOA (θ1 = 0°)

gQK=1=q=1Qxqaq+nq(22) 

첨자 k = 1을 생략하여 선택지수를 단순화하고, pq = 1 또는 0으로 두면, 식 (22)식 (23)과 같이 쓸 수 있다.

gQK=1=q=1Qxqaqpqnq1-pq(23) 

pq = 1을 만족할 경우 q=1Qpq=Q,gqK=1의 신호항이 추출되며, q=1QpqQ인 경우에는 잡음 항이 분리된다. 따라서 gq|k=1식 (24)와 같은 신호-잡음 분해가 가능하다.

gQK=1=gs+gε=q=1Qxqaq+q=1Qxqaqpqnq1-pq(24) 

잡음항 gε는 총 (2Q-1)개의 항으로 구성되며, 각 항은 신호항 gs와 동일한 차수 Q를 갖는다. Q가 증가함에 따라 gε에 포함되는 항의 개수가 급격히 증가하므로, 전체 자기곱에 기여하는 영향도 그에 비례해 커지게 됨을 알 수 있다.

여러 SNR 및 Q에 대해 계산된 빔파워 BQ의 최대부엽수준(MSL: maximum sidelobe level)을 Fig. 5에 보였다. 높은 SNR 조건에서는 Q에 관계없이 낮은 수준을 유지하지만, 잡음이 큰 경우(저 SNR조건) Q가 증가할수록 MSL이 상승하게 된다. 잡음이 존재하는 실제 환경에서 Q를 최대치까지 증가시키는 것이 제한될 수 있음을 의미하나, 적절한 Q를 선택한다면 이 기법은 저 SNR 환경에서도 실용 가능한 수준의 해상도를 제공할 수 있음을 시사한다.

Fig. 5

Maximum sidelobe level (MSL) of BQ as a function of Q and SNR. Results for Q < 5 are not reported as no sidelobe is present in that region

3.2 다중 도래각 시뮬레이션

다중 도래각 시뮬레이션의 목적은 교차항이 유발하는 허위 도래각 성분의 부정적 영향을 확인하기 위함이다. 배열 구성 및 환경은 단일 도래각 상황과 동일하지만, 교차항 유발 허위 도래각을 시각화하기 위해 주파수 대역폭을 |fL, fU| = |0.01fN, 0.6fN|으로 확장하였다.

K = 2, K = 3의 두 가지 경우를 고려하며, 전자의 경우 실 도래각은 θ1 = -45, θ2 = +45에, 후자의 경우 θ1 = -45, θ2 = 0, θ3 = +45에 위치한다. 음원을 xkq = exp(kq)로 모형화하고, 위상 ϕkq는 균등분포(uniform distribution)를 따른다고 가정하였다, ϕkq ~ U|0, 2π). 즉, 단위 진폭을 가진 음원은 서로 비상관적(incoherent)임을 가정한다.

허위 도래각이 주파수에 대해 변하도록 강제하고, 동시에 이를 최대한 억제하기 위해 다음과 같은 주파수 평균화를 적용한다. 주어진 대역 |fL, fU|의 양 끝에서 오프셋(offset) fr만큼 차감하여 얻은 i번째 하위밴드 |fli, fui|에 대해, 식 (20)의 균등 분할방식을 적용한다. 부연하면, i(= 1, …, F)번째 하위밴드는 다시 Q개의 하위-하위밴드(sub-sub-bands)로 균등 분할되며, 이로부터 i-q번째 하위-하위밴드의 주파수 fiq가 결정된다. 수식으로 표현하자면 식 (25), 식 (26)과 같다.

fiq=fli+fli-fliqQ+1,i=1,F,q=1,,Q(25) 
fli=fL+i-1frfui=fU-i-1fr(26) 

여기서 F = (fU-fL)/2fr은 하위밴드의 개수, 즉 주파수 평균횟수이다. 이 시뮬레이션에서는 fr을 10으로 설정하였으며, 이에 F = 511개의 하위밴드가 생성된다. 각 i번째 하위밴드에 대해 스냅샷 L = 20개를 생성하므로, 전체 L×F = 10 200회 평균을 통해 빔파워를 얻는다. 이보다 큰 회수의 평균을 취하더라도, 식 (26)에 보이는 결과보다 더 이상의 주목할만한 변화가 감지되지 않았다.

Fig. 6K = 2(좌측열)와 K = 3(우측열)에 대한 결과를 열로, 세가지 Q에 대한 결과를 행으로 구분하여 보이고 있다. 상단, 중간, 하단 행은 각각 Q = 1(CBF), Q = 2(FS), Q = 3(gFS)에 해당한다. Q ≥ 4인 경우, 합주파수 f가 배열의 나이퀴스트 주파수를 초과하여 엘리어싱에 의한 허위 도래각 성분이 발생하므로 고려하지 않는다.

Fig. 6

Multi-DOAs simulation for K = 2(θ1 = -45°, θ2 = +45° ) (left), K = 3(θ1 =-45°, θ2 = 0°, θ3 =+45° ) (right), and with Q = 1 (top), Q = 2 (mid), Q = 3 (bottom row). The left subfigure in each panel shows |aHgQ|2 for each ith(i = 1, ⋯ F) band, and the right subfigure shows their average along the horizontal axis, i.e., BQ. Dotted black lines indicate the true DOAs, while dashed white lines denote the ghost DOAs computed using Eq. (16)

각 그림은 두 개의 하위 그림으로 구성되며, 좌측은 i번째 하위밴드의 결과를 2D맵으로, 우측은 그 수평축에 대한 평균값인 BQ를 보여준다. 단, 2D맵의 수평축은 하위밴드 인덱스 i 대신 정규화 주파수 α1(= fi1/f)으로 표시하였다.

CBF에서는 교차항이 발생하지 않으므로, 실 도래각 주변으로 정적 성분만이 형성됨을 볼 수 있으며, 빔폭은 넓지만 실 도래각 성분을 식별할 수 있다(상단행). 반면, FS 해석 결과에서는 실 도래각 범위 내에서 동적으로 변하는 허위 도래각이 나타나며(중간행), 식 (16)의 계산결과와 일치함을 흰색 점선을 통해 확인할 수 있다. 허위 도래각은 주파수 평균 이후에도 잔류하며, 예기치 않은 빔파워 증가를 유발한다. 이 예시에 한해 실 도래각을 식별할 수 있으나, 식 (25)외 다른 분할 방식이 사용될 경우 그 구분이 더욱 모호할 수 있다. gFS에서는 허위 성분의 기여도가 지배적이다(하단행). 평균화를 수행하더라도 이 성분들은 θ1θK 사이에 집중되어 여러 국소 피크를 형성한다. 이로 인해, 실제보다 더 많은 도래각이 존재하는 것처럼 보이게 만드는 오해를 초래한다.

Fig. 7은 단일 도래각 시뮬레이션과 같은 저주파 대역 |fL, fU| = |0.01fN, 0.04fN|에 대한 해석결과를 보이고 있다. 고 SNR 조건에서는(상단), Q 증가에 따라 실 도래각 범위내 파워 집중 현상이 더욱 뚜렷해지고, 저 SNR 조건에서는(하단) 그 성능이 더욱 악화됨을 알 수 있다.

Fig. 7

Beampower for multi-DOAs with a low frequency bandwidth [fL, fU] = [0.01fN, 0.04fN]. Vertical dotted lines represent the true DOAs. All cases are with L = 20

3.3 단일 도래각 실험

단일 도래각에 상황에 대한 유효성 검증을 위해, Fig. 8과 같이 관련 실험을 무향실에서 수행하였다. 해당 무향실은 쐐기 끝단 사이의 내부 유효공간이 약 34.7 m3(폭 4.24 m × 깊이 3.67 m × 높이 2.23 m)이며, 250 Hz 이상의 주파수에서 자유음장 조건을 제공한다.

Fig. 8

Experimental setup in an anechoic chamber

저주파 핑크 노이즈를 송출하는 스피커는 배열과 동일한 수평면상에서 배열 중심으로부터 3 m 이격된 위치에 설치되었다(Fig. 8(a)). 참고로, 실험에 사용한 센서배열은 시뮬레이션과 동일하다(d = 0.01 m, M = 20, microphone: SISPIA HJ06)(Fig. 8(b)). 핑크 노이즈 신호 대역은 20 Hz ~ 1 kHz이며, 배열 중심 위치에서의 음압레벨은 약 75 dB SPL로 측정되었다.

이상과 같은 조건에서 다음 두 가지 경우에 대해 실험을 수행하였다. 첫 번째는 스피커가 배열을 정면으로 향하도록 배치한 경우(θ1≃0°), 두 번째는 스피커가 θ1≃-40° 방향을 향하도록 배치한 경우이다. 20 kHz의 샘플링 주파수(신호수집장치: NI-9234)로 마이크로폰 신호를 취득하였고, 스냅샷 당 2048 샘플수로 FFT 변환하였다. 후처리를 위한 나머지 설정은 3.1절에서 기술한 내용과 동일하다. 즉, 주파수 대역폭은, |fL, fU| = |170Hz(= 0.01fN), 680Hz(= 0.04fN)|스냅샷 수는 L = 20이다.

Fig. 9(a)는 정면 방향(θ1≃0°)에서의 빔파워를 보이고 있다. Q = 1인 CBF의 경우, 평탄한 빔 형태로 인해(낮은 해상도로 인해) 도래각 파악이 어렵다. 3.1절에서 설명한 바와 같이, gFS 해석에서는 잡음의 영향으로 인해 Q를 최대차수인 40까지 증가시키는 것은 적절하지 않다. 평면파 가정과 실제 파형간의 불일치, 무향실 내 음파반사 등 다양한 요인이 잡음에 기여하는 것으로 판단된다. 몇 회의 시행착오를 통해 낮은 부엽 수준과 함께 원하는 도래각 성분을 식별할 수 있는 Q값을 7로 선정하였다. 실제 적용에서는 Q를 점진적으로 증가시키면서, 부엽이 뚜렷하게 나타나기 직전의 값을 선택하는 것이 gFS에 대한 실용적 접근이라 할 수 있다. 하단에 θ1≃-40°의 결과를 도시하였고, 측면입사에 가까운 도래각 조건으로 인해 빔폭이 약간 두꺼워지는 점을 제외하면 상단의 결과와 유사한 특성을 보인다.

Fig. 9

Results of the single DOA experiment


4. 결 론

신호파장이 배열크기에 비해 매우 큰 저주파 신호에 대해 고해상도 빔포밍을 구현하는 것은 대단히 어려운 문제이다. 이를 위해, 이 논문에서는 Q차 주파수합 자기곱 gQ를 활용한 일반화된 주파수합 빔포밍(gFS) 기법을 제안하였다.

Q값은 합주파수(f)가 배열의 공간 나이퀴스트 주파수 이하가 되도록 설정되어야 하나, 잡음의 존재는 Q의 이론적 최대값까지 증가시키는 데 제약요인이 된다. 또한, 다중 도래각 상황에서 교차항에 의한 허위 도래각은 실 도래각 범위 내에 존재하므로 단순한 주파수 평균만으로는 억제할 수 없음을 해석적으로 설명하였다. 이처럼 고잡음/다중 도래각 환경에 대한 제약이 따르지만, 제안된 기법은 단일 저주파 음원에 대한 빔포밍 해상도를 높일 수 있는 유의미한 대안으로 평가된다.

Acknowledgments

이 논문은 2025년 ~ 2026년도 국립창원대학교 자율연구과제 연구비 지원으로 수행된 연구결과임.

References

  • Gerstoft, P., Mecklenbräuker, C. F., Seong, W. and Bianco, M., 2018, Introduction to Compressive Sensing in Acoustics, Journal of the Acoustical Society of America, Vol. 143, No. 6, pp. 3731~3736. [https://doi.org/10.1121/1.5043089]
  • Krim, H. and Viberg, M., 1996, Two Decades of Array Signal Processing Research: The Parametric Approach, IEEE Signal Processing Magazine, Vol. 13, No. 4, pp. 67~94. [https://doi.org/10.1109/79.526899]
  • Van Trees, H., 2002, Arrays and Spatial Filters, in Optimum Array Processing : Part IV of Detection, Estimation, and Modulation Theory, Hoboken, John Wiley & Sons, Chapter 2, NJ, United States, pp. 17~89. [https://doi.org/10.1002/0471221104.ch2]
  • Jeon, J. Y., Lee, P. J. and Sato, S.-I., 2009, Use of the Standard Rubber Ball as an Impact Source with Heavyweight Concrete Floors, Journal of the Acoustical Society of America, Vol. 126, No. 1, pp. 167~178. [https://doi.org/10.1121/1.3148193]
  • Širović, A., Hildebrand, J. A. and Wiggins, S. M., 2007, Blue and Fin Whale Call Source Levels and Propagation Range in the Southern Ocean, Journal of the Acoustical Society of America, Vol. 122, No. 2, pp. 1208~1215. [https://doi.org/10.1121/1.2749452]
  • Lee, J.-H., Park, H.-G., Kim, J.-H., Lee, K.-J. and Seo, J.-S., 2014, Reduction of Propeller Cavitation Induced Hull Exciting Pressure by a Reflected Wave from Air-bubble Layer, Ocean Engineering, Vol. 77, pp. 23~32. [https://doi.org/10.1016/j.oceaneng.2013.12.007]
  • Lee, J. and Lee, K., 2018, Balloon for Reducing Propeller Caviation Induced Hull Excitation, Transactions of the Korean Society for Noise and Vibration Engineering, Vol. 28, No. 1, pp. 39~49. [https://doi.org/10.5050/KSNVE.2018.28.1.039]
  • Lee, J.-H. and Lee, K.-J., 2019, Multi-balloons Design for Suppressing Propeller BPF Induced Hull-excitation at Multi-frequencies, Transactions of the Korean Society for Noise and Vibration Engineering, Vol. 29, No. 5, pp. 592~599. [https://doi.org/10.5050/KSNVE.2019.29.5.592]
  • Schmidt, R., 1986, Multiple Emitter Location and Signal Parameter Estimation, IEEE Transactions on Antennas and Propagation, Vol. 34, No. 3, pp. 276~280. [https://doi.org/10.1109/TAP.1986.1143830]
  • Capon, J., 1969, High-resolution Frequency-wavenumber Spectrum Analysis, Proceedings of the IEEE, Vol. 57, No. 8, pp. 1408~1418. [https://doi.org/10.1109/PROC.1969.7278]
  • Xenaki, A., Gerstoft, P. and Mosegaard, K., 2014, Compressive Beamforming, Journal of the Acoustical Society of America, Vol. 136, No. 1, pp. 260~271. [https://doi.org/10.1121/1.4883360]
  • Gerstoft, P., Xenaki, A. and Mecklenbräuker, C. F., 2015, Multiple and Single Snapshot Compressive Beamforming, Journal of the Acoustical Society of America, Vol. 138, No. 4, pp. 2003~2014. [https://doi.org/10.1121/1.4929941]
  • Abadi, S. H., Haworth, K. J., Mercado-Shekhar, K. P. and Dowling, D. R., 2018, Frequency-sum Beamforming for Passive Cavitation Imaging, Journal of the Acoustical Society of America, Vol. 144, No. 1, pp. 198~209. [https://doi.org/10.1121/1.5045328]
  • Abadi, S. H., Song, H. C. and Dowling, D. R., 2012, Broadband Sparse-array Blind Deconvolution Using Frequency-difference Beamforming, Journal of the Acoustical Society of America, Vol. 132, No. 5, pp. 3018~3029. [https://doi.org/10.1121/1.4756920]
  • Douglass, A. S., Song, H. C. and Dowling, D. R., 2017, Performance Comparisons of Frequency-difference and Conventional Beamforming, Journal of the Acoustical Society of America, Vol. 142, No. 3, pp. 1663~1673. [https://doi.org/10.1121/1.5003787]
  • Douglass, A. S. and Dowling, D. R., 2019, Frequency-difference Beamforming in the Presence of Strong Random Scattering, Journal of the Acoustical Society of America, Vol. 146, No. 1, pp. 122~134. [https://doi.org/10.1121/1.5114811]
  • Xie, L., Sun, C. and Tian, J., 2020, Deconvolved Frequency-difference Beamforming for a Linear Array, Journal of the Acoustical Society of America, Vol. 148, No. 6, EL440-EL446. [https://doi.org/10.1121/10.0002927]
  • Worthmann, B. M., Song, H. C. and Dowling, D. R., 2015, High Frequency Source Localization in a Shallow Ocean Sound Channel Using Frequency Difference Matched Field Processing, Journal of the Acoustical Society of America, Vol. 138, No. 6, pp. 3549~3562. [https://doi.org/10.1121/1.4936856]
  • Geroski, D. J. and Dowling, D. R., 2019, Long-range Frequency-difference Source Localization in the Philippine Sea, Journal of the Acoustical Society of America, Vol. 146, No. 6, pp. 4727~4739. [https://doi.org/10.1121/1.5138124]
  • Worthmann, B. M. and Dowling, D. R., 2017, The Frequency-difference and Frequency-sum Acoustic-field Autoproducts, Journal of the Acoustical Society of America, Vol. 141, No. 6, pp. 4579~4590. [https://doi.org/10.1121/1.4985440]
  • Lee, J.-H., Park, Y. and Gerstoft, P., 2024, Direction of Arrival Estimation Using Compressive Frequency Difference Beamforming, Transactions of the Korean Society for Noise and Vibration Engineering, Vol. 34, No. 4, pp. 402~410. [https://doi.org/10.5050/KSNVE.2024.34.4.402]
  • Park, Y., Gerstoft, P. and Lee, J., 2022, Difference-frequency MUSIC for DOAs, IEEE Signal Processing Letters, Vol. 29, pp. 2612~2616. [https://doi.org/10.1109/LSP.2022.3230365]
  • Stanley, R. P., 1997, Enumerative Combinatorics, 2nd Edition, Cambridge University Press, Cambridge, United Kingdom.

Fig. 1

Fig. 1
Single DOA simulation result: BQ. Dotted lines represent the true DOA (θ1 = 0°)

Fig. 2

Fig. 2
Gram matrix G

Fig. 3

Fig. 3
2D map of beampower BQ as a function of the order Q and look angle θ for the true DOAs. The value in parentheses on the horizontal-axis is the d/λ∑ corresponding to each Q. Dotted line shows the maximum Q, below which is aliasing free

Fig. 4

Fig. 4
Single DOA simulation result for a low SNR of 5 dB (L = 20): BQ. Dotted lines represent the true DOA (θ1 = 0°)

Fig. 5

Fig. 5
Maximum sidelobe level (MSL) of BQ as a function of Q and SNR. Results for Q < 5 are not reported as no sidelobe is present in that region

Fig. 6

Fig. 6
Multi-DOAs simulation for K = 2(θ1 = -45°, θ2 = +45° ) (left), K = 3(θ1 =-45°, θ2 = 0°, θ3 =+45° ) (right), and with Q = 1 (top), Q = 2 (mid), Q = 3 (bottom row). The left subfigure in each panel shows |aHgQ|2 for each ith(i = 1, ⋯ F) band, and the right subfigure shows their average along the horizontal axis, i.e., BQ. Dotted black lines indicate the true DOAs, while dashed white lines denote the ghost DOAs computed using Eq. (16)

Fig. 7

Fig. 7
Beampower for multi-DOAs with a low frequency bandwidth [fL, fU] = [0.01fN, 0.04fN]. Vertical dotted lines represent the true DOAs. All cases are with L = 20

Fig. 8

Fig. 8
Experimental setup in an anechoic chamber

Fig. 9

Fig. 9
Results of the single DOA experiment

Table 1

Selection index example for self-term extraction

q
k
1 2 3 Q q=1Qpkq
1 0 0 0 0 0
2 1 1 1 1 Q
3 0 0 0 0 0
: : : : : : :
K 0 0 0 0 0
k=1Kpkq 1 1 1 1 k=1Kq=1Qpkq=Q

Table 2

Selection index example for cross-term extraction

q
k
1 2 3 Q q=1Qpkq
1 0 1 0 0 1 (< Q)
2 1 0 1 0 2 (< Q)
3 0 0 0 0 0 (< Q)
: : : : : : :
K 0 0 0 1 1 (< Q)
k=1Kpkq 1 1 1 1 k=1Kq=1Qpkq=Q