I. 서 론
잡음 상황에서 정확한 정현파 파라미터(주파수)를 추정하는 문제는 다양한 분야에서 기본적으로 요구되는 신호처리 기술이다. 특히 수중 음향신호를 수신하여 특징을 분석하는 소나체계에서는 정확한 주파수 분석이 선행되어야 식별 성능을 보장 받을 수 있다. 레이더 신호처리에서도 빠르게 이동하는 표적을 정확하게 추적하기 위해서 짧은 시간 단위로 갱신되는 도플러 주파수 정보를 정확하게 추정할 수 있어야 한다.
일반적으로 단일 주파수 추정은 DFT(Discrete Fourier Transform)에 의해 변환된 신호의 최대 크기 성분에 해당되는 주파수 빈(Bin)을 찾는 방법으로 해결한다. 하지만, 정확한 주파수 추정을 위해서는 관련된 다양한 요소들의 상호관계를 고려 할 필요가 있다. 신호의 DFT 변환은 주파수 영역에서 정수값을 갖는 빈에 의해 이산화 된다. 따라서 실제 최대값이 주파수 빈 사이에 존재한다면 오차(𝛿)가 발생하게 된다. 이러한 문제를 해결하기 위하여 DFT에 사용되는 데이터 샘플을 증가하여 주파수 해상도를 높일 수도 있으나 연산량이 증가함과 동시에 시간 - 주파수 불확실성(Heisenberg - Gabor uncertainty)에 의해서 시간 영역에서의 해상도가 낮아지게 된다. 한편 주파수 분석주기 불일치에 의한 주파수 누설 현상과 이에 따른 추정 성능 저하를 완화하고자 다양한 윈도우 함수가 제안되었으며 신호 분석 목적 및 신호 환경에 따라 윈도우 함수를 선택하여 사용하고 있다.
주파수 추정에 관계된 요소들이 성능에 미치는 영향력 분석과는 별개로 추정 알고리즘에 대한 연구도 활발하게 진행되고 있다. 관련된 연구 중 가장 대표적인 방법은 DFT 최대 크기를 가지는 주파수 빈과 인근 빈 2개 ~ 3개의 결과값을 이용하여 보간에 의해 실제 주파수를 추정하는 방법이다.[1-4]
또 다른 방법으로 신호의 위상 정보를 주파수 추정에 활용하는 방법이다.[5] 이 방법은 신호의 1차원(주파수 분석) 표현뿐만 아니라, 2차원(시간 – 주파수 영역) 표현(TFR, Time - Frequency Representation)까지 적용될 수 있다. 즉 위상 및 그룹지연(group delay) 정보를 이용하여 시간과 주파수 측정의 정확성을 높이는 재할당 기술로 다양한 분야에서 활용되고 있다.[6,7]
한편, 최근에 신호처리 기술 관련 인터넷 사이트[8]를 통하여 잡음이 없는 상황에서 DFT 값을 이용한 정확한 주파수 계산 해석해(closed form solution) 풀이 방법이 Cedron 에 의해 제안되었다. 비록 잡음이 없는 조건에서의 결과이지만, 추정이 아닌 정확한 해석해를 도출 한다는 특징으로 인하여 관련 분야 연구자들에게 많은 관심을 받고 있다.[11] 그럼에도 불구하고 아직까지 엄밀한 검증을 받지 않은 실험 결과라 이에 대해 전문가 집단의 검증과 확인이 필요한 상황이다.
본 논문에서는 정확한 주파수 추정 알고리즘의 공학문제 적용/응용을 위한 예비연구 목적으로 단일 주파수 신호에 대한 추정 알고리즘들의 성능을 다양한 실험 환경에서 비교, 분석 하였다. 총 4개의 주파수 추정 알고리즘 성능을 DFT에 의한 결과와 비교하였다. 특히 Cedron 알고리즘은 이러한 분석과정을 통하여 관련 파라미터[SNR(Signal to Noise Ratio), 윈도우 함수 종류, 윈도우 길이 등] 가 성능에 미치는 영향을 확인하고 향후 다양한 분야에 활용 될 수 있는 초석을 마련 하고자 하였다.
본 논문의 2장에서는 신호 모델을 정의하고 주파수 추정에 영향을 주는 요소들을 정리하였다. 3장에서는 비교 대상 알고리즘을 소개하고, 4장에서는 성능지표를 소개하였다. 5장에서는 모의실험을 위한 실험환경 및 실험결과를 제시하였으며, 마지막 6장에서 연구 결과를 종합하였다.
II. 신호 모델 및 주파수 추정시 영향요소
본 논문에서 분석하고자 하는 신호는 아래와 같은 모델의 단일 주파수 실수 정현파 이산 신호이다.
| $$s_n=M\;\cos\;(2\pi\;f_{norm}\;n+\varnothing)+w_n$$ | (1) |
여기서 M는 진폭, fnorm은 0에서 0.5까지의 값을 가질 수 있는 정규화 주파수(normalized frequency), 𝜙는 위상이며, wn은 WGN(White Gaussian Noise) 이다. 신호의 총 길이(샘플 수) 는 N이며, n은 0부터 N-1까지의 값을 가지는 시간 인덱스이다.
본 논문에서 논의되는 알고리즘의 목표는 샘플값 sn을 통하여 정확한 fnorm을 찾는 것이다. 일반적으로 Eqs. (2) ~ (4)의 일련의 연산을 통하여 DFT 최대 크기 주파수 빈에 해당하는 주파수를 추정()한다.
| $$Z_k=\sum_{n=0}^{N-1}\;x_ne^{-i2\pi\frac kNn}.$$ | (2) |
| $$k\ast\;=\mathrm{argmax}\vert Z_k\vert\;\;\;\mathrm{for}\;\;\;\forall\;k.$$ | (3) |
| $${\widehat f}_{norm}=k\ast/N.$$ | (4) |
하지만 이러한 방법은 오차(𝛿)를 포함하고 있는 개략적인 주파수를 찾을 때에는 효과적이지만 정확한 주파수를 찾고자 할 때에는 한계를 가진다. 이는 N 값에 의해 결정되는 주파수 해상도와 관련된다. 즉, 정확한 주파수가 주파수 빈 과 빈 사이에 존재할 수 있다는 사실 때문이다. 이 경우 실제 주파수와 N의 부조화(주파수 해상도 비정수배)로 인해 주파수 누설 현상이 발생하여 주파수 측정시 오차가 발생한다(Fig. 1).
물론 분석에 사용되는 데이터양을 증가하여 주파수 해상도를 높이고 추정 오차를 줄일 수 있지만 이에 따라 연산량이 증대된다. 더불어 분석 데이터를 프레임 단위로 나누어 시간과 주파수 영역에서 동시에 분석 결과를 표현하고자 할 때에는 Eq. (5)에 따른 불확실성 제한 조건(Heisenberg - Gabor uncertainty) 으로 시간 영역의 해상도가 낮아진다.
| $$\sigma_t\;\cdot\;\sigma_f\geq\frac1{4\pi},$$ | (5) |
여기서 𝜎t와 𝜎f는 각각 시간과 주파수 추정값의 표준편차이다.
분석 신호 주기 불일치에서 따른 주파수 누설은 강한 부엽(side lobe)에 묻힌 약한 신호를 못 찾게 하거나(amplitude resolution issue), 인접한 두 개의 주파수를 구분하지 못하는 영향(frequency resolution issue) 을 초래할 수 있다. 이러한 현상은 윈도우 함수에 의해 일부 완화될 수 있다. 윈도우 함수는 분석하고자 하는 신호의 특성 및 목적에 맞추어 다양한 함수 중 선택을 할 수 있는데, 본 연구에서는 Rectangular, Hamming, Hanning, Blackman 4가지 종류의 윈도우 함수 영향을 분석하였다. 크기가 다른 인접한 주파수를 가진 테스트 신호, cos(2π303t)+0.3cos(2π315t), 를 SNR = 3 dB 환경에서 합성하고 다양한 윈도우 함수에 의한 DFT 분석 결과를 Fig. 2에 나타내었다. Rectangular 윈도우는 주파수 분해능은 뛰어난 반면 강한 부엽 발생과 정확한 주파수 성분 크기 추정이 안 되는 사실을 확인 할 수 있다. 반면 그 외 윈도우 함수는 다소간의 편차를 가지지만, 주파수 성분 크기 추정이 우수한 사실을 확인 할 수 있다(크기 1을 가진 303 Hz 성분은 -3.01 dB, 0.3 크기를 가진 315 Hz 성분은 -13.47 dB을 나타내야 한다).
III. 주파수 추정 알고리즘
간략히 살펴본 바와 같이 정확한 주파수 추정을 위해서는 고려되어야 할 요소들이 많다. 이러한 요소들의 영향력 분석과는 별도로 정밀한 주파수 추정을 위해 다양한 알고리즘들이 제시되었으며 본 장에서는 다음의 네 가지 알고리즘(Jacobsen, Candan, Reassignment and Cedron)을 소개하고자 한다.
첫 번째와 두 번째 방법은 최대 DFT 결과값을 가지는 주파수 빈과 인근 빈 2개의 DFT 결과값을 이용하여 정밀한 주파수를 추정하는 방법이다. 앞서 Eqs. (2) ~ (4)에서 살펴본 바와 같이 실제 최대값을 가지는 주파수 빈의 위치는 정수 k*와 𝛿만큼의 오차를 가지고 있다. 즉,
| $$k_{peak}=k\ast+\delta.$$ | (6) |
Jacobsen에 의해 제안된 첫 번째 방법은 다음과 같이 δ를 계산하여 주파수를 추정한다.[1]
| $${\widehat\delta}_J=-Re\;\left[\frac{\left(Z_{k^\ast+1}-Z_{k^\ast-1}\right)}{\left(2Z_{k^\ast}-Z_{k^\ast-1}-Z_{k^\ast+1}\right)}\right].$$ | (7) |
Candan은 Jacobsen이 제안한 방법을 수정하여 잔여 bias 제거가 가능한 두 번째 방법을 제안하였다.[3] 즉,
| $${\widehat\delta}_C=\frac{\tan^{-1}\left({\widehat\delta}_C\cdot\tan\;\left(\pi/N\right)\right)}{\pi/N}.$$ | (8) |
세 번째로 주파수 재할당 신호 처리 방법이다. 이 방법은 1차원 음향 신호의 2차원 표현 방법인 스펙트로그램에 적용하여 시간-주파수 격자에서 데이터 위치의 정확성을 높여 주는데 활용되었다.[12]
| $$\overline S_x^h(\tau,w)=\iint_{-\infty}^\infty S_x^h\left(s,\eta\right)\;\delta\left(\tau-{\widehat\tau}_x(s,\eta),w-{\widehat w}_x(s,\eta)\right)\frac1{2\pi}dsd\eta.$$ | (9) |
주어진 신호 x에 대한 주파수 재할당 스펙트로그램 은 Eq. (9)와 같이 나타낼 수 있다. 여기서 은 h 윈도우 함수를 이용한 스펙트로그램을 의미하며, δ(τ,w)은 2차원 Dirac impulse 함수이다. 시간-주파수 격자([s, η])에 위치한 스펙트로그램 데이터는 Dirac impulse 함수와의 적분에 의하여 및 만큼 위치가 이동된다. 두 수치는 Eqs. (10)와 (11)에 의해서 구한다.
| $${\widehat\tau}_x(s,\eta)=s+Re\left(F_x^{th}(s,\eta)/F_x^h(s,\eta)\right).$$ | (10) |
| $${\widehat w}_x(s,\eta)=\eta+Im\left(F_x^{dh}(s,\eta)/F_x^h(s,\eta)\right).$$ | (11) |
여기서 Fhx, Fthx와 Fdhx은 각각 h(t), t.h(t), dh(t)/dt 윈도우 함수를 이용한 STFT(Short-time Fourier Transform) 결과이다. 이러한 재할당 과정을 통하여 각 주파수 빈의 기하학적 중심에서 에너지 중심으로 데이터 측정값이 이동하며 이를 통하여 주파수 측정의 정확도를 향상시킬 수 있다. 본 논문에서는 N개 데이터에 대한 재할당 신호 처리 후 각 프레임에 포함된 시간 정보(τ)는 무시하였다. 재할당 주파수 정보 중 가장 큰 크기를 가지는 주파수를 으로 선정한다.
마지막 방법은 Cedron으로 k* 및 k* ± 1에서의 DFT 값을 이용한 대수조작을 통하여 정확한 주파수 도출이 가능한 해석해를 도출한다.[9] 유도된 최종 식은 Eq. (12)와 같다.
| $${\widehat f}_{norm}=\frac1{2\pi}\cos^{-1}\left(\cos\;\alpha\right)\;,$$ | (12) |
여기서 βk=k.2π/N, Rk = e-jβk이다. 자세한 수식 유도 과정은 부록에 소개하였다.
Eq. (12)의 결과는 Zk를 실수부(Re(Zk))와 허수부(Im (Zk) = yk)로 나누고 선형대수 기법을 이용하여 보다 일반화 할 수 있다.[10] 더욱이 수식 유도과정 중 잡음을 고려하여 앞서 제안한 방법보다 좋은 성능을 가질 수 있도록 하였다.
| $${\widehat f}_{norm}=\frac1{2\pi}\cos^{-1}\left(\frac{\overrightarrow K\cdot\overrightarrow B}{\overrightarrow K\cdot\overrightarrow A}\right),$$ | (13) |
where
| $$\overrightarrow A=\begin{bmatrix}(x_k-x_{k-1})/\sqrt2\\(x_k-x_{k+1})/\sqrt2\\y_{k-1}\\y_k\\y_{k+1}\end{bmatrix}$$ |
일반화된 해석해의 수식 유도과정도 부록에 소개하였다.
IV. 성능지표
잡음이 존재하는 상황에서 추정 파라미터(unknown parameter, θ)의 값은 확률적인 측면에서 고려되어야 한다. 즉, 특정 알고리즘에 의해서 추정된 값 은 참값 θtrue을 얼마나 정확하게 추정하며(bias), 얼마의 오차(Mean-Square Error, MSE)을 가지는가를 고려하고 이를 종합하여 알고리즘을 평가해야 한다.
만약 특정 알고리즘에 의한 가 참값과 비교하여 확률적으로 편중되지 않았다면, 이러한 알고리즘이 가질 수 있는 이론적인 최저 분산값을 크래머 라오 하한(Cramer-Rao Lower Bound, CRLB)라고 한다. 즉, 수식으로 표현하면 의 관계를 갖는다. 이를 이용하면 복수의 알고리즘 성능을 상호 비교하기 용이하다. CRLB는 정칙조건[regularity condition, Eq. (14)]을 만족하는 PDF p(n;θ)에 대해서 Fisher Information matrix을 구함으로 계산할 수 있다.
| $$E\left[\frac{\partial\;\ln\;p(n;\theta)}{\partial\theta}\right]=0\;\;\;\;\mathrm{for}\;\;\;\;\;\forall\theta.$$ | (14) |
WGN 상황에서 정현파에 대한 정보(진폭, 주파수, 위상)가 없는 경우, SNR에 따른 주파수 추정기()의 CRLB 는 아래의 식과 같다.[13]
| $$var\left({\widehat f}_{norm}\right)\geq\frac{12}{(2\pi)^2\eta\;N(N^2-1)},$$ | (15) |
여기서 η는 SNR 이다. Eq. (15)에서 알 수 있듯이 주파수에 대한 CRLB는 SNR이 증가하면 감소하고 데이터 길이에 매우 민감함을 알 수 있다(1/N3에 비례하여 Lower bound 감소).
본 논문에서는 주파수 추정에 있어서 각 알고리즘에 의한 편의(bias) 및 오차(MSE) 정도와 함께 프레임별 추정값의 분산(variance) 정도를 이용하여 성능을 평가하였다.
V. 모의실험
각 알고리즘의 주파수 추정 성능을 비교하기 위하여 성능에 영향을 미치는 다양한 요소들을 변화시키며 실험을 수행하였다. 매 실험 조건마다 1,000번의 Monte-Carlo 시뮬레이션을 수행할 수 있도록 0 < fnorm< 0.5 구간에서 임의로 주파수를 생성(uniform distribution)하였다. 각 주파수 신호는 Eq. (1)에 따라 4,000개의 샘플 데이터(N = 4,000)를 생성하고, 윈도우 함수에 의해 프레임(frame)을 형성한다. 알고리즘의 성능은 주파수별 프레임들의 전체 평균에 의해서 구한다.
우선 잡음이 없는 상황에서 윈도우 종류 및 윈도우 길이(NFFT)에 따른 성능을 비교하였다(Fig. 3). 모든 알고리즘은 비교 기준이 되는 DFT 대비 향상된 성능을 보인다. 특히 Cedron 은 rectangular 윈도우 함수를 사용할 경우 앞서 유도한 바와 같이 정확한 해(주파수)를 구할 수 있다는 사실을 확인할 수 있다[Fig. 3(a)의 Cedron에 대한 MSE, Variance 결과값이 0]. 하지만 그 밖의 윈도우 함수에 대해서는 다른 주파수 추정 알고리즘과 비슷한 성능을 보인다. 이는 윈도우 함수를 적용하는 과정에서 매 프레임의 처음과 끝 부분에 발생한 원신호의 왜곡으로 인하여 Zk*, Zk*-1, Zk*+1 등이 변화되기 때문이다. 알고리즘 및 윈도우 함수에 따른 성능 측정 실험에서 NFFT 가 증가 할수록 각 성능지표가 공통적으로 향상되는 것을 확인 할 수 있다. 편의 지표의 경우, NFFT 증가에 따른 편의 감소 경향을 확인할 수 있으나 일부 경우에 대해서는 경향성과는 상관없이 편의가 0이 되는 경우도 존재하였다. 이는 임의로 선정된 주파수와 주파수 해상도의 정수배 상관관계에 기인한 것으로 추정되나 엄밀한 통계적 확인을 위하여 향후 추가적인 연구가 필요하다. 분산 지표는 신호 자체의 주파수 떨림 특성인 jitter와 구분되기 위해서 가능한 작은 값을 가져야 한다.
잡음이 존재하는 상황에서 NFFT을 256로 고정하고 –10 dB ≤ SNR ≤ 20 dB에 대한 성능 평가를 수행하였다(Fig. 4). SNR이 증가함에 따라 편의와 오차가 감소하는 사실을 확인할 수 있다. 편의의 경우 3 DFT 기반의 4개 알고리즘(DFT, Jacobsen, Candan, Cedron)의 결과 경향이 유사하고 위상 정보를 이용하는 reassignment 방법과 구별이 되었다. 오차(MSE) 의 경우 낮은 SNR 상황에서는 모든 주파수 추정 알고리즘이 DFT 와 유사한 성능을 나타내었다. 하지만 약 -3 dB SNR을 기점으로 DFT는 성능 포화(saturation) 현상을 보이는 반면에 다른 알고리즘들은 지속적으로 향상되는 경향을 나타내거나[Fig. 4(a)], DFT보다 좋은 성능에서 포화현상을 나타내었다 [Fig. 4(b) ~ (d)]. 비교 분석한 알고리즘 중 Cedron 알고리즘이 가장 좋은 성능을 나타내었다.
마지막으로 각 알고리즘에 적용된 윈도우 함수에 따른 성능을 살펴보면, 알고리즘에 상관없이 rectangular 함수가 가장 좋은 성능을 보여준다(Fig. 5). DFT의 경우 약 -3 dB SNR 이상에서는 윈도우 함수에 따른 성능 차이가 나지 않았다.
수행된 모든 경우의 실험결과를 종합하면 Cedron 과 rectangular 윈도우 함수 조합이 가장 좋은 성능을 나타내고 있다. 이를 보다 명확하게 비교하기 위하여 성능이 비교적 안정된 경향을 보이는 SNR 6 dB 상황에서의 MSE 실험결과를 정리하면 Table 1과 같다. 각 수치는 알고리즘과 윈도우 함수가 조합되었을 경우의 실험 결과이다.
반면 연산량 측면에서는 Cedron 이 3 DFT 기반의 알고리즘 중 가장 많은 연산을 필요로 한다. reassignment의 경우 실험에 주어진 윈도우 함수, h (t)뿐만 아니라 t·h (t), dh (t) / dt 윈도우에 의한 DFT 결과를 조합하기 때문에 비교 대상 알고리즘 중 가장 많은 연산을 필요로 한다.
Table 1. Selected MSE result at 6 dB SNR with 256 NFFT. × 10^(-6)
| DFT | Jacobsen | Candan | Reassign | Cedron | |
| Rect. | 1.111 | 0.008 | 0.008 | 0.008 | 0.006 |
| Hamm. | 1.112 | 0.239 | 0.239 | 0.248 | 0.070 |
| Han. | 1.114 | 0.295 | 0.295 | 0.305 | 0.099 |
| Blackman. | 1.164 | 0.423 | 0.423 | 0.432 | 0.184 |
VI. 결 론
본 논문에서는 다양한 조건에서 정현파 신호에 대한 주파수 추정 알고리즘의 성능을 편의와 오차, 분산 측면에서 비교 분석 하였다. DFT 연산이 가지고 있는 단점을 극복하고 고해상도 주파수 정보를 얻기 위한 방법으로 총 4가지 알고리즘이 상호 비교되었으며, 실험 결과 선형대수의 기법을 이용하여 해석해를 제시한 Cedron 알고리즘이 좋은 성능을 보였다. 향후 이번 연구에서 분석하지 못했던 주파수 구간에 따른 성능 분석과 함께 간섭이 존재하는 상황에서 윈도우 함수의 영향력과 이에 따른 알고리즘 성능을 분석 할 예정이다. 더불어 다수의 주파수 성분이 존재하는 복합 신호를 분석하기 위해서 알고리즘에 대한 개선 및 실응용 연구를 진행할 예정이다.
부 록
Cedron의 수식 유도에 앞서 잡음이 없는 순수 정현파 입력신호 xn에 대한 DFT 를 아래와 같이 정의한다.
| $$\begin{array}{l}x_n=M\;\cos\left(2\pi\;f_{norm}\;n+\varnothing\right)\\\;\;\;\;=M\;\cos\left(2\pi\;\frac fNn+\varnothing\right).\end{array}$$ | (A.1) |
| $$\begin{array}{l}where\;f\;is\;\left[cycles/frame\right],\;0\;<\;f\;<\;N/2\;\;and\\\;\;\;\;\;\;\;\;\;\;\;\;N\;is\;\left[samples/frame\right]\end{array}$$ |
이라 정의하면 Eq. (A.1)은 오일러 등식을 이용하여 아래와 같이 정리된다.
| $$x_n=M\frac{e^{i(\alpha n+\varnothing)}+e^{-i(\alpha n+\varnothing)}}2.$$ | (A.2) |
xn에 대한 DFT 는
| $$\begin{array}{l}Z_n=\sum_{n=0}^{N-1}\;x_ne^{-i2\pi\frac kNn}\;\;\;\;(Let\;\beta_k=k\frac{2\pi}N)\\\;\;\;\;\;=\sum_{n=0}^{N-1}\;x_ne^{-i\beta_kn}.\end{array}$$ | (A.3) |
Eq. (A.3)에서 Rk = e-iβk는 복소평면 단위원을 만족하는 k 번째 해와 같다. Eq. (A.2)를 Eq. (A.3)에 대입하여 정리하고, 등비급수(geometric series)을 이용하여 표현하면 (i.e., ) 다음과 같이 표현할 수 있다.
이후 분모, 분자항에 대한 대수정리 및 오일러등식을 통한 간략화를 거쳐 최종적으로 다음과 같이 정리된다.
| $$Z_k=\frac M2\left[\frac{U_e^{i\beta k}-V}{\cos(\alpha)-\cos(\beta_k)}\right],$$ | (A.5) |
여기서 U=cos(αN+Ø)-cos(Ø)이고, V=cos(αN-α+Ø)-cos(-α+Ø) 이다.
주파수(fnorm or f)를 찾는다는 의미는 𝛼를 찾는 것과 동일한 작업이 된다(i.e., α=f(2π/N)). Cedron은 Eq. (A.5)에서 cos 함수 argument로 존재하는 𝛼를 구하기 위하여 k의 인근 주파수 빈, k-1, k, k+1을 이용한 연립방정식을 세운다.[9] 즉,
한편 이고, 유사하게 이다. Eq. (A.6)의 연립방정식을 이용하면 다음 두 개의 식을 얻을 수 있다(미지수 V 의 제거).
| $$\begin{array}{l}2\left[\cos\;\alpha\;(Z_k-Z_{k-1})-(\cos\;\beta_kZ_k-\cos\;\beta_{k-1}Z_{k-1})\right]\\\;\;\;\;\;=M\;Ue^{i\beta_k}(1-R_1),\end{array}$$ | (A.7) |
| $$\begin{array}{l}2\left[\cos\;\alpha\;(Z_k-Z_{k+1})-(\cos\;\beta_kZ_k-\cos\;\beta_{k+1}Z_{k+1})\right]\\\;\;\;\;\;=M\;Ue^{i\beta_k}(1-1/R_1).\end{array}$$ | (A.8) |
다음으로 Eq. (A.7)에서 Eq. (A.8)을 나누어주면 미지수 M 과 U을 제거할 수 있다.
| $$\frac{\cos\;\alpha\;(Z_k-Z_{k-1})-(\cos\;\beta_kZ_k-\cos\;\beta_{k-1}Z_{k-1})}{\cos\;\alpha\;(Z_k-Z_{k+1})-(\cos\;\beta_kZ_k-\cos\;\beta_{k+1}Z_{k+1})}=-R_1.$$ | (A.9) |
마지막으로 Eq. (A.9) 을 cos 𝛼에 정리하면,
을 얻을 수 있다. 우변에 존재하는 항은 이미 값을 알고 있는 항들이기 때문에 𝛼를 구하기 위해서는 cos 역함수를 이용하면 된다. 분석하고자 하는 신호에 잡음이 존재한다면 k는 최대 DFT 크기를 갖는 주파수 빈을 이용한다.
일반화된 Cedron 수식은 Eq. (A.3) DFT 결과를 실수부와 허수부로 구분하는 것을 시작으로 유도된다.[10] 즉, Zk = xk + iyk 이다. 이를 Eq. (A.5)에 적용하면
| $$x_k=\frac M2\left[\frac{U\cos(\beta_k)-V}{\cos(\alpha)-\cos(\beta_k)}\right],$$ | (A.11) |
| $$y_k=\frac M2\left[\frac{U\sin(\beta_k)}{\cos(\alpha)-\cos(\beta_k)}\right].$$ | (A.12) |
분모, 분자를 상호 곱하면
| $$\cos(\alpha)\cdot x_k-\cos(\beta_k)\cdot x_k=\frac M2\left[U\;\cos\;(\beta_k)-V\right],$$ | (A.13) |
| $$\cos(\alpha)\cdot y_k-\cos(\beta_k)\cdot y_k=\frac M2U\;\sin\;(\beta_k),$$ | (A.14) |
동일한 전개를 주파수 빈 k-1과 k+1에도 수행을 하면, DFT 실수부와 관련된 3개의 식[Eq. (A.13) 및 그 변형]을 연립방정식 형태로 얻을 수 있다. 미지수 V을 제거하기 위해 xk 관련 항을 xk-1관련 항과 xk+1 관련 항으로 각각 빼주면 2개의 등식을 얻는다.
| $$\begin{array}{l}\cos(\alpha).(x_k-x_{k-1})-\left[\cos(\beta_k)\cdot x_k-\cos(\beta_{k-1})\cdot x_{k-1}\right]\\\;\;\;\;=\frac M2U\left[\cos(\beta_k)-\cos(\beta_{k-1})\right],\end{array}$$ | (A.15) |
| $$\begin{array}{l}\cos(\alpha).(x_k-x_{k+1})-\left[\cos(\beta_k)\cdot x_k-\cos(\beta_{k+1})\cdot x_{k+1}\right]\\\;\;\;\;=\frac M2U\left[\cos(\beta_k)-\cos(\beta_{k+1})\right].\end{array}$$ | (A.16) |
Eqs. (A.15), (A.16) 및 허수부 관련 3개의 식[Eq. (A.14) 및 그 변형]을 벡터 형태로 표현하게 되면 Eq. (A.17) 이 되고, 이때 , 그리고 는 Eq. (8)에서 정의된 것과 같다.
| $$\cos(\alpha).\overrightarrow A-\overrightarrow B=\frac{MU}2\overrightarrow C.$$ | (A.17) |
Eq. (A.17)을 만족하는 𝛼는 에 직교(orthogonal) 하고 및 에는 직교하지 않는 임의의 벡터 을 선정함으로서 찾을 수 있다. 즉,
이를 통하여 미지수 M 과 U을 제거할 수 있다. 실제 분석 신호는 잡음을 가지고 있기 때문에 와 에 해당하는 항을 각각 로 표현하는 것이 더 적절하다. error 와 관련된 항은 을 조절함으로서 최소화할 수 있는 것이 아니다. 대신 신호 자체()의 크기를 최대화 할 수 있도록 를 설계 할 수 있다(신호와 최대한 동일 방향으로 설계). 앞서 와는 직교해야 한다는 조건을 고려한다면,
| $$\overrightarrow K=\overrightarrow D-w\overrightarrow C\;\;\;(w\;is\;scalar).$$ | (A.19) |
| $$\overrightarrow K\cdot\overrightarrow C=\overrightarrow D\cdot\overrightarrow C-w\overrightarrow C\cdot\overrightarrow C=0.$$ | (A.20) |
| $$w=\frac{\overrightarrow D\cdot\overrightarrow C}{\overrightarrow C\cdot\overrightarrow C},$$ | (A.21) |
여기서 이다. 따라서
Eq. (A.22)를 Eq. (A.18)에 대입하여 최종적으로 cos(𝛼)를 구할 수 있다.
