Pythonの音声処理ライブラリ【LibROSA】で音声読み込み⇒スペクトログラム変換・表示⇒位相推定して音声復元


# LibROSAとは

- LibROSAはPythonの音声処理ライブラリです。
- 様々な音声処理を簡潔に記述できます。
- 今回は以下の音声処理の基本処理をまとめました。
  - 音声の読み込み
    - 周波数を指定して音声を読み込み
    - Notebook上で、音声をプレーヤーで再生
    - 音声波形のグラフを表示
  - スペクトログラムへの変換
    - STFTで音声からスペクトログラムへ変換
    - 強度をdB単位に変換
    - スペクトログラムのカラープロットを表示
  - 音声を復元
    - 逆STFTでスペクトログラムから音声を復元する場合
    - 位相情報を推定して音声を復元する場合

ソースコード：https://github.com/lilacs2039/ColabNotebooks/blob/master/audio/LibROSA%E4%BD%BF%E3%81%84%E6%96%B9.ipynb

# 初期化


```python
import os
from tqdm import tqdm
import numpy as np
import matplotlib.pyplot as plt
import IPython.display
from IPython.display import display
import pandas as pd

import librosa
import librosa.display
```


```python
# pyplotのデフォルト値を設定
plt.rcParams.update({
  'font.size' : 10
  ,'font.family' : 'Meiryo' if os.name == 'nt' else ''  # Colabでは日本語フォントがインストールされてないので注意
  ,'figure.figsize' : [10.0, 5.0]
  ,'figure.dpi' : 300
  ,'savefig.dpi' : 300
  ,'figure.titlesize' : 'large'
  ,'legend.fontsize' : 'small'
  ,'axes.labelsize' : 'medium'
  ,'xtick.labelsize' : 'small'
  ,'ytick.labelsize' : 'small'
  })
# ndarrayの表示設定
np.set_printoptions(threshold=0)  # 可能ならndarrayを省略して表示
np.set_printoptions(edgeitems=1)  # 省略時に１つの要素だけ表示
```

# 音声の確認

## 読み込み
librosaに標準で入っている音声ファイルを読込みます。[^1]


```python
audio_path = librosa.util.example_audio_file(); audio_path
```



```bash:出力
    '（ホームディレクトリなど）\\anaconda3\\lib\\site-packages\\librosa\\util\\example_data\\Kevin_MacLeod_-_Vibe_Ace.ogg'
```

librosa.load() ドキュメント　https://librosa.org/librosa/master/generated/librosa.core.load.html


```python
y, sr = librosa.load(audio_path)  # サンプリング周波数　22.05kHzで読み込み
# y, sr = librosa.load(librosa.util.example_audio_file(), sr=None)  # 元の音声ファイルのサンプリング周波数で読み込む場合
# y, sr = librosa.load(librosa.util.example_audio_file(), sr=4096)  # 約4kHzでリサンプリングして読み込む場合
print([type(y), y.shape], [type(sr), sr])

```

```bash:出力
    [<class 'numpy.ndarray'>, (1355168,)] [<class 'int'>, 22050]
```    

サンプリング周波数`sr`について
- `sr`はSampling Rate（サンプリング周波数）の略です。１秒間の音声を何個のデータで表すかを示します。（↑の例では22050データ/秒）
  - **全データ数`y.size` = `sr` × 再生時間**　が成り立ちます。
- srを指定しない場合、読み込んだファイルのサンプリング周波数に関わらず、**22.05kHzにリサンプリングして読み込みます。　←　【初見注意】**
- `sr=None`　を指定すると、元ファイルのサンプリング周波数で読み込みます。
- `sr=値`　を指定すると、指定したサンプリング周波数で読み込みます。


## 音声をプレーヤーで確認

- IPython.display.display()　ドキュメント https://ipython.readthedocs.io/en/stable/api/generated/IPython.display.html#IPython.display.display
- IPython.display.Audio() ドキュメント https://ipython.readthedocs.io/en/stable/api/generated/IPython.display.html#IPython.display.Audio

音声はBase64エンコーディングしてJupyter Notebookに埋め込まれます。（重いので本記事では置き換えました）


```python
display(IPython.display.Audio(y, rate=sr))
```

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">librosa音声変換記事 <a href="https://t.co/CedmBGtHfS">https://t.co/CedmBGtHfS</a><br><br>Vibe Ace (Kevin MacLeod) / CC BY 3.0 <a href="https://t.co/eMkYzYU3af">pic.twitter.com/eMkYzYU3af</a></p>&mdash; lilacs2039 (@lilacs2039) <a href="https://twitter.com/lilacs2039/status/1279752555782828032?ref_src=twsrc%5Etfw">July 5, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

## 音声波形表示

librosa.display.waveplot() ドキュメント　https://librosa.org/librosa/master/generated/librosa.display.waveplot.html


```python
mpl_collection = librosa.display.waveplot(y, sr=sr)
mpl_collection.axes.set(title="音声波形", ylabel="波形の振幅")
```



```bash:出力
    [Text(0, 0.5, '波形の振幅'), Text(0.5, 1.0, '音声波形')]
```



![output_17_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/ebf09e9f-e1f5-54d3-6f50-dd7c5e28d8fa.png)


# スペクトログラムへの変換

- 音声信号に「短時間フーリエ変換（Short-time Fourier transform : STFT）」を行うことでスペクトログラムが得られます。
- スペクトログラムとは、音声を周波数スペクトルの時間経過で示したものです。

## スペクトログラムの計算・表示


```python
D = librosa.stft(y)  # STFT
S, phase = librosa.magphase(D)  # 複素数を強度と位相へ変換
Sdb = librosa.amplitude_to_db(S)  # 強度をdb単位へ変換
librosa.display.specshow(Sdb, sr=sr, x_axis='time', y_axis='log')  # スペクトログラムを表示
```



```bash:出力
    <matplotlib.axes._subplots.AxesSubplot at 0x1cb80c833c8>
```



![output_21_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/b3873dc8-f797-511d-246b-94c3aa43e8f1.png)


## やっていること

### STFT

- STFTの手順
  - 音声信号のある時間範囲を窓で区切ってフーリエ変換して、その時間範囲の周波数スペクトルを計算します。
  - 窓をずらしながら繰り返すことで、周波数スペクトルの時間経過の行列`D`が得られます。
- STFTの主なパラメータ
  - `n_fft`（窓の長さ）
    - デフォルト値は`2048`
  - `hop_length`（窓の移動幅）
    - デフォルト値は`n_fft/4`

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/3335f616-0719-b044-14b4-d51f74efb3b4.png)

#### 周波数ビン・時刻ビンの計算方法

スペクトログラム計算結果のイメージ

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/6044dbbc-9342-f50c-c9c8-c7418d9197f6.png)


- 「周波数ビンの数」「時刻ビンの数」「周波数ビンの幅」「時刻ビンの幅」は以下のように計算できます。
- サンプリング定理から、音声の周波数はサンプリング周波数/2までしか表現できません。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/afe04726-aecc-9454-0e69-c0110a289847.png)


- 計算例
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/24255c8d-a0aa-cb73-8bbd-8dca0e14a14c.png)


- デフォルト値において、`librosa.stft()`は以下のような仕様になることがわかります。
  - 0Hz～約11kHzの範囲を約1025分割して表現（周波数分解能は約10.76Hz）
  - １秒間を約43分割して表現（時間分解能は約23.26ms）
- `sr`や`n_fft`を調整して目標の分解能やデータサイズが実現できます。

#### STFTの出力の行列`D`の確認

librosa.stft()　ドキュメント https://librosa.org/librosa/master/generated/librosa.core.stft.html


```python
D = librosa.stft(y)  # STFT
print(type(D), D.shape)
print(D)
```

```bash:出力
    <class 'numpy.ndarray'> (1025, 2647)
    [[ 2.5802802e-03-0.j ...  7.2303657e-05-0.j]
     ...
     [-1.7063057e-07-0.j ... -5.1327298e-09-0.j]]
```    

- STFTで得られた`D`は複素数の２次元ndarray ※ であることがわかります。

※　「２次元」ではなく「２階のテンソル」と呼ぶ方が正しいかもしれませんが、本記事では２次元で統一します。

### 複素数を強度と位相へ変換
複素数の行列`D`を、強度の行列`S`と位相の行列`phase`に変換します。
- librosa.magphase() ドキュメント https://librosa.org/librosa/master/generated/librosa.core.magphase.html


```python
S, phase = librosa.magphase(D)  # 複素数を強度と位相へ変換
```

- `librosa.magphase(D)`がやっていること
  - STFTで得られた行列Dの複素数は直交形式ですが、スペクトログラムは極形式（絶対値 r: 強度, 偏角 θ: 位相）が望ましいです。
  - そこで、複素数を「直交形式から極形式へ変換」しています。
  - 参考：複素平面 ～ wikipedia  https://ja.wikipedia.org/wiki/%E8%A4%87%E7%B4%A0%E5%B9%B3%E9%9D%A2
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/2ea153bc-ed33-7e17-ee98-299b521625e2.png)


#### 強度の行列`S`の確認


```python
print(type(S), S.shape)
print(S)
```

```bash:出力
    <class 'numpy.ndarray'> (1025, 2647)
    [[2.5802802e-03 ... 7.2303657e-05]
     ...
     [1.7063057e-07 ... 5.1327298e-09]]
```    


```python
s = pd.DataFrame(S.flatten())
s.hist(bins=20, range=(s.min().values[0],s.quantile(0.9).values[0]) ) # Sの最小値～90%点までの分布を表示
plt.title("強度の分布")
```



```bash:出力
    Text(0.5, 1.0, '強度の分布')
```



![output_42_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/cdeea6e7-30ac-4948-5c17-d9b45e7ecbd0.png)


強度の行列`S`は、２次元ndarrayで各要素は0より大きい実数であることがわかります。

#### 位相の行列`phase`の確認


```python
print(type(phase), phase.shape)
print(phase)
```

```bash:出力
    <class 'numpy.ndarray'> (1025, 2647)
    [[ 1.+0.000000e+00j ...  1.+0.000000e+00j]
     ...
     [-1.+8.742278e-08j ... -1.+8.742278e-08j]]
```    


```python
p = pd.DataFrame(phase.flatten()).sample(n=1000)  # 散布図にするにはデータ数多すぎなので、1000データをランダムサンプリング
mpl_collection = plt.scatter(np.real(p), np.imag(p))
mpl_collection.axes.set(title="複素平面上の位相の散布図", xlabel="実部", ylabel="虚部", aspect='equal')
```



```bash:出力
    [Text(0, 0.5, '虚部'), Text(0.5, 0, '実部'), Text(0.5, 1.0, '複素平面上の位相の散布図'), None]
```



![output_46_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/84a5ec8a-3078-3248-78ab-9fc07ed41e60.png)


位相の行列`phase`は、２次元ndarrayで各要素は絶対値１, 偏角θの複素数であることがわかります。

### 強度をdb単位へ変換

- スペクトログラムの強度が線形スケールで扱いづらい場合は、dB単位に変換して対数スケールで扱えるようにします。
- dB単位では、+20dBで（振幅）スペクトログラムの強度は10倍になります。

librosa.amplitude_to_db() ドキュメント https://librosa.org/librosa/master/generated/librosa.core.amplitude_to_db.html


```python
Sdb = librosa.amplitude_to_db(S)  # 強度をdb単位へ変換
```

#### `librosa.amplitude_to_db()`がやっていること

下式で、振幅スペクトログラムをdB単位に変換します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/255dcdc6-7c98-3783-57a1-52d33f64a036.png)


- 主なパラメータ
  - 最小値`amin`
    - 計算式のlogの中身が0のとき戻り値が-∞とならないように、最小値を指定
    - デフォルト値　1e-5
  - 基準値`ref`
    - スペクトログラムの２乗 = refのとき、0dBとなるように変換する。
    - デフォルト値 1

#### dB単位の強度の行列`Sdb`の確認


```python
print(type(Sdb), Sdb.shape)
print(Sdb)
```

```bash:出力
    <class 'numpy.ndarray'> (1025, 2647)
    [[-33.29261 ... -33.29261]
     ...
     [-33.29261 ... -33.29261]]
```    


```python
sdb = pd.DataFrame(Sdb.flatten())
sdb.hist(bins=20, range=(sdb.min().values[0], sdb.quantile(0.9).values[0]) ) # Sの最小値～90%点までの分布を表示
plt.title("dB単位の強度の分布")
```




```bash:出力
    Text(0.5, 1.0, 'dB単位の強度の分布')
```



![output_58_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/240b965e-cc16-e709-b358-dd9c7ac7446a.png)


### スペクトログラムを表示

- 周波数と強さは線型目盛でも対数目盛でもよく、用途によって使い分けます。

librosa.display.specshow() ドキュメント https://librosa.org/librosa/master/generated/librosa.display.specshow.html

#### 周波数：線形スケール

`librosa.display.specshow()`の引数に、`y_axis='hz'`を指定します


```python
fig,axes = plt.subplots(nrows=1,ncols=2)
axes[0].set_title("強度：線形スケール")
librosa.display.specshow(S, sr=sr, x_axis='time', y_axis='hz', ax=axes[0])
axes[1].set_title("強度：対数スケール（dB単位）")
librosa.display.specshow(Sdb, sr=sr, x_axis='time', y_axis='hz', ax=axes[1])
```




```bash:出力
    <matplotlib.axes._subplots.AxesSubplot at 0x1cb81bed2c8>
```



![output_64_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/a09705bc-29ec-bae9-dccf-bd74f94718c9.png)


#### 周波数：対数スケール

`librosa.display.specshow()`の引数に、`y_axis='log'`を指定します


```python
fig,axes = plt.subplots(nrows=1,ncols=2)
axes[0].set_title("強度：線形スケール")
librosa.display.specshow(S, sr=sr, x_axis='time', y_axis='log', ax=axes[0])
axes[1].set_title("強度：対数スケール（dB単位）")
librosa.display.specshow(Sdb, sr=sr, x_axis='time', y_axis='log', ax=axes[1])
```




```bash:出力
    <matplotlib.axes._subplots.AxesSubplot at 0x1cb81ca9408>
```



![output_67_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/3464c392-a1de-2f2e-a092-d375b117a87d.png)


#### （おまけ）位相スペクトログラム


```python
# 位相`phase`（直交形式）から偏角θを計算してカラープロット
ax = librosa.display.specshow(np.arctan2(np.imag(phase), np.real(phase)), sr=sr, x_axis='time', y_axis='hz')
ax.axes.set_ylabel("Phase")
plt.title("位相スペクトログラム")
```




```bash:出力
    Text(0.5, 1.0, '位相スペクトログラム')
```



![output_69_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/524ddb9d-0b2c-bb0e-28fa-c817d59ff988.png)


ランダムノイズのようにしか見えません。パターンを見出すのは難しそうです。

## 周波数・時間分解能の調整例
- やること
  - 最初の60秒に対して、0~1kHzの範囲を時間分解能10ms程度で確認したい
- 手順
  - 約2kHzでリサンプリング　→　n_fftを調整してSTFT、強度スペクトログラムを表示（周波数は線形スケール、強度はdB単位）


- スペクトログラムのビン形状・分解能の予測値

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/224d3834-a8c6-f330-90db-b04dc3b3d79e.png)


```python
playtime2, sr2, n_fft2, hop_length2 = 60, 2048, 64, 64//4  # パラメータ
y2 = librosa.resample(y[0:sr * playtime2],sr,sr2)
D2 = librosa.stft(y2,n_fft=n_fft2, hop_length=hop_length2)  # hop_lengthはデフォルト値と同じ計算にしたので、指定しなくても同じ。
S2, _ = librosa.magphase(D2)  # 複素数を強度と位相へ変換
Sdb2 = librosa.amplitude_to_db(S2)  # 強度をdb単位へ変換
tmp = librosa.display.specshow(Sdb2, sr=sr2, hop_length=hop_length2, x_axis='time', y_axis='hz')  # スペクトログラムを表示
plt.title("スペクトログラム（周波数：線形スケール, 強度：dB単位）")
plt.colorbar(format='%+2.0f dB')
```




```bash:出力
    <matplotlib.colorbar.Colorbar at 0x1cbb67bec08>
```



![output_73_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/9e4b18c4-d4c1-53c6-c020-be12a98b85a5.png)



```python
print("予測ビン形状 : \t", (33, 128 * playtime2))
print("実際のビン形状 :\t", Sdb2.shape)
```

```bash:出力
    予測ビン形状 : 	 (33, 7680)
    実際のビン形状 :	 (33, 7681)
```    

スペクトログラムのビンの形状は、計算式で予測した値と大体あっていることが確認できます。

# 音声の復元

## 位相情報が使える場合：librosa.istft()

- 手順
  - 下式で、強度`S`,位相`phase`から直交形式の複素数の行列`D`へ変換
  - ISTFTでスペクトログラムを音声に復元します。
- 逆短時間フーリエ変換（Inverse short-time Fourier transform : ISTFT) 
- librosa.istft() ドキュメント https://librosa.org/librosa/master/generated/librosa.core.istft.html?highlight=istft#librosa.core.istft

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/99090399-935b-d20c-65a8-d9722fb6dd78.png)


```python
D = S * np.exp(1j*phase)  # 直交形式への変換はlibrosaの関数ないみたいなので、自分で計算する。
y_inv = librosa.istft(D)
display(IPython.display.Audio(y_inv, rate=sr))
mpl_collection = librosa.display.waveplot(y_inv, sr=sr)
mpl_collection.axes.set(title="位相情報を使って復元した音声波形", ylabel="波形の振幅")
```

```bash:出力
    [Text(0, 0.5, '波形の振幅'), Text(0.5, 1.0, '位相情報を使って復元した音声波形')]
```

![output_80_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/33a6449c-d7c9-359c-847b-ae2e934688f1.png)

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">librosa音声変換記事 <a href="https://t.co/CedmBGtHfS">https://t.co/CedmBGtHfS</a><br><br>Vibe Ace (Kevin MacLeod) / CC BY 3.0 <a href="https://t.co/qaTaR4u7dt">pic.twitter.com/qaTaR4u7dt</a></p>&mdash; lilacs2039 (@lilacs2039) <a href="https://twitter.com/lilacs2039/status/1279753135402053632?ref_src=twsrc%5Etfw">July 5, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

```python
ylim = ax.axes.set_ylim(); ylim  # 次の比較用にy軸の範囲を保存
```


```bash:出力
    (-0.6926389932632446, 0.6926389932632446)
```


## 位相情報が使えない場合：griffinlim()で位相推定

### 位相情報なしで復元した場合

位相をすべてゼロとして直交形式の複素数の行列`D`へ変換した場合です。


```python
D = S * np.exp(np.zeros(S.shape))
y_inv = librosa.istft(D)
display(IPython.display.Audio(y_inv, rate=sr))
mpl_collection = librosa.display.waveplot(y_inv, sr=sr)
mpl_collection.axes.set(title="位相０で復元した音声波形", ylabel="波形の振幅")
mpl_collection.axes.set_ylim(ylim);
```

![output_85_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/4718be84-09d5-0d8b-659d-4b6eab6e2026.png)

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">librosa音声変換記事 <a href="https://t.co/CedmBGtHfS">https://t.co/CedmBGtHfS</a><br><br>Vibe Ace (Kevin MacLeod) / CC BY 3.0 <a href="https://t.co/UzAxBfAFgt">pic.twitter.com/UzAxBfAFgt</a></p>&mdash; lilacs2039 (@lilacs2039) <a href="https://twitter.com/lilacs2039/status/1279753497387253762?ref_src=twsrc%5Etfw">July 5, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>



- 音声の振幅が小さくなっています。聞いてみても何かヘンです。

### griffinlim()で位相推定した場合

- librosa 0.7.x以降のバージョンでは`librosa.griffinlim()`で、強度スペクトログラムのみから位相を推定して音声を復元できます。
  - librosa.griffinlim() ドキュメント　https://librosa.org/librosa/master/generated/librosa.core.griffinlim.html#librosa-core-griffinlim
- librosa 0.7.x　未満を使っているときは、issueから対象のソースコードをコピペして実行できます。（今回はこの方法です）
  - Synthetic phase for inverse transforms #434
    - https://github.com/librosa/librosa/issues/434#issuecomment-291266229


```python
def griffinlim(spectrogram, n_iter = 100, window = 'hann', n_fft = 2048, hop_length = -1, verbose = False):
    if hop_length == -1:
        hop_length = n_fft // 4

    angles = np.exp(2j * np.pi * np.random.rand(*spectrogram.shape))

    t = tqdm(range(n_iter), ncols=100, mininterval=2.0, disable=not verbose)
    for i in t:
        full = np.abs(spectrogram).astype(np.complex) * angles
        inverse = librosa.istft(full, hop_length = hop_length, window = window)
        rebuilt = librosa.stft(inverse, n_fft = n_fft, hop_length = hop_length, window = window)
        angles = np.exp(1j * np.angle(rebuilt))

        if verbose:
            diff = np.abs(spectrogram) - np.abs(rebuilt)
            t.set_postfix(loss=np.linalg.norm(diff, 'fro'))

    full = np.abs(spectrogram).astype(np.complex) * angles
    inverse = librosa.istft(full, hop_length = hop_length, window = window)

    return inverse
```


```python
print("位相推定　開始")
y_inv = griffinlim(S, verbose=True)
print("位相推定　終了")
display(IPython.display.Audio(y_inv, rate=sr))
mpl_collection = librosa.display.waveplot(y_inv, sr=sr)
mpl_collection.axes.set(title="griffinlimで位相推定して復元した音声波形", ylabel="波形の振幅")
mpl_collection.axes.set_ylim(ylim);
```

```bash:出力
    位相推定　開始
    100%|███████████████████████████████████████████████████| 100/100 [01:04<00:00,  1.55it/s, loss=376]
    位相推定　終了
```    

![output_90_4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/186338/6f7fcf2f-f5da-0eb2-8cf0-30eb034dffa4.png)

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">librosa音声変換記事 <a href="https://t.co/CedmBGtHfS">https://t.co/CedmBGtHfS</a><br><br>Vibe Ace (Kevin MacLeod) / CC BY 3.0 <a href="https://t.co/k9M3LXCnGM">pic.twitter.com/k9M3LXCnGM</a></p>&mdash; lilacs2039 (@lilacs2039) <a href="https://twitter.com/lilacs2039/status/1279753735955070977?ref_src=twsrc%5Etfw">July 5, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

- 復元結果
  - 聞いた感じは遜色なく復元できています。
  - 音声波形はよく見ると、元の音声とは少し異なっています。
- 使いどころ
  - 音声認識モデルの学習のためにデータオーグメンテーションでホワイトノイズを載せるなど、スペクトログラムで音声を編集してから、位相を推定して音声を復元することで、編集後の音声を聞いて確認することができます。


# まとめ

- 音声処理ライブラリLibROSAを使って、音声の読み込み　⇒　スペクトログラムへの変換　⇒　音声の復元をしました。

[^1]: Vibe Ace ([Kevin MacLeod](http://freemusicarchive.org/music/Kevin_MacLeod/)) / [CC BY 3.0](http://creativecommons.org/licenses/by/3.0/)