코딩일지(2024-05-10)

계측공학에서 배운걸로 5장 예제를 풀어보자!

---

5장의 9,14,19번 문제를 숙제로 내주셨다. 내일 풀어보는 걸로 하고! 오늘은 일단 예제를 풀면서 몸풀기를 해야겠다!

5장 예제 1번

로직함수 정의와 파라미터 설정 코드>

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['axes.unicode_minus'] = False
plt.rc('font', family='Malgun Gothic')

# 주파수 응답과 위상 지연 데이터를 함수로 정의
def gain(f):
    if f == 1000:
        return 200
    elif f == 2000:
        return 220
    # 고조파가 추가될 경우 여기서 이득을 정의할 수 있다.
    return 1

def phase_lag(f):
    if f == 1000:
        return 0
    elif f == 2000:
        return 20 * (np.pi / 180) # degrees to radians
    # 고조파가 추가될 경우 여기서 위상 지연을 정의할 수 있다.
    return 0

# 사각파 생성
def square_wave(t, f):
    return (4 / np.pi) * np.sum([np.sin(2*np.pi*(2*n-1) * f * t) / (2*n-1) for n in range(1, 10)], axis=0)

# 시간 설정
fs = 100000 # 샘플링 주파수
T = 0.002   # 신호 지속 시간
t = np.linspace(0, T, int(fs*T), endpoint=False)

# 입력 사각파 생성
f1 = 1000 # 기본 주파수
input_signal = square_wave(t,f1)

# 출력 신호 계산
output_signal = np.zeros_like(t)
for harmonic in [1, 2]: # 1000Hz와 2000Hz 성분을 고려
    f_h = harmonic * f1
    g = gain(f_h)
    phi = phase_lag(f_h)
    output_signal += g * (4 / np.pi) * np.sin(2 * np.pi * f_h * t + phi) / harmonic

시각화 코드>

# 시각화

plt.figure(figsize=(12, 6))
plt.subplot(2,1,1)
plt.plot(t, input_signal, label='입력 사각파')
plt.title('입력 사각파 (1000 Hz 기본 주파수)')
plt.xlabel('시간 (s)')
plt.ylabel('진폭')
plt.grid(True)
plt.legend()

plt.subplot(2,1,2)
plt.plot(t, output_signal, label='출력 신호', color='orange')
plt.title('출력 신호 (증폭기 이득 및 위상 지연 적용)')
plt.xlabel('시간 (s)')
plt.ylabel('진폭')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

result>

위 코드에서 np.zeros_like()함수가 어떻게 작동하는지 궁금해져서 다음과 같이 코드를 짜보았다.

t2 = np.linspace(0, 10, 5, endpoint=False)
a = np.zeros_like(t2)
a

result>

array([0., 0., 0., 0., 0.])

위 코드를 통해 np.zeros_like()함수는 주어진 범위만큼 0으로 채워주는 역할을 한다는 것을 알 수 있었다.

또한, 위 코드에서 첫번째 사각파 코드를 보다가, 문뜩, 주파수의 변화에 따라 그 양상이 어떻게 변할까? 궁금해져서 다음과 같이 코드를 구현해보았다.

frequency값을 500, 1000, 2000, 4000Hz로 2배수로 설정해주고 결과를 확인해보았다.

plt.figure(figsize=(8, 16))

input_signal1 = square_wave(t,500)
plt.subplot(4,1,1)
plt.plot(t, input_signal1, label='입력 사각파(1st)')
plt.title('입력 사각파 (500 Hz 기본주파수)')
plt.xlabel('시간 (s)')
plt.ylabel('진폭')
plt.grid(True)
plt.legend()

input_signal2 = square_wave(t,1000)
plt.subplot(4,1,2)
plt.plot(t, input_signal2, label='입력 사각파(2nd)')
plt.title('입력 사각파 (1000 Hz 기본주파수)')
plt.xlabel('시간 (s)')
plt.ylabel('진폭')
plt.grid(True)
plt.legend()

input_signal3 = square_wave(t, 2000)
plt.subplot(4,1,3)
plt.plot(t, input_signal3, label='입력 사각파(3rd)')
plt.title('입력 사각파 (2000 Hz 기본주파수')
plt.xlabel('시간 (s)')
plt.ylabel('진폭')
plt.grid(True)
plt.legend()

input_signal4 = square_wave(t, 4000)
plt.subplot(4,1,4)
plt.plot(t, input_signal4, label='입력 사각파(4th)')
plt.title('입력 사각파 (4000 Hz 기본주파수)')
plt.xlabel('시간 (s)')
plt.ylabel('진폭')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

result>

5장 예제 2번

수식적 해설>

위 수식을 파이썬 넘파이로 구현해보면 다음과 같다.

import numpy as np

# 초기조건 및 상수
T_0 = 75 # T_0: 초기조건
T_inf = 300 # T_inf: 나중조건 (계단형 입력 함수는 75에서 300으로 변한다)
tau = 6 # tau: 시정수
t = 10 # t: 경과시간

# 1차 시스템의 응답 계산
T_t = T_inf + (T_0 - T_inf)*np.exp(-t/tau)
print(f"과정이 시작된 후 10초가 지났을 때의 온도: {T_t:.2f} °F")

result>

과정이 시작된 후 10초가 지났을 때의 온도: 257.50 °F

문제는 풀었지만, 갑자기 호기심이 이를 시각화해보았다. 교재그림 5.13(a): 누진과정을 참고했다.

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['axes.unicode_minus'] = False
plt.rc('font', family='Malgun Gothic')

# 초기 조건 및 상수
T0 = 75  # 초기 온도 (°F)
T_inf = 300  # 최종 온도 (°F)
tau = 6  # 시정수 (초)
t_end = 20  # 관찰 시간의 끝 (초)

# 시간 배열 생성
t = np.linspace(0, t_end, 1000)

# 1차 시스템의 응답 계산
T_t = T_inf + (T0 - T_inf) * np.exp(-t / tau)

# 시각화
plt.figure(figsize=(10, 6))
plt.plot(t, T_t, label='온도 변화', color='orange')
plt.axvline(x=10, color='r', linestyle='--', label='t = 10초')
plt.title('시간에 따른 온도 변화 (300°F에서 75°F로 계단형 변화)')
plt.xlabel('시간 (초)')
plt.ylabel('온도 (°F)')
plt.grid(True)
plt.legend()
plt.show()

result>

5장 예제 3번

(하,,, 왜이렇게 할게 많냐,,, 다른 일처리하다 보니 벌써 새벽3시 7분이네… 이것만 다 풀고 빨리 자야겠다.ㅜㅜ)

이번에도 아까워 똑같은 수식을 가져와서 접근하면 된다. 일단 문제부터 풀어보자.

import numpy as np

# 초기조건:공식풀이의 온도표기를 고려해서 나도 P를 쓰기로 했다.
P_inf = 75
P_A = 300
t = 10
tau = 6

# 1차 시스템의 응답 계산
P = 75 + (P_A - P_inf) * np.exp(-t/tau)
print(f"과정이 시작된 후 10초가 지났을 때의 온도: {P:.2f} °F")

result>

과정이 시작된 후 10초가 지났을 때의 온도: 117.50 °F

위 코드를 시각화해보니 다음과 같더라. 교재그림 5.13(b): 감소과정을 참고했다.

import numpy as np
import matplotlib.pyplot as plt

# 초기 조건 및 상수
P_A = 300  # 초기 온도 (°F)
P_inf = 75  # 최종 온도 (°F)
tau = 6  # 시정수 (초)
t_target = 10  # 관찰 시간의 끝 (초)

# 시간 배열 생성
t = np.linspace(0, t_target*2, 1000)

# 1차 시스템의 응답 계산
P = T_inf + (T0 - T_inf) * np.exp(-t / tau)

# 시각화
plt.figure(figsize=(10, 6))
plt.plot(t, P, label='온도 변화', color='orange')
plt.axvline(x=10, color='r', linestyle='--', label='t = 10초')
plt.title('시간에 따른 온도 변화 (300°F에서 75°F로 계단형 변화)')
plt.xlabel('시간 (초)')
plt.ylabel('온도 (°F)')
plt.grid(True)
plt.legend()
plt.show()

result>

5장 예제 4번

온도 탐침이 어떤 특별한 기체의 흐름을 잴 때 시정수는 10초이다. 기체의 온도는 75°F 와 300°F 사이에서 20초의 주기(0.05 Hz)를 가지고 조화적으로 변한다. 입력되는 온도 는 입력되는 기체 온도의 항으로 표시하면 얼마가 되나? 정확한 온도의 99%정도 되는 값이 읽혀지려면 시정수는 얼마가 되어야 하나?

풀이는 다음과 같다.(GPT형님의 도움을 조금 받았다ㅋㅋ)

파이썬코드로 구현하면 다음과 같다.

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['axes.unicode_minus'] = False
plt.rc('font', family='Malgun Gothic')

# 초기 조건 및 상수
T_min = 75  # 최소 온도 (°F)
T_max = 300  # 최대 온도 (°F)
tau = 10  # 시정수 (초)
T0 = 187.5  # 온도 변화의 평균값 (°F)
A = 112.5  # 온도 변화의 진폭 (°F)
T_period = 20  # 온도 변화의 주기 (초)
omega = 2 * np.pi / T_period  # 각주파수 (rad/s)

# 시간 배열 생성
t = np.linspace(0, 60, 1000)  # 60초 동안의 시간

# 기체 온도 변화
T_in = T0 + A * np.cos(omega * t)

# 1차 시스템의 응답 계산 (주기적 입력 신호에 대한)
T_out = T0 + A / np.sqrt(1 + (omega * tau)**2) * np.cos(omega * t - np.arctan(omega * tau))

# 시각화
plt.figure(figsize=(10, 6))
plt.plot(t, T_in, label='입력 기체 온도', color='blue')
plt.plot(t, T_out, label='출력 온도 탐침 응답', color='orange')
plt.title('기체 온도 변화와 온도 탐침의 응답')
plt.xlabel('시간 (초)')
plt.ylabel('온도 (°F)')
plt.grid(True)
plt.legend(loc="upper right")
plt.show()

result>

이거랑 학생연구원 교육 이수하느라, 그리고 소그룹 활동 비용 처리 문서 진행 하느라 시간을 많이 허비했다ㅜㅠ….

오늘은 늦었으니, 내일 나머지 문제 다 풀자!