CUDAで連続ウェーブレット変換
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
*CUDAで連続ウェーブレット変換 [#rc110058]
計算時間のかかる連続ウェーブレット変換をCUDAで実装してみ...
(CUDAのSDKに離散ウェーブレット変換のサンプルあり)．
**連続ウェーブレット変換 [#a14eee32]
ウェーブレット変換(wavelet transform)はマザーウェーブレッ...
元のデータf(t)を時間(もしくは位置)と周波数(スケーリング)...
#ref(eq_wt.jpg,nolink)
マザーウェーブレット(もしくはウェーブレット関数)の代表的...
-Haar wavelet
#ref(eq_haar.jpg,nolink)
-Mexican hat wavelet
#ref(eq_mexicanhat.jpg,nolink)
などがあります．
**CUDAでの実装 [#p1287064]
1次元の信号f(t)を入力として，そのウェーブレット変換W(a,b)...
入力データfは離散化されたデータf_iとして，配列に格納され...
また，今回は実数データのみで検証しますが，実装は虚数値に...
ウェーブレット関数にはMexican hatを用います．
***メモリ確保とデータ転送 [#ib8ba7ad]
ホストメモリ上の入力データをデバイスメモリに転送します．
#code(C){{
float *dSim, *dSre, *dWre, *dWim;
int size;
// デバイスメモリの確保とホストからの転送
size = ntrans*sizeof(float);
cutilSafeCall(cudaMalloc((void**)&dSre, size));
cutilSafeCall(cudaMalloc((void**)&dSim, size));
cutilSafeCall(cudaMemcpy((void*)dSre, (void*)hFi, size, ...
cutilSafeCall(cudaMemset((void*)dSim, 0, size));
size = ntrans*nscale*sizeof(float);
cutilSafeCall(cudaMalloc((void**)&dWre, size));
cutilSafeCall(cudaMalloc((void**)&dWim, size));
cutilSafeCall(cudaMemset((void*)dWre, 0, size));
cutilSafeCall(cudaMemset((void*)dWim, 0, size));
}}
入力データのサイズをntrans, 周波数方向の変換解像度をnscal...
dSreとdSimがそれぞれ入力データの実部と虚部，dWreとdWimに...
***カーネル [#yeb40233]
変換カーネルは，
#code(C){{
__device__
float MexicanHat(float t)
{
t = t*t;
return MEXICAN_HAT_C*(1.0-t)*exp(-t/2.0);
}
__device__
float MexicanHatIm(float t)
{
return 0.0f;
}
__global__
void cwt(float *src_re, float *src_im, float *wt_re, floa...
int nt, int ns, float t0, float s0, float dt, float ds...
{
int i = blockIdx.x*blockDim.x+threadIdx.x;
int j = blockIdx.y*blockDim.y+threadIdx.y;
if(i < nt && j < ns){
float s = (s0+j*ds); // スケーリング(周波数)
if(s == 0.0f) s = 1e-10;
float t = (t0+i*dt); // 平行移動(位置)
float wr = 0.0;
float wi = 0.0;
float kstep = 1.0/res;
float k;
for(k = 0.0; k < nt; k += kstep){
float T = (k*dt-t)/s;
wr += src_re[(uint)k]*MexicanHat(T)+src_im[(uint)k]*Me...
wi += src_re[(uint)k]*MexicanHat(T)-src_im[(uint)k]*Me...
}
wt_re[j*nt+i] = wr/(sqrt(s)*ivalp);
wt_im[j*nt+i] = wi/(sqrt(s)*ivalp);
}
}
}}
引数として，入力データ，出力先，時間(位置)方向と周波数方...
resはウェーブレット関数合成時の分解数です．
そのまんま実装したので，シェアードメモリやテクスチャメモ...
カーネル呼び出し側は，
#code(C){{
dim3 block(BLOCK_SIZE, BLOCK_SIZE);
dim3 grid((ntrans+block.x-1)/block.x, (nscale+block.y-1)...
cwt<<< grid, block >>>(dSre, dSim, dWre, dWim, ntrans, n...
// カーネル実行エラーのチェック
cutilCheckMsg("Kernel execution failed");
// GPUスレッドが終わるのを待つ
cutilSafeCall(cudaThreadSynchronize());
}}
***ホストへの結果の転送と後片付け [#c2407ecd]
結果をホストに転送し，確保したメモリを解放します．
#code(C){{
// デバイスからホストへ結果を転送
size = ntrans*nscale*sizeof(float);
cutilSafeCall(cudaMemcpy(hWre, dWre, size, cudaMemcpyDev...
cutilSafeCall(cudaMemcpy(hWim, dWim, size, cudaMemcpyDev...
// デバイスメモリ解放
cutilSafeCall(cudaFree(dSre));
cutilSafeCall(cudaFree(dSim));
cutilSafeCall(cudaFree(dWre));
cutilSafeCall(cudaFree(dWim));
}}
***結果 [#q165db72]
CPU用にも実装し(OpenMPで2スレッド並列で実行)，計算時間を...
CPU側の実装は計算の最適化として，ウェーブレット関数値の前...
入力データとして，周波数が異なるsin波の合成波を用い，
データ数は時間方向に256,512,1024，周波数方向に128で計測し...
||256|512|1024|h
|CPU|102|300|694|
|GPU|11|18|42|
単位はミリ秒です．デバイスメモリの確保，転送を含めるとGPU...
気になったので調べてみると，最初のcudaMallocで40-80msぐら...
その後のcudaMemcpyやcudaMallocは1msぐらいなのでだいぶ遅い...
CUDA Programming Guideの3.2章の最初の方に
There is no explicit initialization function for the run...
とあるように，最初にランタイム関数が呼ばれたときに初期化...
今回は1回のみの計測でしたが，何回も計測して平均を取るとき...
変換結果を以下に示します．
CENTER:&ref(wavelet_trans.jpg,nolink);
終了行:
*CUDAで連続ウェーブレット変換 [#rc110058]
計算時間のかかる連続ウェーブレット変換をCUDAで実装してみ...
(CUDAのSDKに離散ウェーブレット変換のサンプルあり)．
**連続ウェーブレット変換 [#a14eee32]
ウェーブレット変換(wavelet transform)はマザーウェーブレッ...
元のデータf(t)を時間(もしくは位置)と周波数(スケーリング)...
#ref(eq_wt.jpg,nolink)
マザーウェーブレット(もしくはウェーブレット関数)の代表的...
-Haar wavelet
#ref(eq_haar.jpg,nolink)
-Mexican hat wavelet
#ref(eq_mexicanhat.jpg,nolink)
などがあります．
**CUDAでの実装 [#p1287064]
1次元の信号f(t)を入力として，そのウェーブレット変換W(a,b)...
入力データfは離散化されたデータf_iとして，配列に格納され...
また，今回は実数データのみで検証しますが，実装は虚数値に...
ウェーブレット関数にはMexican hatを用います．
***メモリ確保とデータ転送 [#ib8ba7ad]
ホストメモリ上の入力データをデバイスメモリに転送します．
#code(C){{
float *dSim, *dSre, *dWre, *dWim;
int size;
// デバイスメモリの確保とホストからの転送
size = ntrans*sizeof(float);
cutilSafeCall(cudaMalloc((void**)&dSre, size));
cutilSafeCall(cudaMalloc((void**)&dSim, size));
cutilSafeCall(cudaMemcpy((void*)dSre, (void*)hFi, size, ...
cutilSafeCall(cudaMemset((void*)dSim, 0, size));
size = ntrans*nscale*sizeof(float);
cutilSafeCall(cudaMalloc((void**)&dWre, size));
cutilSafeCall(cudaMalloc((void**)&dWim, size));
cutilSafeCall(cudaMemset((void*)dWre, 0, size));
cutilSafeCall(cudaMemset((void*)dWim, 0, size));
}}
入力データのサイズをntrans, 周波数方向の変換解像度をnscal...
dSreとdSimがそれぞれ入力データの実部と虚部，dWreとdWimに...
***カーネル [#yeb40233]
変換カーネルは，
#code(C){{
__device__
float MexicanHat(float t)
{
t = t*t;
return MEXICAN_HAT_C*(1.0-t)*exp(-t/2.0);
}
__device__
float MexicanHatIm(float t)
{
return 0.0f;
}
__global__
void cwt(float *src_re, float *src_im, float *wt_re, floa...
int nt, int ns, float t0, float s0, float dt, float ds...
{
int i = blockIdx.x*blockDim.x+threadIdx.x;
int j = blockIdx.y*blockDim.y+threadIdx.y;
if(i < nt && j < ns){
float s = (s0+j*ds); // スケーリング(周波数)
if(s == 0.0f) s = 1e-10;
float t = (t0+i*dt); // 平行移動(位置)
float wr = 0.0;
float wi = 0.0;
float kstep = 1.0/res;
float k;
for(k = 0.0; k < nt; k += kstep){
float T = (k*dt-t)/s;
wr += src_re[(uint)k]*MexicanHat(T)+src_im[(uint)k]*Me...
wi += src_re[(uint)k]*MexicanHat(T)-src_im[(uint)k]*Me...
}
wt_re[j*nt+i] = wr/(sqrt(s)*ivalp);
wt_im[j*nt+i] = wi/(sqrt(s)*ivalp);
}
}
}}
引数として，入力データ，出力先，時間(位置)方向と周波数方...
resはウェーブレット関数合成時の分解数です．
そのまんま実装したので，シェアードメモリやテクスチャメモ...
カーネル呼び出し側は，
#code(C){{
dim3 block(BLOCK_SIZE, BLOCK_SIZE);
dim3 grid((ntrans+block.x-1)/block.x, (nscale+block.y-1)...
cwt<<< grid, block >>>(dSre, dSim, dWre, dWim, ntrans, n...
// カーネル実行エラーのチェック
cutilCheckMsg("Kernel execution failed");
// GPUスレッドが終わるのを待つ
cutilSafeCall(cudaThreadSynchronize());
}}
***ホストへの結果の転送と後片付け [#c2407ecd]
結果をホストに転送し，確保したメモリを解放します．
#code(C){{
// デバイスからホストへ結果を転送
size = ntrans*nscale*sizeof(float);
cutilSafeCall(cudaMemcpy(hWre, dWre, size, cudaMemcpyDev...
cutilSafeCall(cudaMemcpy(hWim, dWim, size, cudaMemcpyDev...
// デバイスメモリ解放
cutilSafeCall(cudaFree(dSre));
cutilSafeCall(cudaFree(dSim));
cutilSafeCall(cudaFree(dWre));
cutilSafeCall(cudaFree(dWim));
}}
***結果 [#q165db72]
CPU用にも実装し(OpenMPで2スレッド並列で実行)，計算時間を...
CPU側の実装は計算の最適化として，ウェーブレット関数値の前...
入力データとして，周波数が異なるsin波の合成波を用い，
データ数は時間方向に256,512,1024，周波数方向に128で計測し...
||256|512|1024|h
|CPU|102|300|694|
|GPU|11|18|42|
単位はミリ秒です．デバイスメモリの確保，転送を含めるとGPU...
気になったので調べてみると，最初のcudaMallocで40-80msぐら...
その後のcudaMemcpyやcudaMallocは1msぐらいなのでだいぶ遅い...
CUDA Programming Guideの3.2章の最初の方に
There is no explicit initialization function for the run...
とあるように，最初にランタイム関数が呼ばれたときに初期化...
今回は1回のみの計測でしたが，何回も計測して平均を取るとき...
変換結果を以下に示します．
CENTER:&ref(wavelet_trans.jpg,nolink);
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
