![[Python] FFT 돌리기](https://img1.daumcdn.net/thumb/R750x0/?scode=mtistory2&fname=https%3A%2F%2Fblog.kakaocdn.net%2Fdn%2F9nFvz%2Fbtq9Lq9uZcv%2FZmWBqjXL7RqrZTXfzjDYXK%2Fimg.png)
0. 관련 글 리스트
MATLAB에서도 FFT를 돌렸으면 Python에서도 FFT를 돌릴 수 있다.
어디서 본 말이 있는데, SciPy와 NumPy만 잘 다뤄도 MATLAB을 완벽대체 할 수 있다고 한다.
음, 내 생각엔 글쎄 .. 완벽대체까지는 아닌 듯.
FFT를 돌리는 flow는 내가 이전에 썼던 MATLAB에서 FFT 돌리기 글을 그대로 따라갈 것이다.
먼저, 그냥 4Hz sinewave signal 하나 만들어서 FFT를 돌려보았다.
Python 에서는 SciPy 혹은 NumPy를 이용해 FFT를 돌릴 수 있는데, 여기서는 NumPy를 이용했다.
<<코드가 보이지 않는다면 오른쪽 아래의 달모양 버튼으로
테마를 다크모드에서 화이트모드로 바꿔주세요>>
# Import Packages
import numpy as np
import matplotlib.pyplot as plt
fin = 4 # Frequency of the Sinewave
Fs = 100 # Sampling Rate
Ts = 1/Fs # Sampling Interval
Nsamp = np.power(2, 13) # Number of Samples (= 2^13 = 8,192)
t = np.linspace(0, Nsamp*Ts, Nsamp) # Time-Axis
y = 2*np.sin(2*np.pi*fin*t) # Sinewave
FFT_data = np.fft.fft(y) # Calculate FFT
# Plot the FFT
fig, ax = plt.subplots(1, 2, figsize = (20,5))
ax[0].plot(t, y)
ax[0].set_xlim(0,1)
ax[1].stem(abs(FFT_data))
plt.show()
역시 MATLAB에서 봤던 현상과 비슷한 게 나타난다.
time-axis에 관련한 정보를 fft 함수에 안 넣었으므로 그건 그렇다 쳐도
amplitude 정보를 보면 역시 들어맞지가 않는다.
따라서, MATLAB에서 했던 방법 그대로 똑같이 PositiveFFT를 만들어보기로 했다.
1. Positive FFT 함수
import numpy as np
def PositiveFFT(x, Fs):
Ts = 1/Fs
Nsamp = x.size
xFreq = np.fft.rfftfreq(Nsamp, Ts)[:-1]
yFFT = (np.fft.rfft(x)/Nsamp)[:-1]*2
return xFreq, yFFT
이러저리 구글을 뒤적여가며 만들어보니까 MATLAB에서 하던 것보다
훨씬 더 간단히 몇 줄의 line들 만으로 구현이 가능했다.
첫째로, MATLAB에서 X-Axis를 뽑아내기위해 했던 짓들을
한 방에 처리할 수 있는 fftfreq라는 함수가 NumPy에 있기 때문이다.
둘째로, rfft 라는게 NumPy에 내장되어 있는데,
negative frequency를 날려주는 함수이다.
즉, ceil을 이용해가며 이런저런 짓을 할 필요가 없었다.
단, negative frequency를 날려주기만 하는 함수이므로
잃어버린 정보를 감안해 x2를 해주어야 한다.
또한, 아래 링크에도 나와있듯이 $N_{samp}$가 짝수 (even)개라면
positive frequency의 마지막 element를 살려놓는데,
사실상 ±$f_{s}/2$, 즉 음과 양의 Nyquist Frequency 지점이 겹쳐진 곳이라
의미있는 data가 아니다.
따라서 이 data는 제외시키기 위해 [:-1]을 이용해
array의 마지막 element를 날려주었다.
numpy.fft.rfft — NumPy v1.21 Manual
Normalization mode (see numpy.fft). Default is “backward”. Indicates which direction of the forward/backward pair of transforms is scaled and with what normalization factor. New in version 1.20.0: The “backward”, “forward” values were added.
numpy.org
위 함수를 .py로 저장하고 같은 폴더에 둔 뒤,
main code에 import하는 방식으로 사용할 수 있다.
이게 귀찮다면 그냥 한 main code에 함수든 뭐든 다 몰아써도 된다.
[PyTorch] 파이토치로 딥러닝 초보자 입문하기 - 설치하기
0. 관련 글 리스트 [PyTorch] 파이토치로 딥러닝 - MLP Network으로 FashionMNIST Dataset Classification 해보기 파이썬(Python)과 딥러닝에 대해 아는게 거의 없다시피한 사람이 PyTorch에 입문하기로 했다 ! 노..
tuttozurich.tistory.com
2. FFT 돌리기
...
xFreq, FFT_data = PositiveFFT(y, Fs) # Calculate FFT
# Plot the FFT
fig, ax = plt.subplots(1, 3)
ax[0].plot(t, y)
ax[0].set_xlim(0,1)
ax[1].stem(xFreq, abs(FFT_data))
ax[2].semilogx(xFreq, 20*np.log10(abs(FFT_data)) )
plt.show()
X-Axis를 보면 frequency가 4Hz로 잘 나오고 있음을 알 수 있으나,
Amplitude를 보면 MATLAB에서 보던 현상과 같은 현상이 나온다.
Peak amplitude도 1.75 근처이고, dB-scale로 보면
skirt도 보이고 noise floor도 -50dB 근처이다.
3. Window 씌워서 FFT 돌리기
import numpy as np
def PositiveFFT(x, Fs, w=None):
Ts = 1/Fs
Nsamp = x.size
xFreq = np.fft.rfftfreq(Nsamp, Ts)[:-1]
if w is None:
yFFT = (np.fft.rfft(x)/Nsamp)[:-1]*2
else:
yFFT = (np.fft.rfft(x*w)/Nsamp)[:-1]*2
return xFreq, yFFT
Window도 input argument로 받기 위해
위와 같이 PositiveFFT 함수를 수정했다.
w=None으로 두면 default가 None이며,
PositiveFFT 함수를 사용할 때 w칸을 비워놓고 쓰면
windowing을 추가로 하지 않게 만들었다.
w가 None이면 x에 w를 곱하지 않고,
None이 아니면 따로 window function을 지정해서 넣은 경우일테니
x와 w를 곱해서 FIR filtering을 적용시켜준다.
아래는 PositiveFFT를 불러와 사용하는 main code이다.
...
w = np.hanning(Nsamp) # Hann Window
xFreq1, FFT_data1 = PositiveFFT(y, Fs) # Calculate FFT without Windowing
xFreq2, FFT_data2 = PositiveFFT(y, Fs, w) # Calculate FFT with Windowing
# Plot the FFT
fig, ax = plt.subplots(1, 3, figsize = (10,5))
ax[0].plot(t, y)
ax[0].set_xlim(0,1)
ax[1].stem(xFreq1, abs(FFT_data1))
ax[1].stem(xFreq2, abs(FFT_data2), 'orange')
ax[1].set_xlim(3.5, 4.5)
ax[1].set_ylim(0, 1)
ax[2].semilogx(xFreq1, 20*np.log10( abs(FFT_data1)) )
ax[2].semilogx(xFreq2, 20*np.log10( abs(FFT_data2)), 'orange')
plt.show()
MATLAB에서 보던 현상과 같은 현상이 나온다.
dB-scale로 보면 skirt가 좀 더 가라앉았으며
noise floor도 -150dB 근처이다.
Linear-scale graph는 잘 안 보이므로
이 친구만 따로 떼내어 놓고 보면 아래와 같다.
plt.stem(xFreq2, abs(FFT_data2), 'orange')
plt.xlim(3.5, 4.5)
plt.ylim(0, 1)
plt.show()
4Hz input에 대해, 4Hz만 보이는게 아니라 주변 값들이 추가로 보이는 것을 알 수 있다.
(MATLAB과는 달리 MatPlotLib은 마우스로 찍어서 임의로 marker를 놓는게 안 된다.
그래서 그냥 4Hz에 인접한 point 3개를 뽑아내는 code를 간단하게 짜보았다.
이 code는 SNDR 계산하는데 사용하면 유용할 듯)
def Find_Nearest(x, value):
idx = (abs(x - value)).argmin()
return idx
idx = Find_Nearest(xFreq, fin)
mag_sum = 0
for i in range(-1,2):
mag = abs(FFT_data)[idx+i]
print('Freq = ', xFreq[idx+i], '/ Magnitude = ', mag)
mag_sum = mag_sum + mag
print(mag_sum)
Find_Nearest 함수는 x라는 array에서 value에 제일 가까운 값의 index를 뽑아내준다.
이걸 이용하면 xFreq array에서 4Hz랑 제일 가까운 index를 뽑아낼 수 있고,
FFT magnitude 중 idx-1번째, idx번째, idx+1번째 값들을 모아 다 더해주는 code이다.
그 결과, 아래와 같은 print 결과를 얻었다.
>>Freq = 3.99169921875 / Magnitude = 0.7072846134339216
>>Freq = 4.00390625 / Magnitude = 0.9503485328902309
>>Freq = 4.01611328125 / Magnitude = 0.3002118064935845
>>1.9578449528177368
합은 1.95로, MATLAB에서 보던 값이랑 같은 값이다.
즉, 이 code가 잘 되고 있다는 것을 의미한다.
4. Coherent Sampling을 통해서 FFT 돌리기
대강의 설명은 MATLAB 글에서 했으니,
여기서는 직접 coherent sampling parameter를 뽑아내는
calculator code를 짜보았다.
def Find_Nearest_Prime_Number(x):
x = np.int64(np.ceil(x))
x_list = []
for num in range(0, x*2):
if num > 1:
for i in range(2, num):
if (num % i) == 0:
break
else:
x_list.append(num)
idx = Find_Nearest(x_list, x)
return x_list[idx]
일단, 어떤 값 (x)이랑 제일 가까운 소수 (Prime Number)를 찾아내는 code를 짰다.
$f_{in}$ 및 $F_{s}$, $N_{samp}$를 이용해 $Mcycle$을 계산해보면 integer로 나오질 않으므로
우선 이걸 ceil로 integer로 자른 뒤 datatype을 강제로 int64로 만들었다.
ceil을 적용하고 나면 datatype이 float가 되기 때문이다.
그 후, 0과 x의 2배로 널찍하게 range를 잡아놓고 소수를 모두 모아다가 list를 만들어 return 한다.
def Coherent_Sampling_Calculator(fin, Fs, Nsamp):
Mcycle = Find_Nearest_Prime_Number(Nsamp * fin / Fs)
fin_new = Fs * Mcycle / Nsamp
Fs_new = fin * Nsamp / Mcycle
return fin_new, Fs_new
그리고 위와 같이 $M_{cycle}$을 계산해서 Prime Number를 구하고,
$f_{in}$과 $F_{s}$를 tuple로 return 한다.
FFT를 돌릴 때는 둘 중 아무거나 하나 잡아다가 적용시켜 돌리면 된다.
print(Coherent_Sampling_Calculator(fin, Fs, Nsamp))
>>(4.04052734375, 98.99697885196375)
위 code를 돌려보면 coherent sampling을 만족하는 $f_{in_new}$ 및 $F_{s_new}$를 return 해준다.
$F_{s_new}$로 반환된 98.99.. 를 복사 붙여넣기 해서 FFT를 다시 돌리니 아래와 같이 나온다.
3개의 amplitude 쌍을 모두 더한 결과는 아래와 같이 출력되었다.
1.999..로, 2에 매우 가까워졌고 이는 MATLAB에서 보던 결과랑 같다는 것을 알 수 있다.
Freq = 3.987915407854985 / Magnitude = 0.469834671708894
Freq = 4.0 / Magnitude = 0.9988255921130113
Freq = 4.0120845921450154 / Magnitude = 0.5303971524806151
1.9990574163025205
5. Spectral Leakage를 조금 더 최소화해서 FFT 돌리기
MATLAB 글에서 언급한 바와 같다. Sample 개수를 늘려주면 된다.
단, 이번에는 Blackman Window도 적용하여 비교를 해보았다.
위 그래프는 DC component를 날려 준 후 plot한 것이다.
파란색은 Hann Window에 213개의 sample 개수,
빨간색은 Hann Window에 216개의 sample 개수,
초록색은 Blackman Window에 216개의 sample 개수를
적용하여 plot한 것이다.
Blackman Window가 sinewave를 표시하는데 성능이 좋음을 알 수 있다.
여기서 Blackman Window는 Blackman-Harris Window와는 다른 window이다.
해당 window를 사용하고 싶다면 MATLAB에서 만들어낸 후
갖고 와서 사용하는 것이 하나의 방법이 되겠다.
Python 환경이 구축되어 있지 않다면
coherent sampling calculator 정도는 인터넷에 떠도는
Online Python Compiler로 쉽게 돌려볼 수 있다.
(구글에 Python Compiler 검색)
