Python의 FFT, 평균화(averaging), 평활화(Smooth)

---

오디오의 FFT

오디오의 주파 수축 정보를 알아보기 위해서 FFT를 사용합니다. 이글에서는 Python을 이용해서 wav 음악 파일의 일부를 읽어서 FFT의 결과를 그래픽으로 표시해 보겠습니다. 동시에 octave 조건을 이용해 보겠습니다. 

---

평활화

주파수 축의 결과가 주파수에 따라서 신호의 크기 변화가 크기 때문에 유의미한 신호의 특성을 얻기가 힘듭니다. 그래서 특정 주파수를 기준으로 주변의 신호들과 평균을 내고 이를 대표해서 표시하는 방법을 사용하는데 이를 평활화(smooth)라고 합니다. 

평균은 다른 글의 내용을 참조해 주세요. 

평활화는 주로 옥타브 단위로 주변의 신호를 평균 내서 화면에 표시합니다. 아래의 코드에서 평활화 없는 FFT 결과, 1 oct로 평활화 그리고 1/3 oct로 평활화한 결과를 보여주겠습니다. 이 중에서 청감적인 특성을 고려한다면 1/3 oct가 유리하며, 전체적인 에너지 분포 특성을 보려면 1 oct가 적합합니다. 

---

Python 코드

다음과 같이 Python 코드를 작성합니다. 이때 scipy numpy matplotlib가 Python을 위해 설치되어 있어야 합니다. 코드 내부에 설명이 포함되어 있습니다.

• 주파수 특성을 알아보기 위해 FFT를 사용할 것입니다.
• 옥타브 주파수를 구합니다.
• wav 파일을 읽고 일부 구간에 대해서만 FFT 처리를 합니다.
• 몇개의 wav 시간 구간에 대해 평균화를 합니다.
• 결과를 그래프로 표시합니다.  

 

'''
octave fft

ref:
    https://techreviewtips.blogspot.com/2017/11/05-02-fft.html
    https://wikidocs.net/92086
'''
import math
import numpy as np
import scipy
import scipy.io as sio
import scipy.io.wavfile
import matplotlib.pyplot as plt
'''
'''
global SAMPLING_RATE
# FFT의 크기 2**N으로 만들어야 한다.
global FFT_SIZE
# wav 신호를 [-1.0~1.0]으로 정규화 하기위한 최대 정수값
global MAX_INT
# 가청대역을 분할하기위한 옥타브 값 64개 밴드로 분할
global OCT_VALUE
# 가청대역의 시작 주파수 
global OCT_START_FREQ
# 가청대역의 옥타브 밴드로 분할된 수량
global BAND_SIZE
# 옥타브로 분할된 대역폭의 낮은 주파수 배열
global f_low
# 옥타브로 분할된 대역폭의 높은 주파수 배열
global f_hi
# 옥타브로 분할된 대역폭의 중심 주파수 배열
global f_cnt
# 옥타브로 분할된 대역폭의 낮은 주파수에 해당하는 FFT 배열을 위치
global fft_lo
# 옥타브로 분할된 대역폭의 높은 주파수에 해당하는 FFT 배열을 위치
global fft_hi
# 옥타브로 분할된 대역폭의 FFT 배열을 할당 수량
global fft_cnt 
# 옥타브로 분할된 대역폭의 평균 신호
global fft_oct
"""
옥타브 밴드를 만든다.
"""
def oct_band_make():
    # 사용될 전역 변수    
    global BAND_SIZE
    global OCT_VALUE
    global SAMPLING_RATE
    global FFT_SIZE
    #주파수 해상도
    df= float(SAMPLING_RATE)/float(FFT_SIZE)    
    # 대역폭의 낮은 주파수, 높은 주파수, 중심 주파수
    f_low = np.arange(BAND_SIZE,dtype='float')
    f_hi = np.arange(BAND_SIZE,dtype='float')    
    f_cnt = np.arange(BAND_SIZE,dtype='float')    
    # 대역폭에 대한 fft 데이터 위치
    fft_lo = np.arange(BAND_SIZE)    
    fft_hi = np.arange(BAND_SIZE)    
    # 대역폭 내에 fft 데이터의 수량
    fft_cnt = np.arange(BAND_SIZE)    
    # 주파수 대역 계산을 위한 상수 계산
    B = 2.**OCT_VALUE
    # 초기 주파수 조건 입력
    f_low[0] = OCT_START_FREQ
    # 옥타브 대역 주파수 계산
    for i in range(0,BAND_SIZE-1):
        f_hi[i] = f_low[i] * B
        f_low[i+1] = f_hi[i]
        f_cnt[i] = math.sqrt(f_hi[i]*f_low[i]) 
    f_hi[BAND_SIZE-1] = f_low[BAND_SIZE-1] * B
    f_cnt[BAND_SIZE-1] = math.sqrt(f_low[BAND_SIZE-1] *f_hi[BAND_SIZE-1])
    # 대역별로 포함된 fft 데이터의 수량을 계산한다.
    for i in range(0,BAND_SIZE):
        fft_lo[i] = math.ceil(f_low[i]/df)
        fft_hi[i] = math.floor(f_hi[i]/df)
        fft_cnt[i] = fft_hi[i] - fft_lo[i] + 1    
    return f_low,f_hi,f_cnt,fft_lo,fft_hi,fft_cnt        
""" 
wav 파일을 읽어온다.
"""
def wav_read(wav_file,start_pos):
    # 사용될 전역 변수    
    global SAMPLING_RATE
    global FFT_SIZE
    global MAX_INT    
    #FFT를 위해 wav 데이터의 일부에 윈도우를 사욯한다.
    window = np.hamming(FFT_SIZE)
    # wav 파일을 읽어 wav_signal에 기록하면서,
    # SAMPLING_RATE를 알아온다.
    SAMPLING_RATE, wav_signal = sio.wavfile.read(wav_file)
    # wav 데이터를 [-1.0~1.0]정규화된 실수값으로 변환한다.
    wav_normal = wav_signal.astype(np.float32)/float(MAX_INT)
    #wav의 데이터 수량를 알아온다.
    wav_data_size = len(wav_signal)
    #wav의 데이터 수량을 이용해서 wav의 재생시간을 계산한다.
    times = np.arange(wav_data_size)/float(SAMPLING_RATE)
    # 일부 구간만 FFT를 구하기 위한 샘플로 가져온다.
    end_pos = start_pos+ FFT_SIZE
    # wav 데이터에서 일부만 추출
    wav_data = wav_normal[start_pos:end_pos]
    # 추출된 데이터에 위도우을 반영
    wav_window = wav_data * window
    return wav_window, wav_normal,times,end_pos
""" 
FFT 연산
"""
def calc_fft(wav_data):
    # FFT는 nyquest 조건에 의해 필요한 주파수 범위에 
    # 두배에 해당하는 대역(double side)을 연산한다.
    fft_result = np.fft.fft(wav_data)/FFT_SIZE
    # FFT 결과에서 nyquest 조건을 고려하여
    # 절반의 주파수 영역(single side)만 가지도록 범위 설정
    single_range_size = math.trunc(FFT_SIZE/2)
    # FFT 결과의 single_side범위를 
    # 절대값을 취해서 FFT의 진폭응답으로 사용한다.
    # fft의 크기를 상대적으로 사용한다면, 아래 식을 사용해도 된다. 
    fft_mag = abs(fft_result[range(single_range_size)])
    # 만약 doble_side의 fft 결과까지 반영하고 싶으면,
    # 아래의 식을 사용한다.
    #fft_mag = 2. * abs(fft_result[range(single_range_size)])    
    return fft_mag
"""
옥타브 대역폭 별로 fft의 크기를구한다.
"""
def calc_fft_oct(fft_mag):
    # 사용될 전역 변수
    global OCT_VALUE
    global BAND_SIZE
    global fft_lo
    global fft_cnt    
    global fft_oct
    #fft를 옥타브 단위로 평균을 계산해서 기록
    fft_oct = np.arange(BAND_SIZE,dtype='float')
    # 밴수 만큼 계산
    for i in range(0,BAND_SIZE,1):      
        # 기록될 버퍼를 초기화
        fft_oct[i] =0
        # 각 밴드에 기록된 FFT 데이터를 모두 합산
        for j in range(0,fft_cnt[i],1):
            fft_oct[i]= fft_oct[i] + fft_mag[fft_lo[i]+j]
        # 합산 결과의 평균
        fft_oct[i] = fft_oct[i] / fft_cnt[i]
    return fft_oct            
""" 
Main function
"""
if __name__ == "__main__":
    #사용될 전연 변수
    global BAND_SIZE
    global f_low
    global f_hi
    global f_cnt
    global fft_lo
    global fft_hi
    global fft_cnt 
    global fft_oct 
    # 전역 변수 초기화
    WAV_FILE = 'smile10_mono.wav'
    FFT_SIZE = 2**15
    MAX_INT = 2**15
    OCT_VALUE = 1./1.
    OCT_START_FREQ = 20.
    BAND_SIZE = math.trunc(1/OCT_VALUE*10.0)
    # wav 데이터를 읽어드릴 위치
    start_pos = 0   
    # wav 파일을 읽어서 데이터를 만든다.
    wav_data,wav_normal,times,end_pos = wav_read(WAV_FILE,start_pos)
    # wav 데이터에서 fft를 계산한다.
    fft_mag = calc_fft(wav_data)    
    # oct 대역폭 정보를 만든다. wav_read()이후에 호촐되어야 한다.
    f_low,f_hi,f_cnt,fft_lo,fft_hi,fft_cnt  = oct_band_make()    
    fft_oct = calc_fft_oct(fft_mag)
    ## 추출된 wav 데이터를 시간축으로 보여준다.
    plt.subplot(3,1,1)
    plt.fill_between(times, wav_normal)
    plt.xlim(times[start_pos],times[end_pos])
    plt.ylim(-1.2, 1.2)
    plt.xlabel('time (s)')
    plt.ylabel('amplitude')
    plt.grid()
    # 추출된 wav 데이터의 FFT 결과를 보여준다.
    plt.subplot(3,1,2)
    plt.semilogx(10)
    plt.plot(20.*np.log10(fft_mag))    
    plt.xlabel('frequency (Hz)')
    plt.ylabel('magnitude(dB)')
    plt.grid()
    # 옥타브 밴드 단위로 FFT 결과를 보여준다.
    plt.subplot(3,1,3)
    plt.semilogx(10)
    f_cnt_x = np.trunc(f_cnt)
    plt.plot(f_cnt_x,20.*np.log10(fft_oct))        
    plt.xlabel('frequency (Hz)')
    plt.ylabel('magnitude(dB)')
    plt.grid()
    plt.tight_layout()
    plt.show()

위의 코드의 실행 결과입니다. 시간축과 주파수 축에서 신호를 볼 수 있습니다. 세 번째 그림은 1 oct 단위로 신호를 평활화(smooth)한 그림입니다.

평활화를 1/3 oct로 수행하면 다음과 같습니다. 사람이 청감적으로 느끼는 신호의 특성에 잘 부합하는 그림입니다.

---

평균화

평균화(averaging)은 시간에 따라 빠르게 변화한는 신호를 시간축에서 구한 샘플들의 평균을 구해서 나타내는 방법입니다. 이는 이동 평균(moving average)를 구하는 내용입니다. 이에 대한 내용은 다음의 글을 참조하시면 됩니다.

이를 FFT에 활용이 가능합니다. 시간에 따라 FFT결과를 출력하면 주파수 변화의 내용을 읽기 어렵습니다. 그래서 주파수 별로 이동평균을 계산해서 주파수별 변화의 변화 추이를 구하면 좀더 해석이 편리해 집니다. 이를 특별히 FFT의 평균화라고 합니다.

실제 스팩트럼 분석기에서는 평활화여 평균화를 선택하는 기능이 기본적으로 포함되어 있습니다.  

---

광고좀 꾹 눌러주시면 고맙겠습니다. 

---

참조 :

https://techreviewtips.blogspot.com/2017/11/05-02-fft.html
https://wikidocs.net/92086

---

위의 내용을 참조용으로만 사용해주세요. 무단 도용이나 무단 복제는 불허합니다.

기타 문의 사항은 gigasound@naver.com에 남겨 주시면 고맙겠습니다.