구형파의 분석 및 합성

[신호 및 시스템]이라는 과목을 공부하는 도중에 푸리에 급수에 대해서 찾아보다가 어떤 교수님께서, 

푸리에 급수에 대한 신호함수를 파이썬 스크립트로 구현한 부분이 있어 매우 흥미로워 가져와 보았습니다.

 

 

※다음은 주기가 1이고, 듀티사이클이 0.5인 구형파의 푸리에 급수 분석과 합성을 실행하는 파이썬 스크립트이다.

import numpy as np
import matplotlib.pyplot as plt

SMALL = 1.e-5
dt = 1/1000.
t = np.arange(0,1.,dt)
T0 = 1.
f0 = 1/T0
w0 = 2*np.pi * f0
xt = 1. * (t < 0.5)


Nmax = 31 #Maximum number of harmonics
# Calculation of harmonics
a = np.zeros(2*Nmax+1).astype('complex')
magResponse = np.zeros(2*Nmax+1).astype('float')

krange = np.arange(-Nmax, Nmax+1)
for k in krange:
    x1 = xt * np.exp(-1j * w0 * t * k)
    a[k-Nmax-1] = dt * f0 * sum(x1)

#Synthesis
xN = np.zeros(len(t)).astype('complex')
for k in krange:
    xN = xN + a[k-Nmax-1] * np.exp(1j*w0*k*t)


magResponse = np.abs(a)
phaseResponse = np.angle(a)


phaseResponse = phaseResponse * (magResponse > SMALL)
phaseResponse = phaseResponse - np.pi \
                * (phaseResponse > np.pi/2)
phaseResponse = phaseResponse + np.pi \
                * (phaseResponse < -np.pi/2)


ax_org = plt.subplot(411)
plt.plot(t, xt)
ax_org.set_ylim(-0.1, 1.1)

ax_mag = plt.subplot(412)
ax_mag.stem(krange, magResponse)
ax_mag.set_xlim(-Nmax, Nmax)

ax_phase = plt.subplot(413)
ax_phase.stem(krange, phaseResponse)
ax_phase.set_xlim(-Nmax, Nmax)
ax_phase.set_yticks([-np.pi/2, 0, np.pi/2])
ax_phase.set_yticklabels(["$-\\frac{\pi}{2}$", \
                          "0", "$\\frac{\pi}{2}$"])

plt.subplot(414)
plt.plot(t, xN.real)
plt.show()

#xt = 2 * t * (t <= 0.5) + (2- 2*t) * (t>0.5)

실행결과>

 

연습문제5.1)

앞의 파이썬 스크립트에서 줄번호 31~35사이를 주석처리하고, 스크립트를 실행해보자. 크기 스펙트럼에는 큰 차이가 없지만 위상 스펙트럼이 무질서하게 나오는 것처럼 보일 것이다. 그 이유가 무엇인지 생각해보자.

import numpy as np
import matplotlib.pyplot as plt

SMALL = 1.e-5
dt = 1/1000.
t = np.arange(0,1.,dt)
T0 = 1.
f0 = 1/T0
w0 = 2*np.pi * f0
xt = 1. * (t < 0.5)


Nmax = 31 #Maximum number of harmonics
# Calculation of harmonics
a = np.zeros(2*Nmax+1).astype('complex')
magResponse = np.zeros(2*Nmax+1).astype('float')

krange = np.arange(-Nmax, Nmax+1)
for k in krange:
    x1 = xt * np.exp(-1j * w0 * t * k)
    a[k-Nmax-1] = dt * f0 * sum(x1)

#Synthesis
xN = np.zeros(len(t)).astype('complex')
for k in krange:
    xN = xN + a[k-Nmax-1] * np.exp(1j*w0*k*t)


magResponse = np.abs(a)
phaseResponse = np.angle(a)


""" <!-- 이 부분을 주석처리합니다. -->
phaseResponse = phaseResponse * (magResponse > SMALL)
phaseResponse = phaseResponse - np.pi \
                * (phaseResponse > np.pi/2)
phaseResponse = phaseResponse + np.pi \
                * (phaseResponse < -np.pi/2)
"""


ax_org = plt.subplot(411)
plt.plot(t, xt)
ax_org.set_ylim(-0.1, 1.1)

ax_mag = plt.subplot(412)
ax_mag.stem(krange, magResponse)
ax_mag.set_xlim(-Nmax, Nmax)

ax_phase = plt.subplot(413)
ax_phase.stem(krange, phaseResponse)
ax_phase.set_xlim(-Nmax, Nmax)
ax_phase.set_yticks([-np.pi/2, 0, np.pi/2])
ax_phase.set_yticklabels(["$-\\frac{\pi}{2}$", \
                          "0", "$\\frac{\pi}{2}$"])

plt.subplot(414)
plt.plot(t, xN.real)
plt.show()

#xt = 2 * t * (t <= 0.5) + (2- 2*t) * (t>0.5)

 

실행결과>

※ 삼각파의 분석 및 합성 (추가)

앞의 파이썬 스크립트 중에서 10번째 줄을 다음과 같이 수정하면 삼각파의 신호에 대해서도 분석과 합성을 할 수 있다.

import numpy as np
import matplotlib.pyplot as plt

SMALL = 1.e-5
dt = 1/1000.
t = np.arange(0,1.,dt)
T0 = 1.
f0 = 1/T0
w0 = 2*np.pi * f0
"""
xt = 1. * (t < 0.5)
"""
xt = 2 * t * (t <= 0.5) + (2- 2*t) * (t>0.5)

Nmax = 31 #Maximum number of harmonics
# Calculation of harmonics
a = np.zeros(2*Nmax+1).astype('complex')
magResponse = np.zeros(2*Nmax+1).astype('float')

krange = np.arange(-Nmax, Nmax+1)
for k in krange:
    x1 = xt * np.exp(-1j * w0 * t * k)
    a[k-Nmax-1] = dt * f0 * sum(x1)

#Synthesis
xN = np.zeros(len(t)).astype('complex')
for k in krange:
    xN = xN + a[k-Nmax-1] * np.exp(1j*w0*k*t)


magResponse = np.abs(a)
phaseResponse = np.angle(a)


"""
phaseResponse = phaseResponse * (magResponse > SMALL)
phaseResponse = phaseResponse - np.pi \
                * (phaseResponse > np.pi/2)
phaseResponse = phaseResponse + np.pi \
                * (phaseResponse < -np.pi/2)
"""


ax_org = plt.subplot(411)
plt.plot(t, xt)
ax_org.set_ylim(-0.1, 1.1)

ax_mag = plt.subplot(412)
ax_mag.stem(krange, magResponse)
ax_mag.set_xlim(-Nmax, Nmax)

ax_phase = plt.subplot(413)
ax_phase.stem(krange, phaseResponse)
ax_phase.set_xlim(-Nmax, Nmax)
ax_phase.set_yticks([-np.pi/2, 0, np.pi/2])
ax_phase.set_yticklabels(["$-\\frac{\pi}{2}$", \
                          "0", "$\\frac{\pi}{2}$"])

plt.subplot(414)
plt.plot(t, xN.real)
plt.show()

 

실행결과>

참고 및 출처:

http://contents.kocw.net/KOCW/document/2014/Hallym/parkseophyeng/14.pdf

 

 

(아직 완전히 이해를 하지는 못했지만, 매우 흥미로웠습니다. 좀 더 공부를 해봐야겠네요^^)