python〜ファミマの入店音をフーリエ変換してみる

こんにちは、はっこんです。

皆さん、フーリエ変換はご存知ですか？

知らない方も多いと思いますので、簡単に説明すると、

信号などの時間領域のデータを周波数領域へ変換するものです。

実際のフーリエ変換の式はこんな感じ、

F(\omega)=\int_{-\infty}^{\infty} f(t)e^{-i \omega t}dt

この式を用いれば、時間依存の信号を周波数に分解することができます。

今回フーリエ変換する素材として、ファミマの入店音を選びました。

ちゃらららららー、ちゃららららーってやつです。

音声のプロット

素材のデータを読み取り、波形をプロットするとこんな感じ、

コードがこんな感じ、

import sys
import scipy.io.wavfile
import numpy as np
import matplotlib.pyplot as plt



args = sys.argv
wav_filename = args[0]
rate, data = scipy.io.wavfile.read('sample.wav')


data = data / 32768 #normalization
time = np.arange(0, data.shape[0]/rate, 1/rate)  
#plot
plt.figure(1)
plt.plot(time, data[:,0])
plt.savefig('famima')

16bitの音声ファイルを正規化しました。

肝心のフーリエ変換ですが、ちょっと難しいです。

そのままフーリエ変換を試みると失敗してしまうので、その前に1つの操作が必要です。

その操作が窓関数をかけるというものです。

窓関数とは？

窓関数(ハニング関数)とは次のような関数です。

w(x)=0.54-0.46\cos(2\pi x)~~(0\leq x \leq 1)

このコードでプロットできます、

import numpy as np
import matplotlib.pyplot as plt
window = np.hamming(100)
plt.plot(window)
plt.show()

フーリエ変換は切り出した範囲が周期関数であることを仮定して行われるものですが、

適当に切り出すと両端がうまく繋がらず、周期関数にならないという問題があります。

窓関数をかけることによって両端が繋がり、周期関数としての扱いを許すことができるとのこと。

フーリエ変換

それではフーリエ変換していきましょうか！

結果はこんな感じです。

コードがこちら、

import numpy as np
import matplotlib.pyplot as plt
import wave

def ReadWavFile(FileName = "sample.wav"):

    try:
        wr = wave.open(FileName, "r")
    except FileNotFoundError:
        print("[Error 404] No such file or directory: " + FileName)
        return 0
    data = wr.readframes(wr.getnframes())
    wr.close()
    x = np.frombuffer(data, dtype="int16") / float(2**15)
    return x


sampling_rate = 44100 #sampling rate
data = np.array(ReadWavFile("sample.wav")) 

NFFT = 1024 # size of frame
OVERLAP = NFFT / 2 # half shift
frame_length = data.shape[0] 
split_number = len(np.arange((NFFT / 2), frame_length, (NFFT - OVERLAP))) 

window = np.hamming(NFFT)  

spec = np.zeros([split_number, NFFT // 4]) 

pos = 0

for fft_index in range(split_number):
    frame = data[int(pos):int(pos+NFFT)]
    if len(frame) == NFFT:
        windowed = window * frame  
        fft_result = np.fft.rfft(windowed)
        fft_data2 = np.real(fft_result[:int(len(fft_result)/2)]) 
        fft_data2 = np.log(fft_data2** 2)  

        for i in range(len(spec[fft_index])):
            spec[fft_index][-i-1] = fft_data2[i]

        pos += (NFFT - OVERLAP) # move window

# plot
plt.imshow(spec.T, extent=[0, frame_length, 0, sampling_rate/2], aspect="auto")
plt.xlabel("time[s]")
plt.ylabel("frequency[Hz]")
plt.colorbar()
plt.show()

窓関数の位置をずらしながらフーリエ変換を行っていきます。

結果のグラフを見ても確かにファミマの入店音っぽいですね(どこがww)

今日は疲れているのでここらへんにしておきます。

(たくさんのサイトを参考にさせていただきました、ありがとうございます)