[참조]

공학자를 위한 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()