現在音楽ファイルで最もよく入手できるフォーマットはサンプリングレート44.1kspsか48ksps、ビット幅は16bitです。フリーソフトのAudacityを使用すれば、サンプリング変換やbit幅などのフォーマット変換を簡単にやってくれます。
どの程度の精度で変換できるかわからないので、自分でfirフィルタを作って検討したところ、計算に時間がかかりすぎ、それを早くしようとPCに最初から搭載されているゲーム用のGPUを使って高速化しようとして、泥沼にはまってしまいました。格闘の末「 Out of memory error 」を回避して、無事に高速化できました。
大規模な畳み込み演算を検討されている方、同じ悩みをかかえていると思ったので、具体的なプログラムも公開します。（重要なポイントのみ）
音楽に色々とフィルタリングをかけてどのような音になるかを試したくなり、自由にフィルタ設計できるFIRフィルタを検討することにしました。周波数5Hz以上のフィルターを設計するには、48kspsでは48000/5=9600タップのFIRフィルターにしないといけません。また、補間値を計算して２倍のサンプリングにするにはSinc関数を使いますが、ΣSinc(i*(t-0.5))/t(Σはi=±Ntap/2)の値が１．０との誤差として0.006%ぐらいになります。
----
RATE=48000
FT=4800 #tap数はこの倍9600
tc = np.arange(-(FT-0.5)/RATE, (FT)/RATE,1/RATE, dtype='float64')
sc=np.sin(tc*RATE*np.pi)/tc/RATE/np.pi
print(np.sum(sc)
------
ということで、実際の8分09秒の音楽ファイルmp3(44.1ksps16bit2ch)を読みこんで、tap数9600タップでnumpy.convolve(aout,sc)で計算しました。
#音声データをarray.arrayに読み込む
adata=audio_data.get_array_of_samples()
#読み込んだデーターをL,Rに分離
bdata=np.frombuffer(adata, dtype="int16").reshape(-1,2).transpose()/(2**15)
#２倍サンプリング補間の例
ad0 =np.convolve(bdata[0],sc,'same')
ad1 =np.convolve(bdata[1],sc,'same')
実行結果、最後の２行のコンボルーション演算だけで、64秒かかりました。
何とか、早くできないかと調べていたら、NumpyをCupyに置き換えればよいと書いてありました。
import numpy as np
に下記の行を加え、
import cupy as cp
pool = cp.cuda.MemoryPool(cp.cuda.malloc_managed)
cp.cuda.set_allocator(pool.malloc)
np.をcp.の記述に変えるだけ
#２倍サンプリング補間の例
ad0_cp =cp.convolve(cp.asarray(bdata[0]),cp.asarray(sc),'same')
ad1_cp =cp.convolve(cp.asarray(bdata[1]),cp.asarray(sc),'same')
を実行しました。すると、
OutOfMemoryError: Out of memory allocating 172,818,432 bytes (allocated so far: 1,723,058,688 bytes).
メモリエラーが出てしまいました。このノートPCは、CPUがi7-9750H 2.6GHzで、GPUとしてNVIDIA GEFORCE MX-250(cudaコア384)がついています。すんなり高速化できるのにと期待していたのにうまくいきません。
さらに調べると、GPUで使用できるメモリはあまり大きくないので、演算を分断してやる工夫が必要という解説文を見つけました。早速、分断してconvolution演算をしようと思い、アルゴリズムとして少し長めのデータでconvolution演算を行い、フィルター側のタップ数の半分の数の両端のデーターをちょん切ってフィルムのようにつなぎあわす方法を考えました。しかし、プログラムを書こうとするとすごく複雑になり、さっきの事をすぐに忘れてしまう年寄りの脳ではプログラムを書くとバグだらけになりそうです。
プログラムが書けるようにアルゴリズムをわかりやすく図に整理して、それを見ながらコーディング、こうすれば誤りを極端に減らすことができます。
アルゴリズムの方針としては、ブロック毎にconvlove計算を行っては、del でメモリを解放する
計算した結果の両端のデータをちょん切ってうまくつなぎ合わせる。
------------
#iwave0(音声データfloat64 numpy),filrfs(フィルタータップfloat64 numpy)
def cpconvolve(iwave0,filrfs):
#音楽ファイルの長さ
inn=len(iwave0)
#フィルターの半分の長さ
tfr=len(filrfs)//2
#分断するサイズ5M要素(float64ならx8 Bytesのサイズ)
#使用するGPUのメモリ規模に合わせて設定する
itn=5000000
#分断したブロック数
nn=(inn-2*tfr)//(itn-2*tfr)
#GPUメモリに転送
filrfs_cp=cp.asarray(filrfs)
iwn=tfr
if inn<=itn:
iwave_cptmp=cp.asarray(iwave0)
cwaves = cp.asnumpy(cp.convolve(iwave_cptmp,filrfs_cp,'same'))
del iwave_cptmp
elif inn<=itn+tfr:
iwave_cptmp=cp.asarray(iwave0[0:itn])
cwaves = cp.asnumpy(cp.convolve(iwave_cptmp,filrfs_cp,'same'))
del iwave_cptmp
iwave_cptmp=cp.concatenate([cp.asarray(iwave0[itn-2*tfr : inn]), cp.zeros(2*itn-inn-2*tfr)])
cwave = cp.asnumpy(cp.convolve(iwave_cptmp,filrfs_cp,'same'))
del iwave_cptmp
cwaves = np.concatenate([cwaves[:itn-tfr],cwave[tfr : inn-itn+2*tfr]])
else :
iwave_cptmp=cp.asarray(iwave0[0:itn])
cwaves = cp.asnumpy(cp.convolve(iwave_cptmp,filrfs_cp,'same'))[:itn-tfr]
del iwave_cptmp
iwn=iwn + (itn-2*tfr)
while iwn+itn-2*tfr< len(iwave0):
iwave_cptmp=cp.asarray(iwave0[iwn-tfr:iwn-tfr+itn])
cwave = cp.asnumpy(cp.convolve(iwave_cptmp,filrfs_cp,'same'))
cwaves= np.concatenate([cwaves,cwave[tfr: itn-tfr]])
del iwave_cptmp
iwn=iwn + (itn-2*tfr)
if iwn+itn-2*tfr != len(iwave0):
iwave_cptmp=cp.concatenate([cp.asarray(iwave0[iwn-tfr:inn]),cp.zeros(iwn+itn-tfr-inn)])
cwave = cp.asnumpy(cp.convolve(iwave_cptmp,filrfs_cp,'same'))
del iwave_cptmp
cwaves = np.concatenate([cwaves[:iwn],cwave[tfr : inn-iwn+tfr]])
del filrfs_cp
return cwaves
-----------
mainでは下記のように記述
ad0 = cpconvolve(bdata[0],sc)
ad1 = cpconvolve(bdata[1],sc)
これで計算すると、メモリエラーも発生せずに、わずか4秒で計算が終了しました。
64秒/4秒=16倍 高速化成功！
長い音声データーにフィルタをかける場合はこの方法でGPUを使用すると大幅に速度改善できます。
特にタップが1万tap規模のような積和計算のときGPUの能力が発揮できます。
この関数cpconvolve()は全くの自作でまだどこかにバグがあるかもしれませんので、詳しい方、使ってみて是非ともフィードバックお願いいたします。