Pythonとresample_polyでオーバーサンプリング

scipy.signal.resample_polyを使ってオーバーサンプリングを実装してみました。

resample_polyのデフォルトの設定では、パラメタβ=5.0のカイザー窓によるFIRローパスフィルタが使われています。

注意すべき箇所として、フィルタ長はresample_polyの引数では設定できません。仮に、アップ／ダウンサンプリング周波数の比率up/down の大きい値のほうを n_over、FIRフィルタ長をwin_lengthとすると、

win_length = n_over * 2 * 10 + 1

として自動的に決定しています。つまり、4倍オーバーサンプリングであればフィルタ長は81サンプルとなります。

以下に実際のコードを引用します。

scipy/signaltools.py at v1.0.0 · scipy/scipy · GitHub

# Design a linear-phase low-pass FIR filter
max_rate = max(up, down)
f_c = 1. / max_rate  # cutoff of FIR filter (rel. to Nyquist)
half_len = 10 * max_rate  # reasonable cutoff for our sinc-like function
h = firwin(2 * half_len + 1, f_c, window=window)

また、カットオフ周波数は元信号のナイキスト周波数がセットされますが、実用的にはもう少し低い周波数に設定したいことが多々あります（参考）。従って、事前にscipy.signal.firwinで任意のフィルタ長でフィルタ係数を生成し、resample_polyの引数windowに渡す方が良さそうです。

以下に、サンプリング周波数 48 kHzの白色雑音に対して、カイザー窓の形状パラメタとタップ数を変えてオーバーサンプリングを実行した例を示します。実行後、元信号との誤差を周波数領域で可視化しました。

オーバーサンプリング処理はそれなりに演算量がかかるためフィルタ長を短くしたいですが、トレードオフとしてフィルタ長を短くすると誤差が大きくなります。その比較の参考になれば幸いです。

カイザー窓のパラメタ（beta）を変えた場合

FIRフィルタ長を変えた場合

カイザー窓のパラメタは5.0に固定しています。

ソースコード

ライセンス：MIT

クリックで展開

import soundfile as sf
import numpy as np
import scipy.signal
import librosa
import matplotlib.pyplot as plt


filepath = "wn.wav"
x, sr = sf.read(filepath)
n_over = 4
sr_over = sr * 4

f_offset = 500.0
fc = (sr/2.0 - f_offset) / (sr_over/2.0) 


windows = [
    scipy.signal.firwin(1+2*4*n_over, fc, window=("kaiser", 0.5)),
    scipy.signal.firwin(1+2*4*n_over, fc, window=("kaiser", 1.0)),
    scipy.signal.firwin(1+2*4*n_over, fc, window=("kaiser", 2.0)),
    scipy.signal.firwin(1+2*4*n_over, fc, window=("kaiser", 4.0)),
    scipy.signal.firwin(1+2*4*n_over, fc, window=("kaiser", 5.0)),
    scipy.signal.firwin(1+2*4*n_over, fc, window=("kaiser", 8.0)),
    scipy.signal.firwin(1+2*4*n_over, fc, window=("kaiser", 16.0)),
]

windows_name = [
    "('kaiser',  0.5)          ",
    "('kaiser',  1.0)          ",
    "('kaiser',  2.0)          ",
    "('kaiser',  4.0)          ",
    "('kaiser',  5.0) [Default]",
    "('kaiser',  8.0)          ",
    "('kaiser', 16.0)          ",
]


windows = [
    scipy.signal.firwin(511, fc, window=("kaiser", 5.0)),
    scipy.signal.firwin(257, fc, window=("kaiser", 5.0)),
    scipy.signal.firwin(129, fc, window=("kaiser", 5.0)),
    scipy.signal.firwin(65, fc, window=("kaiser", 5.0)),
    scipy.signal.firwin(33, fc, window=("kaiser", 5.0)),
    scipy.signal.firwin(17, fc, window=("kaiser", 5.0)),
    scipy.signal.firwin(9, fc, window=("kaiser", 5.0)),
]

windows_name = [
    "FIR 511 taps             ",
    "FIR 257 taps             ",
    "FIR 129 taps             ",
    "FIR  65 taps             ",
    "FIR  33 taps             ",
    "FIR  17 taps             ",
    "FIR   9 taps             ",
]



print(f"sr:      {sr}")
print(f"sr_over: {sr_over}")
errors = []
n_fft = 2048
n_bins = n_fft//2+1
freqs = np.arange(0,n_bins) / n_bins * (sr/2)

for k in range(len(windows)):
    y     = scipy.signal.resample_poly(x, n_over, 1, window=windows[k])
    # y = your_processing_with oversampled(y)
    x_rec = scipy.signal.resample_poly(y, 1, n_over, window=windows[k])
    
    err   = np.mean((x-x_rec)**2) 
    diff = np.abs(librosa.stft(x, n_fft=n_fft) - librosa.stft(x_rec, n_fft=n_fft))
    err_f = np.mean(diff**2, axis=1)
    errors.append({
        "name":  windows_name[k],
        "err":    f"{err:0.5}"
    })
    
    plt.plot(freqs, err_f, label=windows_name[k])
plt.legend()
plt.xscale("log")
plt.yscale("log")
plt.xlabel("Freq.")
plt.ylabel("Error")
plt.show()

print(errors)