06-5. 엑셀을 이용한 DFT 계산(2/4)

이번 페이지에서는 정현파에 대한 신호값을 생성하고, 이 신호값에 대한 DFT를 엑셀을 이용해서 수작업으로 계산해볼 것이다. (수작업이라고 표현한 이유는, 엑셀에 내장되어 있는 FFT함수를 써서 계산하지 않고, DFT의 수식을 있는 그대로 사용해서 하나씩 계산해본다는 의미)

 

엑셀 자체의 기능보다는, 신호값에 대해서 어떻게 DFT 수식이 적용되는지에 집중해서 보면 되겠다. 기본적인 방법은 앞 페이지에서 4개의 값에 대한 DFT 구하는 것과 동일.

 

---

[예제 2] 단순 정현파에 대한 DFT 계산

이 예제에 대한 엑셀은, 이 페이지 하단부에 첨부된 "DFT예제_2.xlsm"의 Sheet "2"에서 확인할 수 있음

 

가장 단순한 정현파 신호를 만들자. 하나의 코사인 파형이면, 가장 단순하다고 할 수 있겠다. 

 

• 코사인 파형값을 만들기 위한 기본 수식: $\cos{(2\pi f t)}$  --> f:주파수,   t: 시간
• 주파수(f): 1/64  --> 1초에 $\frac{1}{64}$개 주기의 신호. 즉, 64초마다 1개 주기의 코사인 신호 (주기가 64초)
• 시간 샘플링 단위: 1초 --> t의 값을 1초 간격으로 0~63초까지의 값을 구해놓으면, 코사인 1주기 파형값이 생성될 것임
•  

(1) 정현파 값 y(n) 생성 

 

위와 같은 파형이 생성되도록 y(n)값을 생성한다. (ROUND 함수를 써서 소수점 2자리까지만 표시되도록 했다. 숫자가 길어지면 가독성이 떨어져서 줄여버린 것이고, ROUND 처리하지 않아도 무방)

• 수식:  =ROUND(COS(2*PI()*(1/64)*B14),2)

 

 

(2) W(n,k) 값 계산

 

W(n,k)의 수식은 $e^{-i\frac{2\pi}{N}kn}$이다. 

 

엑셀에서 "E14"셀에 수식을 넣고, k=63까지, n=63까지 자동채우기로 "E14"수식을 적용한다.

 

• E14셀 수식: =CTrim(IMEXP(COMPLEX(0,-2*PI()*E$13*$B14/$C$11)))

엑셀에서 자동채우기를 할 때, 절댓값 기호 $를 어떻게 부여해야 하는지 헷갈릴 수 있다. Row에만 붙일지 Column 주소에만 붙여야 할지 등.

쉬운 방법은, 자동채우기를 할 때 "변하지 않을 것이 무엇인가"를 생각하는 것이다.

예를 들어 E14셀의 수식에서 k의 값이 들어 있는 E13에 대해서, E14셀을 오른쪽으로 쭈욱 긁어나갈 때 {E13, F13, G13, ...}의 값들이 사용돼야 할 것이다. 그렇다면 여기서 변하지 않는 것은 무엇일까를 생각하면 된다는 의미. 여기서 변하지 않는 것은 "13"다. 따라서 "13"만 절대 주소로 처리해서 "E$13"로 처리하면 되는 것이다.

n의 값인 B14에 대해서는 "E14"의 값을 아래로 자동채우기할 때 변하게 되는데, {B14, B15, B16, ... }이렇게 적용되야하고, 이때 변하지 않는 것은 "B"이기에 "$B14"로 하면 되는 것.   

(3) Y(k)의 계산

 

Y(k)는 y(n)과 W(n,k)의 합성곱이다. 여기서 $W(n,k)= e^{-i\frac{2\pi}{N}kn}$

 

$$ Y(k) = \sum _{n=0}^{N-1}{y(n)W(n,k)}$$

 

엑셀에서 y(n)과 W(n,k)값을 구했으니, 두 값에 대한 합성곱을 계산하면 되겠다. 

합성곱은 이전 페이지에서 작성했던 CSumProduct 함수를 이용하면 되겠다. 

 

• "E78"의 수식: =CSumProduct($C14:$C77,E14:E77)

수식에서 y(n) 값의 범위에서 "C열"을 고정하기 위해 절대 주소기호 "$" 쓴 것에 유의.

"E78"의 수식을, 오른편으로 자동 채우기 해서 Y(63)까지 구한다.

 

(4) Y(k)에 대한 실수부, 허수부, 크기 구하기

 

복소수에 대한 실수부는 IMREAL 함수를, 허수부 값은 IMAGINARY 함수를 이용해서 뽑아낼 수 있다. 

크기는 실수부의 제곱과 허수부의 제곱에 대한 제곱근이다.

 

• "E79" 셀:  =IMREAL(E78)
• "E80" 셀: = RIREAL(E78)
• "E81" 셀: =ROUND(SQRT(E79^2+E80^2),2)

Y(k) 값과 마찬가지로, 오른쪽으로 자동채워넣기해서, Y(63)에 대한 실수부, 허수부, 크기를 구한다.

 

(5) Y(k)에 대한 값을 그래프로 그려보기

 

Y(k)값만 봐서는 한눈에 파악이 안 된다. Y(k)에 대해 그래프를 그려보면 값에 대한 분포 파악이 가능하다. 여기서 중요한 것은 특히 Y(k)의 크기 값이다.

 

(6) DFT 값인 Y(k)에 대한 해석

 

신호값인 y(t)가 주파수 $f=\frac{1}{64}$인 정현파였고, 1초 단위로 시간 샘플링을 한 값이었다. 즉, 원래의 주파수는 $\frac{1}{64}Hz$이고, 샘플링 주파수 $f_s=1$이다. 따라서, n의 간격은 1초이고, n=63까지의 값을 구한 것이기에, y(n)은 0~63초까지의 신호 값을 구해서 사용했다.

 

$f=\frac{1}{64}$이고 $f_s=1$이기에 매우 촘촘하게 샘플링한 것이다. (1초에 1/64 주기만 진행하는, 즉 64초에 1주기인 신호에 대해서, 1초마다 샘플링을 한 것이기에, 1주기 신호를 64번 샘플링한 셈이다.)

 

주파수에 대한 이산 샘플링 주기값을 나타내는 $N=64$이다. 그리고 n 사이 간격은 1초이기에, k의 간격은 $\frac{1}{64}$이 된다. 따라서, $k=1$일 때의 DFT값인 $Y(1)$은 실제는 $\frac{1}{64}Hz$에 대한 값이 된다. 

 

위 그림 중 Y(k)_magnitude에 대한 그래프를 보면, k=1일 때 값이 크게 나타남을 볼 수 있다. 즉, 실제 신호에서의 $\frac{1}{64}Hz$에서의 값이 가장 큰 것이고, 이 신호는 주파수가 $\frac{1}{64}Hz$인 신호임을 유추할 수 있다.

이 부분이 쉽게 이해가 안될 수 있다. (나의 경우 그랬다.)

시간 t에 대한 신호에서, 이를 시간에 대한 이산 샘플링을 해서 y(n)을 구하고, 또한 주파수에 대한 이산 샘플링을 해서 k 간격으로 나타내다 보니, 시간 t일 때의 주파수인 Hz와, 이산 시간 n과 이산 주파수 단위 k에 대한 변환이 쉽게 이해되지 않을 수 있다. 특히, N의 값에 의해서 변하는 k 간격 단위가 쉽게 이해되지 않는다.

그래서, DFT에 대한 공식만을 유도해서 넘어가는 것이 아니고, 이전 페이지에서 Y(k) 값을 손으로 구해보고, 이 페이지에서 엑셀로도 풀어보고 하는 것이다. 

다음 페이지의 예제 하나를 더 따라가 보면, 이해될 수 있을 것이다.

   

---

이 페이지에서 사용된 엑셀 파일:

DFT예제_2 - update.xlsm

0.10MB

---

2022.08.10 수정사항

 

"Leo"님이 지적하신 CSumProduct 함수의 결괏값표출 오류 수정했습니다.

 

오류 발생 코드

• z.re=31.99, z.im=-0.0000000000001인 경우 im=0인데 반해 z.im < 0 임
• 따라서 CString = re & "" & im & "i" 루틴으로 가게되고, 이것은 31.99 & "" & 0 & "i" = 31.99i 가 되어, 이상하게 변하게 됨. 즉, 원래는 "31.99 + 0i" 인데, 이상하게 표출되는 것
•  아래와 같이 수정하여 해결함

---

-끝-

이전글: 06-4. 엑셀을 이용한 DFT 계산(1/4)

다음글: 06-6. 엑셀을 이용한 DFT 계산(3/4)

다음다음글: 06-7. 엑셀을 이용한 DFT 계산(4/4)