01.短時間フーリエ解析

1. イントロダクション

   音声解析には歴史があり、いろいろな手法が既に知られています。 ここでは、まず最初にフーリエ変換を利用して解析を行います。

   フーリエ変換では、信号が無限にくり返していることが仮定されています。 しかし、実際の音楽や音声信号では、音は一定で鳴り続けるのではなく、時間が経つと変化してしまいます。 こういった信号に対して、フーリエ解析をかける手法として、短時間フーリエ変換があります。 この手法は、信号の一部分を窓関数で切り出してきて、その一部の信号が無限にくり返しているものとして フーリエ解析をかけるといったものです。 こうすることで、信号に含まれる周波数成分の時系列変化が解析できます。

2. pythonを使った音声解析

   それでは、実際に短時間フーリエ変換を用いて解析を試してみます。 ここではpythonを用いることとします。

   python3にnumpy, matplotlib, soundfileを入れて、次のコードを試します。

   
   #!/usr/bin/env python
   
   import numpy as np
   import matplotlib.pyplot as plt
   from matplotlib.colors import LogNorm
   import soundfile as sf
   
   NFFT = 4096
   overlap = int(NFFT * 3 / 4)
   
   window = np.hamming(NFFT)
   
   def readwav(filename):
       spec = []
       with sf.SoundFile(filename, 'r') as f:
           print('sound file %d samples %f sec'%(f.frames, f.frames/f.samplerate))
           time_x = []
           fftfreq = np.fft.rfftfreq(NFFT, 1/f.samplerate)
           for pos in range(0, f.frames, NFFT-overlap):
               f.seek(pos)
               data = f.read(NFFT)
               if data.shape[0] %lt; NFFT:
                   break
               if data.ndim > 1:
                   windowed = window * data[:,1]
               else:
                   windowed = window * data
   
               fft_result = np.fft.rfft(windowed)
               fft_data = np.abs(fft_result)
   
               spec.append(fft_data)
               time_x.append(pos/f.samplerate)
   
       spec = np.array(spec).T
       time_x = np.array(time_x)
       X,Y = np.meshgrid(time_x, fftfreq)
   
       plt.pcolormesh(X, Y, spec, norm=LogNorm())
       plt.xlabel("time[s]")
       plt.ylabel("frequency[Hz]")
       plt.yscale('log')
       plt.ylim((20, 20000))
       plt.colorbar(orientation='horizontal')
       plt.show()
   
   if __name__ == '__main__':
       import sys
   
       argv = sys.argv
       argc = len(argv)
   
       if argc < 2:
           print('Usage: python %s input.wav'%argv[0])
           quit()
   
       readwav(argv[1])
   

3. 解析の例

   • 複数のsin波が混合している波形

     

     300Hz, 440Hz, 880Hzを混合しています。

     
     #!/usr/bin/env python
       
     import soundfile as sf
     import numpy as np
     import matplotlib.pyplot as plt
     
     sin = [(440, 1, 0),(880, 0.5, 0.2),(300, 0.3, 0.5)]
     length = 2.0
     pre = 0.5
     post = 0.5
     samplefq = 44100
     
     
     t = np.linspace(0, length, num=length*samplefq)
     buf = np.zeros(t.shape[0])
     
     for fq, amp, phase in sin:
         buf += amp*np.sin(t*fq*2*np.pi+phase*2*np.pi)
     
     wav = np.concatenate([np.zeros(int(pre*samplefq)), buf, np.zeros(int(post*samplefq))], axis=0)
     
     #plt.plot(wav)
     #plt.show()
     
     wav /= wav.max()
     wav *= 0x4000
     wav = wav.astype(np.int16)
     
     with sf.SoundFile('sin2.wav', 'w', samplefq, 1) as sfile:
         sfile.write(wav)
     

   • スイープ信号

     

     110Hzから1760Hzまでスイープする信号です。

     
     #!/usr/bin/env python
       
     import soundfile as sf
     import numpy as np
     import matplotlib.pyplot as plt
     from scipy.signal import chirp
     
     length = 5.0
     pre = 0.5
     post = 0.5
     samplefq = 44100
     
     swstart = 110
     swend = 1760
     amp = 1
     
     const=1 
               
     t = np.linspace(0, length, num=length*samplefq)
     buf = np.zeros(t.shape[0])
              
     buf += amp*chirp(t, swstart, length, swend)
     
     wav = np.concatenate([np.zeros(int(pre*samplefq)), buf, np.zeros(int(post*samplefq))], axis=0)
     
     #plt.plot(wav)
     #plt.show()
         
     wav /= wav.max()
     wav *= 0x7000
     wav = wav.astype(np.int16) 
         
     with sf.SoundFile('sweep.wav', 'w', samplefq, 1) as sfile:
         sfile.write(wav)
     

   • 楽曲

     

     ようこそジャパリパークへ(どうぶつビスケッツ×PPP) 作詞/作曲:大石昌良 (amazon link) の冒頭部分です。

4. 短時間フーリエ変換の特性

   短時間フーリエ変換を用いた解析では、解析に用いるサンプル数(ここではNFFT=4096)と窓関数(ここではハミング窓) によって得られるスペクトラムの見え方が変わります。

   • 窓関数を変えることにより、メインローブの大きさとサイドローブの形が変わります。 できるだけメインローブの幅が細くて、サイドローブが小さい窓関数が望ましいのですが、相反する条件のため解析に応じた窓関数を選ぶ必要があります。

   • 解析に用いるサンプル数を変えると、時間分解能および周波数分解能が変わります。 Nを小さくすると、時間分解能が上がる代わりに周波数分解能が落ちます。特に、Nサンプル以下の周波数をもつ低い波を検出することができなくなります。 逆にNを大きくすると、周波数分解能が上がる代わりに時間分解能が落ちます。特に、Nサンプル内で周波数が変化したとしてもその変化を捉えることができません。

   NFFT=4096という大きめの解析幅を用いているので、混合sin波やスイープ信号はうまく捉えられています。 しかしながら、楽曲の解析においては低音領域においてまだ分離が十分ではなく、4096/44100=90msの時間分解能もまた十分とは言えず、ぼんやりとした結果が得られています。

   このことを解決しようとするのが次の手法です。