연속 시간 푸리에 급수 실습 with MATLAB

 연속 시간 푸리에 급수로 주기 신호를 분해해보자. 직접 해보면 좀 더 이해가 쉽다. 이 번에는 사각파의 푸리에 급수를 유도하고 많이 쓰이는 공학 계산 툴인 MATLAB을 이용한다.

 

목차

• 1. 푸리에 급수의 복소 표현
• 2. 사각파의 표현
• 3. 연속 시간 푸리에 급수의 유도
• 4. 주기 함수의 평균
• 5. MATLAB 코드
• 6. 코드 실행 결과
• 7. 심심하니까 평균 이동

1. 푸리에 급수의 복소 표현

 푸리에 급수의 복소 표현은 아래와 같다. 자세한 내용은 예전 글 '연속 시간 푸리에 급수'에서 확인해보자.

$$\begin{align} f(t) &= \sum_{n=-\infty}^{\infty} F_n e^{j2\pi nf_0t} \cdots (1)\\ F_n &= \frac{1}{T_0} \int_{\frac{-T_0}{2}}^{\frac{T_0}{2}} f(t) e^{-j2\pi nf_0t} dt \cdots (2) \end{align}$$

 

2. 사각파의 표현

 상한은 값 A를 갖고 하한은 값 B를 갖는 사각파를 생각해보자. 여기에서 A가 0.5, B가 -0.5이면 평균값 0을 지나는 사각파가 그려진다. 이 사각파를 그림 1에 나타냈다.

$$\begin{align} f(t) = A \qquad& 0 < t < \frac{T_0}{2} \\ f(t) = B \qquad& \frac{T_0}{2} < t < T_0 \\ \end{align}$$

 

$$T_0 = 2.0 \qquad f_0 = \frac{1}{T_0} = 0.5 \qquad A = 0.5 \qquad B = -0.5 \qquad$$

그림 1. 주기 2를 갖는 사각파

 

3. 연속 시간 푸리에 급수의 유도

 푸리에 급수는 식 (1)로부터 다음과 같이 구할 수 있다.

$$\begin{align} F_n &= \frac{1}{T_0} \int_{0}^{T_0} f(t) e^{-j2\pi nf_0t} dt \\ &= \frac{1}{T_0} \int_{0}^{\frac{T_0}{2}} A e^{-j2\pi nf_0t} + \frac{1}{T_0} \int_{\frac{T_0}{2}}^{T_0} B e^{-j2\pi nf_0t} \\ &= -\frac{A}{j2\pi nT_0f_0} \left[ e^{-j2\pi nf_0t} \right]_{0}^{\frac{T_0}{2}} -\frac{B}{j2\pi nT_0f_0} \left[ e^{-j2\pi nf_0t} \right] _{\frac{T_0}{2}}^{T_0} \\ &= -\frac{A}{j2\pi n} \left[ e^{-j\pi n} - 1\right] -\frac{B}{j2\pi n} \left[ e^{-j2\pi n} - e^{-j\pi n} \right] \end{align}$$

 

 따라서 사각파의 푸리에 급수 $F_n$은 다음과 같이 정리할 수 있다.

$$F_n = \frac{A \left( 1 - e^{-j\pi n} \right) + B \left(e^{-j\pi n} - e^{-j2\pi n} \right)}{j2\pi n}$$

 

 여기에서 $exp(-j2\pi n) = 1$이므로 이러쿵 저러쿵 정리하면 아래와 같다.

$$F_n = \frac{ \left( A-B \right) \left( 1 - e^{-j\pi n} \right)}{j2\pi n}$$

 

 이제 푸리에 급수의 계수 $F_n$을 구했으니 다시 식 (1)에 대입해 급수를 구한다.

$$f(t) = \sum_{n=-\infty}^{\infty} F_n e^{j2\pi nf_0t}$$

 

 $f(t)$의 우변은 복소수 $F_n$에 복수 지수함수 $e^{j2\pi nf_0t}$를 곱했으니 $f(t)$의 치역은 복소수가 된다. 결과가 복소수라니 그럼 어떻게 그래프를 그려야 할까? 결론부터 말하면 $f(t)$의 실수부만 취하면 된다. 따라서 우리가 구하려는 $f(t)$는 복소수 $f(t)$의 실수부 $\operatorname{Re} f(t)$가 된다.

 

 푸리에 급수를 복소 지수 함수로 표현할 때 오일러 공식을 이용해 아래와 같은 관계식을 유도했다.

$$F_n = \frac{1}{2} \left( a_n - jb_n \right)$$

 

 하나의 항만 놓고 생각해보면 아래와 같이 정리된다.

$$\begin{align} f_n(t) &= F_n e^{j2\pi nf_0t} \\ &= \frac{1}{2} \left( a_n - jb_n \right) \left[ \cos(2\pi nf_0t) + j\sin(2\pi nf_0t) \right]\\ &= \frac{1}{2} \left[ a_n \cos(2\pi nf_0t) + b_n sin(2\pi nf_0t) \right] + \frac{j}{2} \left[a_n\sin(2\pi nf_0t) -b_n\cos(2\pi nf_0t) \right] \end{align}$$

 

 그러므로 우리가 원하는 푸리에 급수는 아래와 같다.

$$\operatorname{Re} f(t) = \frac{1}{2} \sum_{n=-\infty}^{\infty} \left[ a_n \cos(2\pi nf_0t) + b_n \sin(2\pi nf_0t) \right]$$

 

 그런데 $F_{-n} = F_{n}$인 것을 알고 있으니 우리는 절반만 구해서 2배를 해주면 계산량도 줄일 수가 있다.

$$\operatorname{Re} f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left[ a_n \cos(2\pi nf_0t) + b_n sin(2\pi nf_0t) \right]$$

 

 사실 이 함수는 기함수이기 때문에 $\sin$ 급수만 있고 $\cos$ 급수는 모두 0이 되어 사라진다. 이런 단순한 경우에는 복소 지수 함수를 이용하는 것보다 그냥 푸리에 적분에서 구한 푸리에 사인 변환을 이용하는 것이 더 빠르다. 그러나 일반적인 경우에는 복소 지수 함수를 이용하는 것이 편리하다.

 

4. 주기 함수의 평균

푸리에 급수의 첫 번째 항은 $\frac{a_0}{2} = F_0$이며 주기 함수의 평균이다. 이 $F_0$를 구하기 위해 위에서 구한 $F_n$에 0을 대입해보자. 분모가 $0$이 되어서 계산이 안 된다는 것을 알 수 있다. $F_0$는 $F_n$에서 $n$을 $0$으로 보내는 극한을 통해 구할 수 있다. $F_n$의 분자와 분모가 각각 $n\rightarrow0$일 때 모두 $0$으로 향하므로 로피탈의 정리를 사용할 수 있다. 분자 분모를 $n$에 대해 미분해서 구해 보면 $F_0$는 아래와 같다. 과연 평균이라는 것을 알 수 있다.

$$F_0 = \frac{1}{2} \left( A + B \right)$$

 

5. MATLAB 코드

 이 푸리에 급수를 MATLAB 코드로 작성해보자. 간단하게 아래처럼 반복문으로 작성했다. 평균 $F_0$는 따로 구하고 $f(t)$의 초기값으로 부여 한 뒤 $n$ 번째 $F_n$을 계산하고 다시 $f_n(t)$를 계산했다. 마지막으로 복소수 $t$ 벡터의 실수부만을 취해 각 항의 $f(t)$를 구했다. 이를 누적하여 항을 늘려 가면서 과연 급수가 원함수 사각파에 가까워지는지 확인해보자.

 

%% Clear
clear all;
close all;
clc;

%% Source Parameter
A = 0.5;
B = -0.5;
T_0 = 2.0; % Fundamental Period
freq_0 = 1/T_0; % Fundamental Frequency
T_s = 0.001; % Sampling increment

%% Source Plot
t = 0:T_s:T_0; % Time vector

f_t = zeros(size(t)); % zero initialization
t_half = find(t==T_0/2, 1); % find index of half period
f_t(1:t_half) = A;
f_t(t_half+1:end) = B;

%% Continus Time Fourier Series
F_0 = 0.5 * (A + B); % Average of signal
f2_t = F_0;

figure(2)
for n=1:1000
    F_n = (A-B) * (1-exp(-1i*pi*n)) / (1i*2*pi*n);
    f_nt = F_n .* exp(1i*2*pi*n*freq_0*t); % complex f_n(t)
    f2_t = f2_t + 2 * real(f_nt); % cumulative f(t)
    
    plot(t, f_t);
    hold on
    plot(t, f2_t);
    hold off
    axis([-0.1, T_0*1.1, min(A, B)-0.1, max(A, B)*1.2])
    grid on
    
    strN = num2str(n);
    str = strcat('n=[-', strN, ',', strN, ']');
    text(T_0/2, max(A, B)*1.1, str);
    
    text(T_0*0.55, A, 'A');
    text(T_0*0.45, B, 'B'); 
    
    pause(0.00001)
end

 

6. 코드 실행 결과

 그림 2를 보면 뾰족한 불연속부에서 삐죽 튀어나와 잘 수렴하지 못하는 Gibbs 현상을 확인할 수 있다. 항을 늘려가면서 점차 줄어들어면서 수렴하는 듯하지만 다시 나타났다가 사라졌다가 하는 것도 알 수 있다. 아무튼 푸리에 급수로 만든 사각파가 원래 함수를 쫓아가는 것을 확인했다.

 

그림 2. 사각파의 푸리에 급수

 

7. 심심하니까 평균 이동

 사각파를 위로 조금 옮겨서 평균이 0이 아닌 경우를 살펴보자. $A=1, B=0$로 두면 평균이 $0.5$가 된다. 이 사각파의 푸리에 급수를 그래프로 나타내 보면 그림 3과 같다.

 

 당황스럽게도 시간이 $0$일 때 함숫값이 $0$도 아니고 $1$도 아니다. 그래프를 보면 $0.5$에서 시작하고 끝난다는 것을 알 수 있다. 이게 머선129? 아까 평균이 $0$일 때는 마치 $\sin$같아서 깜빡 속아서 신경을 안 썼다. 다시 보니 여기도 불연속부에서 값이 평균값이었던 것이다.

 

  푸리에 급수 편에서 설명했듯이 이런 불연속부에서는 푸리에 급수의 결과가 좌극한과 우극한의 평균이 된다. 이 사각파에서 $t=0^{-}$일 때 $0$이고 $t=0^{+}$은 $1$이므로 평균인 $0.5$가 나오는 것도 확인해볼 수 있다. 다시 그림 2를 보면 항을 더해나갈 때마다 항상 불연속부에서 평균을 지나는 것을 알 수 있다. 사실 푸리에 급수는 원래 함수가 조각 연속이기만 하면 모든 점에서 좌극한과 우극한의 평균이지만 매끄러운 연속점에서는 좌극한과 우극한이 일치하고 평균이 원래 함수 값이 되는 것이다.

 

그림 3. 평균이 0이 아닌 사각파