NumpyのFFT入門(4) :NumPyを用いた音声信号解析

# NumPyを用いた音声信号解析: FFTによる周波数分析

音声信号解析は、オーディオ処理や音声認識の分野で重要な役割を果たします。この解説では、`NumPy`ライブラリを使用して、音声信号の周波数成分を解析する方法を紹介します。具体的には、音声データに対してフーリエ変換（FFT: Fast Fourier Transform）を適用し、信号の周波数特性を視覚化します。



## 1. FFTによる周波数分析

### コード例

**音源ソース**
https://www.ne.jp/asahi/music/myuu/wave/carenginestart1.wav

```python
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wavfile

# WAVファイルの読み込み
rate, data = wavfile.read("./sound/carenginestart1.wav")

# ステレオ音声の場合は、片方のチャネルを選択
if data.ndim > 1:
    data = data[:, 0]

data = data / 32768.0  # 16ビット音声を-1から1の範囲に正規化

# データの一部を選択
N = min(2048, len(data))  # サンプル数
audio_sample = data[:N]

# FFT実行
fft_result = np.fft.fft(audio_sample)

# 周波数軸の計算
freqs = np.fft.fftfreq(N, 1/rate)

# 音声データとFFT結果のプロット
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(audio_sample)
plt.title('Audio Signal')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.stem(freqs[:N//2], np.abs(fft_result)[:N//2], 'b', markerfmt=" ", basefmt="-b")
plt.title('FFT of Audio Signal')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.grid(True)
plt.tight_layout()
plt.show()

```



音声信号解析は、オーディオ処理や音声認識の分野で重要な役割を果たします。この解説では、`NumPy`ライブラリを使用して、音声信号の周波数成分を解析する方法を紹介します。具体的には、音声データに対してフーリエ変換（FFT: Fast Fourier Transform）を適用し、信号の周波数特性を視覚化します。

### 必要なライブラリ

- `numpy`: 数値計算を効率的に行うためのライブラリです。
- `matplotlib`: データの可視化に用います。
- `scipy.io.wavfile`: WAV形式の音声ファイルを読み込むためのモジュールです。

```python
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wavfile
```

### ステップ1: 音声ファイルの読み込み

まず、`wavfile.read`関数を用いてWAV形式の音声ファイルを読み込みます。音声データとサンプリングレートを取得します。

```python
rate, data = wavfile.read("./sound/carenginestart1.wav")
```

#### チャネル選択

音声データがステレオである場合、片方のチャネルを選択します。これにより、処理を簡略化し、モノラル信号として扱います。

```python
if data.ndim > 1:
    data = data[:, 0]
```

#### 正規化

16ビット音声データを-1から1の範囲に正規化します。これにより、計算が安定しやすくなります。

```python
data = data / 32768.0
```

### ステップ2: フーリエ変換の適用

選択したデータの一部に対してFFTを実行します。ここでは、2048サンプルを使用しますが、音声データの長さが2048サンプル未満の場合には、音声データ全体を使用します。

```python
N = min(2048, len(data))
audio_sample = data[:N]
fft_result = np.fft.fft(audio_sample)
```

### ステップ3: 周波数軸の計算

FFTの結果をプロットするために、対応する周波数軸を計算します。

```python
freqs = np.fft.fftfreq(N, 1/rate)
```

### ステップ4: 結果のプロット

音声信号とそのFFT結果をプロットします。これにより、時間領域と周波数領域での信号の特性を視覚的に確認できます。

#### 音声信号のプロット

まず、音声信号の波形をプロットします。

```python
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(audio_sample)
plt.title('Audio Signal')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.grid(True)
```

#### FFT結果のプロット

FFTの結果を周波数領域でプロットし、信号の周波数成分を確認します。`stem`を使用してスペクトルを示します。

```python
plt.subplot(2, 1, 2)
plt.stem(freqs[:N//2], np.abs(fft_result)[:N//2], 'b', markerfmt=" ", basefmt="-b")
plt.title('FFT of Audio Signal')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.grid(True)
plt.tight_layout()
plt.show()
```



![image-20240807104811443.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/3e2a4806-7920-395f-b491-9842ce8ccf01.png)






## 2. NumPyを用いた音声信号の周波数成分解析とノイズフィルタリング



FFTを用いて得られた周波数成分を分析し、音声信号の特性を理解します。

1. **周波数ピークの特定**: 周波数スペクトルから主要な周波数成分（ピーク）を検出します。この情報は音声の基音や倍音の特定に役立ちます。
2. **ノイズの検出とフィルタリング**: スペクトルからノイズを特定し、信号の質を向上させるためにフィルタリングを行います。



### コード例

```python
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wavfile
from scipy.signal import find_peaks, butter, filtfilt

# WAVファイルの読み込み
rate, data = wavfile.read("./sound/carenginestart1.wav")

# ステレオ音声の場合は、片方のチャネルを選択
if data.ndim > 1:
    data = data[:, 0]

data = data / 32768.0  # 16ビット音声を-1から1の範囲に正規化

# データの一部を選択
N = min(2048, len(data))  # サンプル数
audio_sample = data[:N]

# FFT実行
fft_result = np.fft.fft(audio_sample)

# 周波数軸の計算
freqs = np.fft.fftfreq(N, 1/rate)

# 周波数ピークの特定
magnitude = np.abs(fft_result[:N//2])
peaks, _ = find_peaks(magnitude, height=np.max(magnitude) * 0.1)
peak_freqs = freqs[peaks]
peak_magnitudes = magnitude[peaks]

# ノイズの検出とフィルタリング
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = filtfilt(b, a, data)
    return y

# ここではカットオフ周波数を1000Hzと仮定
filtered_audio = lowpass_filter(audio_sample, 1000, rate)

# プロット
plt.figure(figsize=(12, 8))

# 元の音声信号
plt.subplot(3, 1, 1)
plt.plot(audio_sample)
plt.title('Original Audio Signal')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.grid(True)

# FFT結果とピーク
plt.subplot(3, 1, 2)
plt.stem(freqs[:N//2], magnitude, 'b', markerfmt=" ", basefmt="-b")
plt.plot(peak_freqs, peak_magnitudes, "r*", label='Peaks')
plt.title('FFT of Audio Signal with Peaks')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.legend()
plt.grid(True)

# フィルタリング後の音声信号
plt.subplot(3, 1, 3)
plt.plot(filtered_audio)
plt.title('Filtered Audio Signal (Low-pass)')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.grid(True)

plt.tight_layout()
plt.show()

```



### 周波数ピークの特定

信号のスペクトルから主要な周波数成分（ピーク）を特定します。

```python
magnitude = np.abs(fft_result[:N//2])
peaks, _ = find_peaks(magnitude, height=np.max(magnitude) * 0.1)
peak_freqs = freqs[peaks]
peak_magnitudes = magnitude[peaks]
```

この手順では、ピークの高さが最大値の10%を超える周波数をピークとして検出しています。



### ノイズの検出とフィルタリング

ノイズの影響を軽減するためにローパスフィルタを適用します。ここでは、バターワースフィルタを用いてフィルタリングを行います。

```python
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = filtfilt(b, a, data)
    return y

filtered_audio = lowpass_filter(audio_sample, 1000, rate)
```

フィルタのカットオフ周波数は1000Hzと仮定していますが、信号特性に応じて調整が可能です。



### 結果のプロット

最後に、元の音声信号、FFT結果とピーク、フィルタリング後の音声信号をプロットして可視化します。

```python
plt.figure(figsize=(12, 8))

# 元の音声信号
plt.subplot(3, 1, 1)
plt.plot(audio_sample)
plt.title('Original Audio Signal')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.grid(True)

# FFT結果とピーク
plt.subplot(3, 1, 2)
plt.stem(freqs[:N//2], magnitude, 'b', markerfmt=" ", basefmt="-b")
plt.plot(peak_freqs, peak_magnitudes, "r*", label='Peaks')
plt.title('FFT of Audio Signal with Peaks')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.legend()
plt.grid(True)

# フィルタリング後の音声信号
plt.subplot(3, 1, 3)
plt.plot(filtered_audio)
plt.title('Filtered Audio Signal (Low-pass)')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.grid(True)

plt.tight_layout()
plt.show()
```


![image-20240807104934460.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/23c0435f-4217-b691-ab28-dbabcb5a5206.png)




## まとめ

この解説では、`NumPy`のFFT機能と`matplotlib`を用いて音声信号を解析し、周波数成分を視覚化する方法を紹介しました。この手法は音声信号の特徴抽出やスペクトラム解析において強力なツールです。FFTを活用することで、音声の周波数特性を効率的に分析でき、基音や倍音を特定することが可能です。

さらに、主要な周波数ピークの特定により、音声信号の重要な成分を明確にし、ローパスフィルタを用いたノイズフィルタリングを実施しました。これにより、不要な周波数成分を除去し、信号の質を向上させることができます。これらの技術は音響解析や音声認識の分野で重要な役割を果たし、音響体験をより良いものにするために活用されます。

# 関連資料

https://qiita.com/kotai2003/items/05b4f17948c20c31e15a

https://qiita.com/kotai2003/items/364cb155f25a0646af17

https://qiita.com/kotai2003/items/577c89b20909b989dbf7