第4回　スペクトログラムの理解と可視化演習

この回の資料

• 📊 スライドを開く（HTML）
• 📄 スライド（PDF）
• 形式：オンライン同期回（Zoom・録画を後日公開）

---

本日の到達目標

1. フーリエ変換の直感的な意味を説明できる
2. STFTとスペクトログラムの仕組みを理解する
3. Pythonでスペクトログラムを描画・読み取れる
4. 有声音・無声音・無音の違いをスペクトログラムで識別できる
5. 発展：FFTの周波数ビン間隔 \( \Delta f = f_s/N \) を導き、窓長・FFT長・ホップ長が分析に与える影響を区別できる

進化の視点：周波数を分ける蝸牛と、スペクトログラム

蝸牛では、基底膜の場所によって応答しやすい周波数が異なる。時間とともに変化する音を周波数ごとに捉える点で、スペクトログラムとの類似を考えられる。ただしSTFTは窓を掛けた信号のフーリエ変換であり、内耳の非線形な応答や神経の活動をそのまま再現したものではない。¹

もう一つ、フォルマントには発声器官の情報が乗る。形状を保って声道を長くすると、共鳴周波数は低くなる。アカゲザルでは、フォルマントの平均間隔（formant dispersion）と声道長・体サイズの相関が報告されている（Fitch 1997, JASA）。この種・標本での結果を、そのままヒトの体格推定に適用できるわけではない。²

技術への接続：フォルマントは声道長だけでなく、舌・顎・唇の配置でも変わる。音声認識では、話者間の違いの一部を声道長正規化などで補正する。一方、話者識別ではその違いが手がかりになる。どの変動を残すかは、目的によって変わる。

問い：同じ量が、ある分野では信号で、別の分野では雑音になる。この授業で他にそういう量はあるか。

---

1. 音とは何か？—「音はサイン波の重ね合わせ」

直感的理解から始めよう

フーリエ解析では、複雑な波形を振幅・周波数・位相の異なる正弦波や余弦波の成分で表します。有限長の標本列はDFTで分解・再構成できます。連続信号の表現には、周期性や収束についての条件があります。

例えばピアノの「ラ（440Hz）」は：

• 基音：およそ440Hzの成分（最も強い成分とは限らない）
• 倍音：およそ880Hzの成分
• 倍音：およそ1320Hzの成分
• …（高次の倍音が続く）

各成分の振幅・位相と時間変化が音色に関わります。実際のピアノには、弦の性質による倍音の整数倍からのずれや打撃音もあります。次の式は有限個の成分による模式的な表現です。

アナロジー：光とプリズム

白色光をプリズムに通すと、赤・橙・黄・緑・青・藍・紫に分かれる。 音をフーリエ変換すると、音を構成する周波数（音程）の成分に分かれる。

\[ x(t) = A_0 + \sum_{\ell=1}^{L} A_\ell \sin(2\pi f_\ell t + \phi_\ell) \]

---

2. DFT（離散フーリエ変換）の計算方法

連続 → 離散へ

コンピュータは連続信号を扱えないため、一定間隔でサンプリングされたデジタル信号（離散信号）を処理します。

DFTの定義式

\[ X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j \frac{2\pi}{N} k n}, \quad k = 0, 1, \ldots, N-1 \]

• \(x[n]\)：時間領域の離散信号（\(N\)個のサンプル）
• \(X[k]\)：周波数領域の複素数（\(k\)番目の周波数成分の大きさと位相）
• \(|X[k]|\)：振幅スペクトル（この定義では未正規化。波形の振幅と直接同じ数値ではない）

直感的な意味

\(X[k]\) は「信号の中に \(k\) 番目の周波数成分がどれだけ含まれているか」を測る。

非負の周波数側での対応は \(f_k = k f_s/N\)。実数信号では負の周波数側と共役対称になるため、片側だけを表示できる。np.fft.fft の後半は負の周波数に対応し、その対応は np.fft.fftfreq で得られる。偶数長ではナイキスト周波数のビンも特別に扱う。³

導出①：DFT は「その周波数の波と、どれだけ似ているか」を測っている

オイラーの式 \( e^{-j\theta} = \cos\theta - j\sin\theta \) を定義式に入れると、

\[ X[k] \;=\; \sum_{n=0}^{N-1} x[n]\cos\!\left(\frac{2\pi k n}{N}\right) \;-\; j \sum_{n=0}^{N-1} x[n]\sin\!\left(\frac{2\pi k n}{N}\right) \]

右辺に現れるのは、どちらも 和の形の掛け算である。信号を \( N \) 個の数を並べたベクトル \( \boldsymbol{x} \) と見て、同じ長さのベクトル

\[ \boldsymbol{c}_k = \left[\cos\frac{2\pi k \cdot 0}{N},\; \cos\frac{2\pi k \cdot 1}{N},\; \dots \right], \qquad \boldsymbol{s}_k = \left[\sin\frac{2\pi k \cdot 0}{N},\; \sin\frac{2\pi k \cdot 1}{N},\; \dots \right] \]

を用意すれば、実部は \( \boldsymbol{x}\cdot\boldsymbol{c}_k \)、虚部は \( -\,\boldsymbol{x}\cdot\boldsymbol{s}_k \) にほかならない。DFT は、信号と「周波数ビン \( k \) に対応する波」との内積を、\( k \) を変えながら並べたものである。 \(k\) 自体の単位はHzではない。

内積は「2つのものがどれだけ同じ向きを向いているか」を測る量である。統計学Ⅰ第5回の相関係数も、標準化した2変数の積の和として同じ形をしている。似ている度合いを積の和で測るという発想は共通している。

余弦と正弦の両方と内積を取るのはなぜか。信号の位相がずれていると、余弦だけとの内積はたまたま 0 になりうる。そのとき正弦との内積が大きくなる。2つを組にして

\[ |X[k]| = \sqrt{(\boldsymbol{x}\cdot\boldsymbol{c}_k)^2 + (\boldsymbol{x}\cdot\boldsymbol{s}_k)^2} \]

とすれば、その複素係数の大きさを測れる。ビンに一致する純音（DC・ナイキストの境界を除く）なら、位相を変えても大きさは変わらない。一般の有限窓では漏れや複数成分の干渉もあるため、位相の影響が常に消えるとは限らない。スペクトログラムでは通常、複素係数の位相を表示せず、\( |X[k]| \) やその2乗を描く。

\(k\) と周波数の対応

\( \boldsymbol{c}_k \) は \( N \) 点の区間にちょうど \( k \) 周期入る波である。区間の長さは \( N/f_s \) 秒なので、その周波数は \( k \div (N/f_s) = k f_s / N \)。上の \( f_k \) の式はここから出ている。

ここで示したのは「DFT が何を計算しているか」までである。どんな信号もこの形に分解し切れる（完全性）ことや、連続信号のフーリエ変換が収束する条件は、この授業で扱う数理の範囲外である。

使った数学：オイラーの式（数III・複素数）、ベクトルの内積（線形代数I）、三平方の定理。

FFT（高速フーリエ変換）

DFTの計算量は \(O(N^2)\)。FFTはアルゴリズムの工夫により \(O(N \log N)\) に短縮。

• \(N = 1024\) なら DFT は約100万回の演算、FFTは約1万回
• Numpy の np.fft.fft() は自動的にFFTを使用

PythonでFFTを計算してみよう

import numpy as np
import matplotlib.pyplot as plt

# サンプリング設定
fs = 44100       # サンプリング周波数（Hz）
duration = 0.1   # 信号の長さ（秒）
t = np.linspace(0, duration, int(fs * duration), endpoint=False)

# 440Hz と 880Hz のサイン波を合成
x = np.sin(2 * np.pi * 440 * t) + 0.5 * np.sin(2 * np.pi * 880 * t)

# FFTで周波数解析
X = np.fft.fft(x)
N = len(x)
freqs = np.fft.fftfreq(N, d=1/fs)

# 正の周波数だけ取り出す
positive = freqs > 0
plt.figure(figsize=(10, 4))
plt.plot(freqs[positive], np.abs(X[positive]))
plt.xlabel("周波数 (Hz)")
plt.ylabel("|X[k]|（未正規化）")
plt.title("FFTによる周波数スペクトル（440Hz + 880Hz）")
plt.xlim(0, 2000)
plt.grid(True)
plt.tight_layout()
plt.show()

---

3. 短時間フーリエ変換（STFT）—なぜ「短時間」なのか

問題点：全区間の振幅スペクトルだけでは、変化の時刻を読み取りにくい

全区間をまとめた振幅スペクトルは、各成分がいつ現れたかを直接示しません。ただし、複素数のDFT係数をすべて保持すれば逆変換で元の波形を復元できます。時間情報を完全に失うことと、図から局所的な変化を読み取りにくいことは別です。

例：「あいうえお」という発話を分析する場合

• FFT：音声全体の周波数成分 → 「あ」「い」「う」「え」「お」が混在
• STFT：各時刻における周波数成分 → 「今この瞬間はどの母音か」がわかる

STFTのアイデア

\[ \text{STFT}\{x\}(m, k) = \sum_{n=0}^{N-1} x[n + mH] \cdot w[n] \cdot e^{-j\frac{2\pi}{N}kn} \]

• \(w[n]\)：窓関数（後述）
• \(m\)：フレーム番号（時間位置）
• \(H\)：ホップサイズ（フレームをどれだけずらすか）
• \(N\)：この式では窓長とFFT長を等しくした点数

この式はフレームの左端を mH とする表現。実習のlibrosaは既定の center=True を使い、フレーム中心が mH に来るよう信号端をパディングする。win_length と n_fft は別々に指定できる。⁴

処理の流れ

全音声信号
  ↓ 短い区間（フレーム）に切り出す（例：25ms）
フレーム1 → FFT → 周波数スペクトル1
フレーム2 → FFT → 周波数スペクトル2
フレーム3 → FFT → 周波数スペクトル3
  ...
  ↓ 横に並べる
スペクトログラム（時間 × 周波数の2次元マップ）

フレーム長とホップサイズの設定例

パラメータ	典型的な値	意味
フレーム長	約25ms（44.1kHz なら 1,024サンプルで ≒23.2ms）	各分析窓の長さ
ホップサイズ	約10ms（44.1kHz なら 512サンプルで ≒11.6ms）	フレームのずらし幅
オーバーラップ	50〜75%	フレーム間の重複率

導出②：FFTのビン間隔と、窓長による分解能を区別する（核心）

導出①より、\( k \) 番目の成分が表す周波数は \( f_k = k f_s / N \) である。隣り合う \( k \) と \( k+1 \) の差を取ると、FFTを評価する周波数の刻み（ビン間隔）が出る。

\[ \Delta f \;=\; f_{k+1} - f_k \;=\; \frac{f_s}{N} \]

一方、窓長とFFT長が等しく、ともに \(N\) 点の場合、窓の時間長は

\[ T \;=\; \frac{N}{f_s} \]

である。この2つを掛けると、\( N \) も \( f_s \) も消える。

\[ \Delta f \times T \;=\; \frac{f_s}{N} \times \frac{N}{f_s} \;=\; 1 \qquad\Longleftrightarrow\qquad \Delta f = \frac{1}{T} \]

この積が1になるのは、窓長とFFT長を等しくしたときのビン間隔についてである。実際に近い2つの成分を分けて見られるかは、窓長・窓形状・成分の振幅差などで変わる。長い窓は一般に周波数の分離に有利だが、その窓内の時間変化を平均してしまう。短い窓は変化の時刻を追いやすいが、周波数のピークが広がる。

ゼロパディングとの区別：窓を \(L\) 点のまま、末尾などに0を補って \(N>L\) 点FFTを使うと、\(\Delta f=f_s/N\) は細かくなるが、観測時間は \(T=L/f_s\) のままである。スペクトルの描画点が増えるだけで、元の短い窓で重なった2成分を新たに分離できるわけではない。ホップ長も表示時刻の間隔であり、窓による時間分解能とは区別する。⁴

数値で確かめる（\( f_s = 16{,}000 \) Hz の場合）：

窓長 \( T \)	点数 \( N \)	周波数の刻み \( \Delta f \)
5 ms	80	200.00 Hz
25 ms	400	40.00 Hz
50 ms	800	20.00 Hz

\( f_s = 44{,}100 \) Hz で \( N = 1024 \) なら \( \Delta f = 44{,}100 / 1024 = 43.07 \) Hz である。

短い窓と長い窓は、別のものを見せる

例えばF0が100 Hzの有声音では、倍音はおよそ100 Hz間隔に並ぶ。50 msの窓では調波の横縞を分けて見やすい。5 msの窓では各調波を分けにくいが、十分短いホップで表示すると、声門の周期的な励振に対応する縦縞を見やすい。同じ音声でも窓長・窓形状・ホップで模様が変わる。前者を狭帯域、後者を広帯域のスペクトログラムと呼ぶ。

一般の時間・周波数の不確定性関係は、時間幅・周波数幅を分散などで定義した積の下限についての関係である。上のビン間隔の等式そのものではなく、その定式化はこの授業の範囲外とする。

使った数学：導出①の \( f_k \) の式、単位の約分。

---

4. 窓関数の役割

なぜ窓関数が必要か？

フレームをそのまま切り出すと（矩形窓）、端で信号が急に途切れるため、本来存在しない周波数成分（スペクトル漏れ）が発生します。

スペクトル漏れ（Spectral Leakage）の問題

矩形窓で切り出した場合：
  |────────────|    信号が端で急に0になる
  ↓FFT
  → 本来の周波数の周りに不要な成分が広がる（漏れ）

窓関数：端の重みを小さくし、漏れを抑える

ハニング窓（Hann窓）：

\[ w[n] = 0.5 \left(1 - \cos\left(\frac{2\pi n}{N-1}\right)\right) \]

これは対称形のHann窓。librosaの window='hann' は通常、スペクトル解析用の周期形（分母が \(N\)）を使う。Hamming窓の端は0にはならない。⁵

その他の主な窓関数（同じ窓長での主な傾向）：

窓関数	特徴	用途
矩形窓	主ローブが狭い・サイドローブが高い	ビンに一致する純音の分析など
ハニング窓	バランスが良い	一般的な音声処理
ハミング窓	ハニングに近い	音声処理（古典的）
ブラックマン窓	漏れ小・分解能低	精密スペクトル測定

導出③：なぜ端を 0 にすると漏れが減るのか

DFT は \( N \) 点だけを見る。このとき裏側で何が起きているかというと、その \( N \) 点が延々と繰り返されている信号を分析したことになる。導出①の \( \boldsymbol{c}_k \) が \( N \) 点でちょうど整数周期になる波しか使っていないので、そうならざるをえない。

すると、切り出した区間を周期的に延長する見方ができる。周期境界の不連続は、成分が広がって見えることを理解する手がかりになる。より一般には、時間領域で窓を掛けることで周波数領域では窓のスペクトルとの畳み込みが生じる。値の段差だけが漏れの原因ではない。

純音の周波数がちょうど \( f_k = k f_s/N \) に一致していれば、\(N\) 点の区間に整数周期が収まり、矩形窓でもそのDFTビン以外への漏れは出ない。ビンの間にある純音では複数のビンへ値が広がる。これは純音についての説明であり、一般の音声に単純な整数倍の条件を当てはめるものではない。

対称形のHann窓では、上の式で両端の値を計算すると、

\[ w[0] = 0.5\left(1 - \cos 0\right) = 0, \qquad w[N-1] = 0.5\left(1 - \cos 2\pi\right) = 0 \]

どちらもちょうど0になる。端の重みを抑えることでサイドローブが低くなる。連続的に式を見ると中央 \(n=(N-1)/2\) で \(w=1\)。ただし偶数長の対称窓では中央が標本点の間にあり、実際の配列に1は含まれない。

ただし漏れは無料では消えない。窓を掛けることは信号を変形することなので、1本の純音のピークが太くなる。表の「漏れ小・分解能低」（ブラックマン窓）と「漏れ大・分解能高」（矩形窓）は、この代償の両端を言っている。窓の選択もまた、導出②と同じくトレードオフである。

各窓のピークの太さや漏れの減り方を定量的に設計する議論は、この授業で扱う数理の範囲外である。

使った数学：余弦関数の値（数II）、導出①の \( f_k \) の式。

Pythonで窓関数を可視化

import numpy as np
import matplotlib.pyplot as plt

N = 64
n = np.arange(N)

# 各窓関数
rect   = np.ones(N)
hann   = 0.5 * (1 - np.cos(2 * np.pi * n / (N - 1)))
hamming = 0.54 - 0.46 * np.cos(2 * np.pi * n / (N - 1))

plt.figure(figsize=(10, 4))
plt.plot(rect,    label="矩形窓", linestyle='--')
plt.plot(hann,    label="ハニング窓")
plt.plot(hamming, label="ハミング窓")
plt.xlabel("サンプル番号")
plt.ylabel("重み")
plt.title("窓関数の比較")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

---

5. スペクトログラムの読み方

スペクトログラムとは

STFTの結果を2次元カラーマップとして可視化したもの。

縦軸（Y軸）：周波数（Hz）
  ↑高い音
  |
  |
  ↓低い音
  ─────────────────→ 横軸（X軸）：時間（秒）

色（濃さ）：そのフレーム・その周波数のエネルギー（強さ）
  明るい色（黄・白）→ 強い
  暗い色（紺・黒） → 弱い（またはない）

縦軸と色のスケールは別に設定する

• 周波数軸：線形Hz軸、対数周波数軸、メル軸など。低域を広く表示したいときは対数周波数軸などを使う。
• 色：振幅・パワーの線形値、または基準値との比をdBにした値。dB表示は強弱の広い範囲を見やすくする。
• この実習の amplitude_to_db(np.abs(D), ref=np.max) は、表示する振幅の最大値を0 dBとする。既定では最大値から80 dB下で表示値を切り詰める。マイクの校正がないため、音圧レベル（dB SPL）ではない。⁶

デシベル（dB）とは

\[ P_{dB} = 10 \log_{10} \frac{P}{P_{\text{ref}}} \]

• 0dB：基準レベル
• 60dB：基準の100万倍のエネルギー
• -20dB：基準の1/100のエネルギー

---

6. 有声音・無声音・無音のスペクトログラムの違い

3種類の音の特徴

種類	例	スペクトログラムの特徴
有声音	母音の多く、/m/・/n/など	狭帯域表示では調波構造。広帯域表示では周期的な縦縞も見える
無声音	「さ」の子音/s/、「し」の子音[ɕ]、/t/・/k/の閉鎖・破裂など	摩擦音では雑音状、破裂では短いバースト、閉鎖では低エネルギー
無音・背景音のみ	発話の間	声はないが環境雑音が残る場合がある

「さ」全体は子音と母音を含むため、一律に無声音とは呼ばない。背景雑音と無声音の区別も、色の明暗だけでは決められない。

有声音の特徴：調波構造（Harmonic Structure）

• 基本周波数（F0）：準周期的な声帯振動の繰り返し周波数。話者・年齢・発声条件で変わる
• 倍音（高調波）：F0の整数倍（2F0, 3F0, 4F0, ...）に強いエネルギー
• 十分長い窓では、線形Hz軸上におおよそ等間隔の調波として見える。F0が変化すれば縞も曲がる

周波数
  ↑
4F0 ─────────────  (倍音)
3F0 ─────────────  (倍音)
2F0 ─────────────  (倍音)
 F0 ─────────────  (基本周波数)
  ─────────────────────→ 時間

無声音の特徴

• 声帯振動なし → 調波構造がない
• 摩擦音には持続する雑音状成分、破裂音には閉鎖と短いバーストがある
• /s/の摩擦は高域に強い成分を持ちやすいが、分布は調音・話者・録音帯域に依存する

---

7. 日本語の母音「あいうえお」のスペクトログラム比較

フォルマントとは

声道（口や喉の形）により強調される周波数帯域をフォルマントといいます。

• F1（第1フォルマント）：おおむね舌が低く、開口度が大きいほど高くなる
• F2（第2フォルマント）：おおむね前舌母音で高く、後舌・円唇母音で低くなる

これらは傾向であり、声道形状との一対一の対応ではない。フォルマントはスペクトル包絡の共鳴の山で、F0の整数倍に並ぶ個々の調波とは区別する。

実習の合成母音に使うフォルマント設定例

以下は違いを可視化するための模式的な値であり、日本語話者の平均値・判定基準ではない。実際の値は話者・文脈・発声で変わる。

母音	F1 (Hz)	F2 (Hz)	口の形
あ /a/	800	1200	大きく開く
い /i/	300	2500	横に広げる
う /u/	400	1400	唇を強く突き出さない。標準語では非円唇寄りの[ɯ]など7
え /e/	500	1900	中間
お /o/	500	800	丸く突き出す

スペクトログラムで見分けるポイント

• 「あ」：F1が高い（800Hz）→ スペクトログラムの低〜中域に強い縞
• 「い」：F2が高い設定（2500Hz）→ その付近の調波群が強調される
• 「お」：F2が低い（800Hz）→ 全体的に低域に集中

---

8. Python実習—スペクトログラムの描画

環境設定

Google Colab でセルを実行してください。

# ===== セル1：必要なライブラリのインストール =====
!pip install librosa soundfile matplotlib numpy

# ライブラリのインポート
import librosa
import librosa.display
import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf
from IPython.display import Audio, display

# ===== セル2：サンプル音声の読み込みとスペクトログラム描画 =====

# librosa の組み込みサンプル音声を使用
audio_path = librosa.example('libri1')    # LibriSpeechの朗読音声。初回はダウンロード
y, sr = librosa.load(audio_path, sr=None) # y：波形データ, sr：サンプリングレート

print(f"サンプリングレート: {sr} Hz")
print(f"音声の長さ: {len(y)/sr:.2f} 秒")
print(f"サンプル数: {len(y)}")

# 音声を再生（Google Colab）
display(Audio(y, rate=sr))

# 波形を表示
plt.figure(figsize=(12, 3))
librosa.display.waveshow(y, sr=sr)
plt.title("波形（時間領域）")
plt.xlabel("時間 (秒)")
plt.ylabel("振幅")
plt.tight_layout()
plt.show()

# ===== セル3：STFTを計算してスペクトログラムを描画 =====

# STFTのパラメータ
n_fft   = 2048   # FFT長。この例ではwin_lengthも同じ値
hop_length = 512 # ホップサイズ（フレームのずらし幅）
window = 'hann'  # 窓関数

# STFTを計算
D = librosa.stft(y, n_fft=n_fft, hop_length=hop_length, window=window)
# D の形状：(n_fft//2 + 1, フレーム数) の複素数配列（モノラル入力）
# center=True（既定）：フレームの中心が m*hop_length に対応

# 振幅スペクトログラムをdBスケールに変換
S_db = librosa.amplitude_to_db(np.abs(D), ref=np.max)

print(f"スペクトログラムの形状: {S_db.shape}")
print(f"  周波数ビン数: {S_db.shape[0]}")
print(f"  フレーム数: {S_db.shape[1]}")

# 描画
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 上段：線形スケール
S_linear = np.abs(D)
img1 = librosa.display.specshow(S_linear, sr=sr, hop_length=hop_length,
                                  x_axis='time', y_axis='hz', ax=axes[0])
axes[0].set_title("スペクトログラム（線形スケール）")
fig.colorbar(img1, ax=axes[0], format='%+.0f')

# 下段：対数スケール（dB）
img2 = librosa.display.specshow(S_db, sr=sr, hop_length=hop_length,
                                  x_axis='time', y_axis='log', ax=axes[1],
                                  cmap='magma')
axes[1].set_title("スペクトログラム（対数スケール、dB）")
fig.colorbar(img2, ax=axes[1], format='%+.0f dB')

plt.tight_layout()
plt.show()

# ===== セル4：母音「あいうえお」のスペクトログラム比較 =====
# ※ここでは合成音声でシミュレート（実際の授業では録音データを使用）

import numpy as np
import matplotlib.pyplot as plt

def make_vowel(f0, formants, duration=0.5, sr=22050):
    """調波の振幅を共鳴の包絡で重み付けした、教育用の定常母音モデル。"""
    if not 0 < f0 < sr / 2 or duration <= 0:
        raise ValueError("f0とdurationを確認してください")
    t = np.arange(int(sr * duration)) / sr
    # ナイキスト周波数未満の調波だけを生成
    k = np.arange(1, int(np.ceil((sr / 2) / f0)))
    freqs = k * f0
    envelope = np.full_like(freqs, 0.03, dtype=float)
    for formant, width in zip(formants, [80, 110, 150]):
        envelope += np.exp(-0.5 * ((freqs - formant) / width) ** 2)
    amplitudes = envelope / k  # 音源の1/k減衰 × 共鳴の包絡
    signal = np.sum(amplitudes[:, None] *
                    np.sin(2 * np.pi * freqs[:, None] * t), axis=0)
    # 端のクリックを和らげ、各音を個別にピーク正規化
    signal *= np.minimum(1, t / 0.01) * np.minimum(1, (duration - t) / 0.01)
    return signal / max(np.max(np.abs(signal)), 1e-12)

# 母音ごとの模式的なフォルマント設定（測定値・平均値ではない）
vowels = {
    'あ (a)': [800, 1200, 2600],
    'い (i)': [300, 2500, 3200],
    'う (u)': [400, 1400, 2400],
    'え (e)': [500, 1900, 2700],
    'お (o)': [500,  800, 2400],
}

f0 = 150   # 比較用に全母音で同じF0を使う。性別を指定する値ではない
sr = 22050
n_fft = 1024
hop_length = 256

fig, axes = plt.subplots(2, 5, figsize=(18, 8))

for col, (vowel_name, formants) in enumerate(vowels.items()):
    signal = make_vowel(f0, formants, duration=0.5, sr=sr)

    # 波形
    axes[0, col].plot(np.arange(len(signal)) / sr, signal, linewidth=0.5)
    axes[0, col].set_title(f"波形\n{vowel_name}", fontsize=10)
    axes[0, col].set_xlabel("時間 (秒)")
    axes[0, col].set_xlim(0, 0.5)
    if col == 0:
        axes[0, col].set_ylabel("振幅")

    # スペクトログラム
    D = librosa.stft(signal, n_fft=n_fft, hop_length=hop_length)
    S_db = librosa.amplitude_to_db(np.abs(D), ref=np.max)
    img = librosa.display.specshow(S_db, sr=sr, hop_length=hop_length,
                                    x_axis='time', y_axis='hz',
                                    ax=axes[1, col], cmap='magma')
    axes[1, col].set_title(f"スペクトログラム\n{vowel_name}", fontsize=10)
    axes[1, col].set_ylim(0, 5000)

fig.suptitle("日本語母音「あいうえお」のスペクトログラム比較（合成音声）", fontsize=13)
plt.tight_layout()
plt.show()

# ===== セル5：Colabのブラウザで3秒録音する =====
from google.colab import output

# マイク許可後に録音。Promiseの完了を待ってPythonへデータを返す
record_js = r'''
(async () => {
  const stream = await navigator.mediaDevices.getUserMedia({audio: true});
  try {
    const recorder = new MediaRecorder(stream);
    const chunks = [];
    return await new Promise((resolve, reject) => {
      recorder.ondataavailable = e => {
        if (e.data.size > 0) chunks.push(e.data);
      };
      recorder.onerror = e => reject(e.error || new Error("録音エラー"));
      recorder.onstop = () => {
        const blob = new Blob(chunks, {type: recorder.mimeType});
        const reader = new FileReader();
        reader.onloadend = () => resolve(reader.result);
        reader.onerror = () => reject(reader.error);
        reader.readAsDataURL(blob);
      };
      recorder.start();
      setTimeout(() => {
        if (recorder.state === "recording") recorder.stop();
      }, 3000);
    });
  } finally {
    stream.getTracks().forEach(track => track.stop());
  }
})()
'''
print("マイクを許可したら、3秒間「あいうえお」などと話してください。")
audio_data = output.eval_js(record_js)
print("録音完了。次のセルで変換・分析します。")

# ===== セル6：録音した音声をWAVへ変換し、可視化 =====
import base64
import shutil
import subprocess
import tempfile
from pathlib import Path
import librosa
import librosa.display
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Audio, display

if shutil.which("ffmpeg") is None:
    raise RuntimeError("別セルで !apt-get install -y ffmpeg を実行してください")

# ブラウザの録音形式は環境によって異なる。内容をFFmpegに判別させる
_, encoded = audio_data.split(",", 1)
with tempfile.TemporaryDirectory() as tmp:
    source = Path(tmp) / "recording_input.bin"
    target = Path(tmp) / "recording.wav"
    source.write_bytes(base64.b64decode(encoded))
    subprocess.run(
        ["ffmpeg", "-y", "-i", str(source), "-ac", "1", "-ar", "22050",
         "-c:a", "pcm_s16le", str(target)],
        check=True, stdout=subprocess.DEVNULL, stderr=subprocess.PIPE)
    y, sr = librosa.load(target, sr=None, mono=True)

print(f"録音長: {len(y)/sr:.2f} 秒, サンプリングレート: {sr} Hz")
display(Audio(y, rate=sr))
fig, axes = plt.subplots(3, 1, figsize=(12, 10))
librosa.display.waveshow(y, sr=sr, ax=axes[0])
axes[0].set_title("録音した音声の波形")

D = librosa.stft(y, n_fft=2048, hop_length=512)
S_db = librosa.amplitude_to_db(np.abs(D), ref=np.max)
img = librosa.display.specshow(S_db, sr=sr, hop_length=512,
                              x_axis="time", y_axis="log",
                              ax=axes[1], cmap="magma")
axes[1].set_title("スペクトログラム（対数周波数軸・相対dB）")
fig.colorbar(img, ax=axes[1], format="%+.0f dB")

S_mel = librosa.feature.melspectrogram(y=y, sr=sr, n_mels=128,
                                     n_fft=2048, hop_length=512)
S_mel_db = librosa.power_to_db(S_mel, ref=np.max)
img2 = librosa.display.specshow(S_mel_db, sr=sr, hop_length=512,
                               x_axis="time", y_axis="mel",
                               ax=axes[2], cmap="magma")
axes[2].set_title("メルスペクトログラム（次回の先取り）")
fig.colorbar(img2, ax=axes[2], format="%+.0f dB")
plt.tight_layout()
plt.show()

マイクを使えない環境では、手元のWAVをColabへアップロードし、y, sr = librosa.load("ファイル名.wav", sr=None) として分析部分を実行する。録音セルはColabのブラウザ機能に依存する。上の合成母音は共鳴による調波の強調を学ぶための簡略モデルで、自然音声の完全な合成ではない。各母音・各図を個別に最大値で正規化しているので、色の明暗から母音間の絶対的な音量を比べない。

---

9. 演習問題と考察

演習課題（授業中）

課題1（個人）：スペクトログラムを読み取ろう

librosa のサンプル音声を使い、以下を答えなさい：

1. 有声音と無声音の区間を時間で特定せよ
2. 基本周波数（F0）をスペクトログラムから目視で推定せよ
3. 窓関数を「矩形窓」と「ハニング窓」に変えて、スペクトログラムの違いを観察せよ

課題2（任意・できる範囲で）：自分の声を録音して比較

1. 「あ」と「い」を約1秒ずつ発声して録音
2. スペクトログラムを描画し、F1・F2の位置の違いをスクリーンショットに矢印で示す
3. やれる場合はZoomチャットに貼り付けて共有（メールで送ってもよい）

---

10. まとめと次回予告

本日のまとめ

概念	ポイント
フーリエ変換	音 = サイン波の重ね合わせ。FFTで効率的に計算
STFT	「短時間」ごとにFFT → 時間変化を追える
窓関数	サイドローブの漏れを抑えるが、主ローブは広がる
スペクトログラム	横=時間、縦=周波数、色=強度の2次元マップ
有声音	窓と軸の設定に応じて調波・周期性を観察できる
無声音	摩擦・閉鎖・破裂など、生成の仕組みによって見え方が異なる

次回（第5回・動画回）予告：MFCC

スペクトログラムの情報を人間の聴覚特性に合わせて圧縮する。 従来の音声認識で広く使われてきた特徴量 MFCC を学びます。対数メルスペクトログラムとの違いも確認します。

予習課題（やれる範囲で）

• 「メル尺度」「ケプストラム」という言葉を調べておく
• 第5回の動画教材を視聴してからテストに取り組む

---

参考文献・資料

• Müller, M. (2015). Fundamentals of Music Processing. Springer.
• librosa 公式ドキュメント: https://librosa.org/doc/
• McFee, B. et al. (2015). librosa: Audio and music signal analysis in Python. Proceedings of the 14th Python in Science Conference.

内容確認・訂正：2026年10月20日。

---

1. NIDCD, How Do We Hear?; Slaney (1998), Auditory Toolbox（STFTと蝸牛モデルを区別）。 ↩

2. Fitch, W. T. (1997). Vocal tract length and formant frequency dispersion correlate with body size in rhesus macaques. JASA, 102, 1213–1222. 原論文。 ↩

3. NumPy, Discrete Fourier Transform（定義・周波数の順序・正規化）。 ↩

4. librosa 0.11.0, stft（窓長・FFT長・ホップ・center）。 ↩↩

5. SciPy, hann（周期形と対称形）。 ↩

6. librosa 0.11.0, amplitude_to_db（基準値と表示下限）。 ↩

7. 東京外国語大学言語モジュール、日本語の母音「う」。 ↩