メディアネットワーク演習II（サウンドプログラミング）

科目名：　メディアネットワーク演習II（2021年～）
対象：　メディアネットワークコース3年目
日時：　4月19日（月）8:45～10:15
場所：　オンライン
レポート提出締切：　4月26日（月）13:00
レポート提出先：　メールにファイルを添付し、aoki@ime.ist.hokudai.ac.jpまで提出すること．
連絡先：　青木 直史（情報エレクトロニクス棟6階6-07）（Tel: 011-706-6532）（E-mail: aoki@ime.ist.hokudai.ac.jp）

目的

　音はマルチメディアコンテンツを構成する重要な要素である．本演習は，Pythonによる具体例を通して，サウンドプログラミングに対する理解を深めることを目的としている．

1．はじめに

　本演習は，Jupyter Notebookを利用し，ブラウザを使ってプログラムを実行しながら進めるものとする．なお，音を確認する場合は各自のイヤフォンまたはヘッドフォンを使うこと．環境のインストールはつぎの手順のとおり．

(1) Anacondaのインストール

https://www.anaconda.com/distribution/
から「Python 3.8」のインストーラをダウンロードして実行．

(2) Jupyter Notebookの起動

Anaconda Promptを開く．

jupyter notebook
と入力して実行．

(3) ブラウザ（Chrome推奨）を利用してプログラミングを行う．

「New」ボタンをクリックして「Python 3」を選択．

音ファイルや画像ファイルを利用するときは「Upload」ボタンをクリックし，必要なファイルをJupyter Notebookにアップロードすること．

2．音をつくる

　「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．サイン波の音を耳で聞いて確かめなさい．

(1)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

fs = 8000
length_of_s = int(fs * 1)

s = np.zeros(length_of_s)
for n in range(length_of_s):
    s[n] = 0.5 * np.sin(2 * np.pi * 1000 * n / fs)

(3)

for n in range(length_of_s):
    s[n] = (s[n] + 1.0) / 2.0 * 65536.0
    if s[n] > 65535.0:
        s[n] = 65535.0
    elif s[n] < 0.0:
        s[n] = 0.0;
    s[n] = (s[n] + 0.5) - 32768

wavfile.write('sine.wav', fs, s.astype(np.int16))

(4)

Audio('sine.wav')

3．音をながめる

　「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．サイン波の音を目で見て確かめなさい．

(1)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def Hanning_window(N):
    w = np.zeros(N)
    if N % 2 == 0:
        for n in range(N):
            w[n] = 0.5 - 0.5 * np.cos(2 * np.pi * n / N)
    else:
        for n in range(N):
            w[n] = 0.5 - 0.5 * np.cos(2 * np.pi * (n + 0.5) / N)
    return w

(3)

np.random.seed(0)

fs, s = wavfile.read('sine.wav')
s = s.astype(np.float)
length_of_s = len(s)
for n in range(length_of_s):
    s[n] = s[n] + 32768.0 + ((np.random.rand() * 2.0) - 1.0)
    s[n] = s[n] / 65536.0 * 2.0 - 1.0

(4)

t = np.zeros(length_of_s)
for n in range(length_of_s):
    t[n] = n / fs

(5)

plt.figure()
plt.plot(t, s)
plt.xlim([0, 0.01])
plt.ylim([-1, 1])
plt.xlabel("time [s]")
plt.ylabel("amplitude")
plt.savefig('waveform.png')

(6)

N = 1024

w = Hanning_window(N)
x = np.zeros(N)
for n in range(N):
    x[n] = s[n] * w[n]

(7)

X = np.fft.fft(x)
X_abs = np.abs(X)

(8)

f = np.zeros(N)
for k in range(N):
    f[k] = fs * k / N

(9)

plt.figure()
plt.plot(f[0 : int(N / 2)], X_abs[0 : int(N / 2)])
plt.xlim([0, int(fs / 2)])
plt.ylim([0, 200])
plt.xlabel("frequency [Hz]")
plt.ylabel("amplitude")
plt.savefig('frequency_characteristics.png')

4．スペクトログラム

　スペクトログラムは，窓関数を少しずつずらしながら周波数分析を行い，横軸を時間，縦軸を周波数として，周波数特性を濃淡表示した2次元の画像である．「New」ボタンをクリックし，新しくウィンドウを作成しなさい．つづいて，以下のプログラムを順番にセルに貼りつけ，実行しなさい．サイン波の音のスペクトログラムを表示しなさい．

(1)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from IPython.display import display, Audio

(2)

def Hanning_window(N):
    w = np.zeros(N)
    if N % 2 == 0:
        for n in range(N):
            w[n] = 0.5 - 0.5 * np.cos(2 * np.pi * n / N)
    else:
        for n in range(N):
            w[n] = 0.5 - 0.5 * np.cos(2 * np.pi * (n + 0.5) / N)
    return w

(3)

np.random.seed(0)

fs, s = wavfile.read('sine.wav')
s = s.astype(np.float)
length_of_s = len(s)
for n in range(length_of_s):
    s[n] = s[n] + 32768.0 + ((np.random.rand() * 2.0) - 1.0)
    s[n] = s[n] / 65536.0 * 2.0 - 1.0

(4)

N = 512
shift_size = 64
number_of_frame = int((length_of_s - N) / shift_size)

(5)

x = np.zeros(N)
w = Hanning_window(N)
S = np.zeros((int(N / 2 + 1), number_of_frame))

(6)

for frame in range(number_of_frame):
    offset = shift_size * frame
    for n in range(N):
        x[n] = s[offset + n] * w[n]

    X = np.fft.fft(x)
    X_abs = np.abs(X)

    for k in range(int(N / 2 + 1)):
        S[k, frame] = 20 * np.log10(X_abs[k]);

(7)

plt.figure()
xmin = (N / 2) / fs
xmax = (shift_size * (number_of_frame - 1) + N / 2)/ fs
ymin = 0
ymax = fs / 2
plt.imshow(S, aspect = 'auto', cmap = 'Greys', origin='lower', vmin = 0, vmax = 20, extent = [xmin, xmax, ymin, ymax])
plt.xlim([0, length_of_s / fs])
plt.xlabel("time [s]")
plt.ylabel("frequency [Hz]")
plt.savefig('spectrogram.png')

5．レポートについて

　下記の課題について，レポートを作成しなさい．

（課題1） n倍音の振幅が1/nで定義される波形をノコギリ波と呼ぶ．ノコギリ波をつくり，音を確認したのち，波形，周波数特性，スペクトログラムを表示しなさい．これらを見比べ，音の高さを割り出すには，どのような特徴に着目すればよいか考察しなさい．

（課題2） 基本音の初期値を500Hzとし，1秒ごとに音の高さが500Hzずつ階段状に高くなっていくノコギリ波をつくりなさい．音を確認したのち，スペクトログラムを表示し，音を正しくつくることができているかどうか確認しなさい．標本化周波数が8kHzの場合，音の高さは最大でいくつになるか考察しなさい．

Last Modified: April 1 12:00 JST 2021 by Naofumi Aoki
E-mail: aoki@ime.ist.hokudai.ac.jp