[03. 패키지 - SciPy] (4) 푸리에 변환: np.fft, np.ifft, 진동수영역 해석
작성자 : kim2kie
(2023-02-19)
조회수 : 12746
[참조]
공학자를 위한 Python, 조정래, 2022: 5. SciPy
https://wikidocs.net/15636
Dookie Kim, Python: From Beginning to Application, 2022
https://www.dropbox.com/s/oa86j9ap62esmtz/Python.pdf?dl=0
SciPy는 과학기술계산을 위한 라이브러리이다.
(4) 푸리에 변환: np.fft, np.ifft, SDOF EOM
1) fft vs ifft
2) discrete Fourier transform
3) 단자유도계(SDOF) 진동수영역 해석
(4) 푸리에 변환: np.fft
numpy.fft에서 DFT(discrete Fourier transform)를 지원한다.
fft(x), ifft(x), fftfreq(n), fftshift(x) 등의 함수가 있다.
1) fft vs ifft
fft:
ifft:
Ex)
import matplotlib.pyplot as plt
import numpy as np
plt.ion() # interactive mode on; plt.ioff()
x = np.array([1.0,0.8,0.1,0.2,0.5,0.1,0,0.2]) # even number(짝수)
xf = np.fft.fft(x)
(Note: 여기서 fft함수를 사용하기 위해서 x는 1D array여야 한다.
x는 (8,) array이므로, 1D array이다.
참고로 (8,1) array는 2D array이므로 fft를 사용할 수 없다.)
plt.subplot(3,1,1)
plt.grid(True)
plt.stem(x)
plt.ylabel('x')
plt.title('Even number real signal')
plt.subplot(3,1,2)
plt.grid(True)
plt.stem(np.abs(xf))
plt.ylabel('abs(X)')
plt.subplot(3,1,3)
plt.grid(True)
plt.stem(np.angle(xf))
plt.ylabel('angle(X)')
plt.tight_layout()
2) discrete Fourier transform
60 Hz 와 120 Hz의 사인 곡선이 중첩된 신호가 있다.
샘플링 주파수는 1,000 Hz이고 1,500개의 데이터가 있다.
signal:
Ex)
import matplotlib.pyplot as plt
import numpy as np
fs = 1000 # sampling frequency 1000 Hz
dt = 1/fs # sampling period
N = 1500 # length of signal
t = np.arange(0,N)*dt # time = [0, dt, ..., (N-1)*dt]
x = 0.7*np.sin(2*np.pi*60*t) + np.sin(2*np.pi*120*t)
xr = x + 2*np.random.randn(N) # random number Normal distr, N(0,2)... N(0,2*2)
plt.ion() # interactive mode on; plt.ioff()
plt.subplot(2,1,1)
plt.plot(t[0:51],x[0:51],label='x')
plt.plot(t[0:51],xr[0:51],label='xr')
plt.legend()
plt.xlabel('time')
plt.ylabel('x(t)')
plt.grid()
# Fourier spectrum
df = fs/N # df = 1/N = fmax/N
f = np.arange(0,N)*df # frq = [0, df, ..., (N-1)*df]
xf = np.fft.fft(xr)*dt
plt.subplot(2,1,2)
plt.plot(f[0:int(N/2+1)],np.abs(xf[0:int(N/2+1)]))
plt.xlabel('frequency(Hz)')
plt.ylabel('abs(xf)')
plt.grid()
plt.tight_layout()
3) 단자유도계 운동방정식(SDOF EOM) 진동수영역 해석
시스템: 질량 1, 고유주기 2 sec, 감쇠비 0.1
하중: Half-sine 함수, 1,024개의 데이터, 지속시간 40 sec
Ex)
import matplotlib.pyplot as plt
import numpy as np
# Signal properties
N=1024; T=40; dt =T/N
t = np.arange(0,N*dt,dt) # [0,dt, ..., (N-1)*dt]
# System properties
Tn = 2.;
xi = 0.1;
wn = 2*np.pi/Tn # natural (angular) frequency
m = 1
c = 2*m*wn*xi
k = wn*wn
# Excitation force
p0 = 1; td = 2;
p = np.zeros(N)
p[t]>
# Frequency response
df = 1/T
f = np.arange(0,N*df,df) #[0,df,...,(N-1)*df]
w = 2*np.pi*f
pf = np.fft.fft(p)*dt
H = 1/(-w*w*m+c*m*1j+k) # transfer function
uf = pf*H
for i in range(0,int(N/2)):
uf[N-i-1] = np.conjugate(uf[i+1])
ut = np.real(np.fft.ifft(uf)/dt)
plt.ion() # interactive mode on; plt.ioff()
plt.subplot(2,2,1)
plt.plot(t,p)
plt.title('Excitation p(t)')
plt.xlabel('time[sec]')
plt.grid()
plt.subplot(2,2,2)
plt.plot(f[0:int(N/2+1)],abs(pf[0:int(N/2+1)]))
plt.title(r'Excitation p($\omega$)')
plt.xlabel('frq [Hz]')
plt.grid()
plt.subplot(2,2,3)
plt.plot(t,ut)
plt.title('Response u(t)')
plt.xlabel('time[sec]')
plt.grid()
plt.subplot(2,2,4)
plt.plot(f[0:int(N/2+1)],np.abs(uf[0:int(N/2+1)]))
plt.title(r'Reponse u($\omega$)')
plt.xlabel('frq [Hz]')
plt.grid()
plt.tight_layout()