山内セミナーⅠ(2021/07/14)

関連サイトと資料

サンプルプログラム(11)

list36.py
import numpy as np
import matplotlib.pyplot as plt
  
if __name__ == '__main__':
  fs = 128                    # number of samples
  amp0, amp1 = 0.8, 0.15      # amplitude
  freq0, freq1 = 2.0, 15.0    # freq
  
  t = np.linspace(0, 2*np.pi, fs)
  base = amp0*np.sin(freq0*t) # 基本波
  nois = amp1*np.sin(freq1*t) # ノイズ
  
  plt.plot(t, base)           # 基本波 表示
  plt.plot(t, nois)           # ノイズ 表示
  plt.show()
  
  signal = base + nois        # 基本波 + ノイズ
  plt.plot(t, signal)         # 基本波 + ノイズ 表示
  plt.show()
  
  f = np.fft.fft(signal)      # フーリエ変換(FFT)
  f_abs = np.abs(f)
  plt.plot(f_abs[:int(fs/2)+1])
  plt.show()
  
  fc = 14                     # cut off freq.
  f[fc:]=0                    # low pass filter
  f_abs = np.abs(f)
  plt.plot(f_abs[:int(fs/2)+1])
  plt.show()
  
  f_ifft = np.fft.ifft(f*2)   # 逆フーリエ変換(IFFT)
  plt.plot(t, f_ifft)
  plt.show()
    

list37.py
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
  
if __name__ == '__main__':
  nframes = 1024          # number of samples
  dt = 0.001              # サンプリング周期 [s]
  freq0, freq1 = 3, 50    # freq [Hz]
  amp0, amp1 = 2.0, 0.8   # amplitude
  
  t = np.arange(0, nframes*dt, dt)  # 時間 [s]
  base = amp0*np.sin(2*np.pi*freq0*t)
  nois = amp1*np.sin(2*np.pi*freq1*t)
  
  ylim=max([amp0, amp1])+1
  plt.ylim([-ylim, ylim])
  plt.plot(t, base)       # 基本波 表示
  plt.plot(t, nois)       # ノイズ 表示
  plt.show()
  
  in_signal = base + nois
  plt.plot(t, in_signal)  # 入力波形 表示
  plt.show()
  
  taps = 63               # number of filter taps
  fc = 30                 # cut off freq.
  k = signal.firwin(numtaps = taps, cutoff = fc, fs = 1/dt)
  
  # flat fir
  out_signal = np.zeros(nframes)
  for i in range(in_signal.size):
    out_signal[i] = 0
    for j in range(k.size):
      if (i-j) >= 0:
        out_signal[i] += k[j]*in_signal[i-j]
  
  plt.ylim([-ylim, ylim])
  plt.plot(t, out_signal)
  plt.show()
    

list38.py
import sys
import os
import wave
import numpy as np
from scipy import signal
  
if __name__ == '__main__':
  target = 'sample1.wav'
  bn, ext = os.path.splitext(target)
  savename = bn + '_lowpass' + ext
  
  #音声ファイル読み込み
  in_wav = wave.open(target, 'r')

  p = in_wav.getparams()
  if p.nchannels != 1:
    print('error: input must monaural!')
    exit()
  
  data = in_wav.readframes(p.nframes)
  tmp_data = np.frombuffer(data, dtype='int16')
  in_signal = tmp_data.copy().astype(dtype='float32')
  
  taps = 63               # number of filter taps
  fc = 5000.0             # cut off freq.
  nyq = p.framerate / 2.0 # ナイキスト周波数
  fc_nyq = fc / nyq       # ナイキスト周波数が1になるように正規化
  k = signal.firwin(numtaps = taps, cutoff = fc_nyq)
  
  out_signal = signal.lfilter(k, 1, in_signal)
  
  #出力ファイル
  out_wav = wave.Wave_write(savename)
  out_wav.setparams(p)
  
  x = out_signal.astype(dtype='int16')
  
  #書き込み
  out_wav.writeframes(x)
  
  in_wav.close()
  out_wav.close()
    

サンプルプログラム(12)

list39.py
import wave
import numpy as np
import matplotlib.pyplot as plt
  
if __name__ == '__main__':
  target = 'sample1.wav'
  
  #音声ファイル読み込み
  in_wav = wave.open(target, 'r')
  
  p = in_wav.getparams()
  tdata = in_wav.readframes(p.nframes)
  data = np.frombuffer(tdata, dtype='int16')
  
  #表示
  plt.specgram(data, Fs=p.framerate)
  
  plt.xlabel('time [sec]')
  plt.ylabel('frequency [Hz]')
  plt.show()
  
  in_wav.close()
    

list40.py
import wave
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
  
if __name__ == '__main__':
  target = 'sample1.wav'
  
  #音声ファイル読み込み
  in_wav = wave.open(target, 'r')
  
  p = in_wav.getparams()
  tdata = in_wav.readframes(p.nframes)
  data = np.frombuffer(tdata, dtype='int16')
  
  #表示
  cm = mpl.cm.binary
  plt.specgram(data, Fs=p.framerate, cmap=cm, NFFT=256, mode='psd');
  
  plt.xlabel('time [sec]')
  plt.ylabel('frequency [Hz]')
  plt.show()
  
  in_wav.close()