1 はじめに
インテル® アドバンスト・ベクトル・エクステンション (インテル® AVX) は、インテル® ストリーミング SIMD 拡張命令 (インテル® SSE) の 256 ビット命令セットの拡張で、浮動小数点演算を多用するアプリケーション向けに設計されています。これらの命令は、3D 幾何学、ビデオ処理、イメージ処理、空間 (3D) オーディオなどの、浮動小数点演算に強く依存するアプリケーションを高速化します。
この記事では、 Wiener フィルターについて説明し、インテル® AVX を使用して最適化したコードのサンプルも紹介します。この記事とコードは、参考文献 [4] の「Wiener Filtering Using Streaming SIMD Extensions (AP-807)」で紹介されたものに基づいています。オリジナルの記事にはストリーミング SIMD 拡張命令を使用した最適化が含まれていますが、ここでは 256 ビット拡張命令セット (例えば、インテル® AVX) へのコードの移行について説明します。
2 Wiener フィルターのアルゴリズム
Wiener フィルター (LMS フィルターとも呼ばれる) は、イメージから不要なノイズを除去します。このアルゴリズムの説明は、[1] に記述されています。アルゴリズムの入力は 4 つの (フーリエ変換された) ベクトルで、それぞれ、元のイメージ (Image)、劣化イメージ (Guv)、ノイズ・イメージ・スペクトル (Noise)、劣化補正 (Huv) を表します。各入力は、複素数行列のベクトルです。複素数は、数値の実数部と虚数部を 2 つの浮動小数点数で表現します。追加のパラメーター、ガンマも計算に含まれます。ガンマが 1.0 の場合、フィルターはノンパラメトリックです。フィルターのパラメーターは、フィルター処理されたイメージが最適な状態になるまで調整できます。
2.1 Wiener フィルターを使用するアプリケーション
Wiener フィルターは、再構築されたイメージからノイズを除去する際にイメージ処理アプリケーションで一般的に使用されるフィルターで、特にぼかしイメージの復元に使用されます。Wiener フィルターは適応フィルターにおいて重要で、ウェーブレット変換をはじめ、通信やその他のデジタル信号処理 (DSP) 関連のアプリケーションでも使用されます。信号処理の分野では、フーリエ変換が重要な要素の 1 つです。フーリエ変換の実装についての詳細は、[3] を参照してください。
2.2 Wiener フィルターの実装
前述したように、関数の入力は複素数の 4 つの配列です。イメージの各要素について、複素数演算を使用して以下の演算を実行します。
複素変数 D および Hs は、計算中に使用する中間変数です。関数 Complex_conj は、共役複素数を出力します。
除算を行う場合、(if 文を使用して) 分母がゼロでないことを確認する必要があります。分母がゼロの場合、結果をゼロに設定すべきです。
1. Complex Noise = gamma * (Noise * Complex_conj ( Noise ) )
2. Complex D = Image * Complex_conj ( Image )
3. Complex D = Noise / D
4. Complex Hs = Huv * Complex_conj ( Huv )
5. Complex Num = Complex_conj ( Huv ) * Guv
6. Complex Image = Num / (Hs + D)
3 128 ビット SIMD を使用したコードのベクトル化
この記事で引用したソースコード (C バージョンと 128 ビット SIMD バージョンの両方) は、[4] を参照してください。以下の説明は、最初の 128 ビット SIMD スカラーコードの実装について述べたものです。まず、多くの演算に数値とその共役複素数の乗算が含まれているため、これを最適化します。結果には虚数部コンポーネントがないため、セクション 2.2 の多くの演算が単純化されます。
最適化後の C コードはセクション 5 を参照してください。コードをインテル® SSE 向けに最適化する前に、SIMD 実行に適した形式に変換する必要があります。元の C コードの 4 つの反復をまとめて、新しいループの 1 つの反復で処理します。新しいループの各パスで、元の 4 つの反復が処理されます。セクション 2.2 の演算のリストは、反復ごとに、ステップ 3 で 1 回 (虚数部がゼロで、虚数部コンポーネントを除算する必要がないため)、ステップ 6 で 2 回、計 3 回の除算を行う必要があります。後者の 2 回の除算の 1 つは、どちらも同じ分母であることから削除できます。分母がゼロかどうか確認する if 文は、マスキングを使用することで削除できます。すべての除算を逆数近似に置換することで、さらにパフォーマンスを向上させることができます。
以下に、使用した手法を示します。コードを SIMD 形式に変換した後、非ゼロ分母要素のマスクを作成し、マスクと除算結果の AND 演算を行うことでゼロ除算が発生する要素を排除して、ゼロ分母の確認 (if 文) を削除します。この手法では、MXCSR レジスターの SIMD 浮動小数点例外をマスクして、SNAN ではなく QNAN が生成されると仮定しています。そのため、浮動小数点のゼロ除算例外は発生しません。例えば、量 ( N / D ) を計算すると仮定します (N と D は浮動小数点数)。一般的なコードシーケンスは以下のようになります。
If ( D != 0 )
Result = N / D;
Else
Result = 0.0;
(and ( div ( N, D ), cmp_neq ( D, 0 ) ) ) として計算されるため、if 文は必要ありません。
組込み関数を使用すると、以下のように表現できます。
mm_and_ps ( mm_div_ps ( N, D ), mm_cmpneq_ps ( D, zero ) )
この手法はベクトル化されたコード特有で、組込み関数とコードのアセンブリー言語バージョンで実装されます。
ニュートン-ラフソン法は近似を求める古典的な手法です。逆数の最初の「推定」は、rcpps 命令を使用して計算されます。続いて、ニュートン-ラフソン法を使用して「推定」が改善されます。結果は divide 命令で計算されるものほど正確ではありませんが、非常に高速に結果を得ることができます。プログラマーは、アプリケーションで精度の低い解を許可するかどうか判断する必要があります。この手法の詳細および説明は、[2] を参照してください。以下に、この Wiener フィルターで使用したコードシーケンスを示します。ゼロ除算が発生しないように分母を確認する必要があります。
RC = _mm_rcpps( D );
RECIP = _mm_sub( _mm_add( RC, RC ), _mm_mul( RC, _mm_mul( RC, D ) ) );
4 256 ビット SIMD を使用したコードのベクトル化
128 ビットでベクトル化されたコードは、256 ビット SIMD に簡単に移植できます。256 ビットのコードは 1 つのループで 8 回分の反復を処理します (128 ビットのインテル® SSE では 1 つのループで 4 回分の反復を処理します)。load/store 命令だけでなく、rcp/mul/add/sub 命令も 8 つのデータポイントを処理します。上記で説明した分岐なし手法は、インテル® AVX コードでも使用されています。コードから分かるように、対応するインテル® AVX 組込み関数 (例えば、_mm_mul_ps の代わりに _mm256_mul_ps) を使用するだけで移植は完了です。完全な 256 ビットのインテル® AVX コードは、セクション 6.3 を参照してください。
5 入出力配列のグループ化
入出力配列をシーケンシャルにグループ化することで、256 ビット SIMD コードのパフォーマンスをさらに向上させることができます。キャッシュ/メモリーでシーケンシャルなメモリーアクセスを行うと、潜在的なキャッシュウェイの競合が減ります。CPU ハードウェア・プリフェッチャーに単純なアクセスパターンが提供されることで、データがより正確になります。
6 まとめ
以下に、128 ビット SIMD コードと 256 ビット SIMD コードで大量の反復を処理した場合のパフォーマンス結果 (CPU クロックサイクル) を示します。
|
|
インテル® AVX (256 ビット) |
インテル® SSE (128 ビット) |
インテル® AVX のパフォーマンス向上率 (対インテル® SSE) |
|
Wiener フィルター |
45871 |
66933 |
1.46 倍 |
|
Wiener フィルターとグループ配列 |
42464 |
64473 |
1.51 倍 |
インテル® AVX におけるパフォーマンス向上率は、インテル® SSE の 1.46 倍でした。Wiener フィルターのアルゴリズムにインテル® AVX を使用した場合、128 ビット SIMD コードに比べパフォーマンスが向上することが分かります。入出力配列をグループ化すると、インテル® AVX におけるパフォーマンス向上率は、インテル® SSE の 1.51 倍になりました。
パフォーマンスの向上は、インテル® AVX を使用すること (コードのベクトル化) と、インテル® SSE およびインテル® AVX のマスク操作を使用して条件分岐命令 (if 文) を削除することにより得られたものです。結果の数値の精度が低下してもかまわない場合は、(ニュートン-ラフソン法を使用して) 除算を逆数近似に置換することで、さらにパフォーマンスを向上させることができます。この記事は、既存の浮動小数点コードをインテル® AVX に簡単に移植できることも示しています。
7 コードサンプル
/*
* Wiener Filter (also known as the Least Mean Square filter)
*
* Reference: The Pocket Handbook of Image Processing Algorithms in C
* by Harley R Myler & Arthur R. Weeks
* 1993 Prentice-Hall, ISBN 0-13-642240-3 p260-3.
*
* The data is several arrays of complex floats in row major order.
* The description for the algorithm from p260 states:
*
* The algorithm computes a parametric Wiener filter on the
* Fourier transform of a degraded image, Guv, with noise
* spectra N, degradation function Huv, and original image Img.
* The computation is in place, so that the filtered version of
* the input is returned in the original image variable. The
* original and noise images are either estimations form some
* predictive function or ad hoc approximations. If the noise
* image is zero, the process reduces to the inverse filter.
*
* The Weiner parameter gamma is passed to the algorithm.
* If this parameter is 1.0, the filter is non-parametric.
* Methods exist in the literature to derive the parameter value;
* however, it is sometimes determined from trial and error.
*
*NOTE!!!! The code on page 263 has an error. In cxml, the complex
* multiply routine, the imaginary part of the computation should be
* a*d + b*c, not a*d - b*c.
*
*NOTE! (another error) The *complex* array length is rows*cols, so the
* *float* array length should be 2*rows*cols. Also, note that the
* algorithm operates on one component of the pixel.
*/
void wiener_filter ( float *Img,
float *Huv,
float *No,
float *Guv,
float gamma,
int rows,
int cols)
{
int i, sz;
float numr, numi, dr, hsr;
sz = 2 * rows * cols;
for (i = 0; i < sz; i += 2)
{
/* Compute (in place) the noise spectral density with Wiener gamma*/
No[i] = (float) ( gamma * ( No[i]*No[i] + No[i+1]*No[i+1] ) );
No[i+1] = (float) 0.0;
/* Compute image spectral density */
dr = (float) ( Img[i]*Img[i] + Img[i+1]*Img[i+1] );
/* Compute denominator spectral density term */
if (dr != 0.0)
dr = (float) (No[i] / dr) ;
/* Compute degradation power spectrum */
hsr = (float) ( Huv[i]*Huv[i] + Huv[i+1]*Huv[i+1] );
/* Compute numerator term */
numr = (float) ( Huv[i]*Guv[i] + Huv[i+1]*Guv[i+1] );
numi = (float) ( Huv[i]*Guv[i+1] - Huv[i+1]*Guv[i ] );
/* Final computation */
if ( (hsr + dr) != 0.0 )
{
Img[i] = (float) (numr / (hsr + dr));
Img[i+1] = (float) (numi / (hsr + dr));
}
else
{
Img[i] = (float) 0.0;
Img[i+1] = (float) 0.0;
}
}
} /* wiener_filter */
7.2 128 ビット組込み関数コード
#include <iostream.h>
//#define MM_FUNCTIONALITY
#include <xmmintrin.h>
#include <assert.h>
void intrin_wiener_rcp_sse( float *Img,
float *Huv,
float *No,
float *Guv,
float gamma,
int rows,
int cols )
{
int i, sz;
__m128 first2, next2, nor4, noi4, nr4, inr4, ini4, dr4;
__m128 hr4, hi4, hsr4, gr4, gi4, numr4, numi4;
__m128 rc, denom;
__m128 zero = _mm_set_ps1 (0.0);
sz = 2 * rows * cols;
assert( (sz > 3) & !(sz & 3) );
assert( !( ((int)Img) & 15 ) ); /* Assume alignment */
assert( !( ((int)Huv) & 15 ) );
assert( !( ((int)No) & 15 ) );
assert( !( ((int)Guv) & 15 ) );
for (i = 0; i < sz; i += 8)
{
* Compute (in place) the noise spectral density with Wiener gamma
*
* complex Noise = gamma * (Noise * complex conj Noise)
*
* No[i] = (float) ( gamma * ( No[i]*No[i] + No[i+1]*No[i+1] ) );
* No[i+1] = (float) 0.0;
*/
first2 = _mm_load_ps ( &No[i] );
next2 = _mm_load_ps ( &No[i+4] );
nor4 = _mm_shuffle_ps( first2, next2, 0x88 );
noi4 = _mm_shuffle_ps( first2, next2, 0xdd );
nr4 = _mm_mul_ps ( _mm_set_ps1( gamma ) ,
_mm_add_ps ( _mm_mul_ps( nor4 , nor4 ),
_mm_mul_ps( noi4 , noi4 ) ) );
_mm_store_ps( &No[i ], _mm_unpacklo_ps ( nr4, zero ) );
_mm_store_ps( &No[i+4], _mm_unpackhi_ps ( nr4, zero ) );
/*
* Compute image spectral density
*
* Complex D = Image * complex conj Image
*
* dr = (float) ( Img[i]*Img[i] + Img[i+1]*Img[i+1] );
*/
first2 = _mm_load_ps ( &Img[i] );
next2 = _mm_load_ps ( &Img[i+4] );
inr4 = _mm_shuffle_ps( first2, next2, 0x88 );
ini4 = _mm_shuffle_ps( first2, next2, 0xdd );
dr4 = _mm_add_ps ( _mm_mul_ps( inr4 , inr4),
_mm_mul_ps( ini4 , ini4) );
/*
* Compute denominator spectral density term
*
* Complex D = noise / D
*
* if (dr != 0.0)
* dr = (float) (No[i] / dr) ;
*
* Do that reciprical division thing!
*/
rc = _mm_rcp_ps(dr4);
rc = _mm_sub_ps( _mm_add_ps( rc, rc),
_mm_mul_ps( rc, _mm_mul_ps( rc, dr4) ) );
dr4 = _mm_and_ps ( _mm_mul_ps ( nr4 , rc ),
_mm_cmpneq_ps( dr4, zero ) );
/*
* Compute degradation power spectrum
*
* Complex Hs = Huv * complex conj Huv
*
* hsr = (float) ( Huv[i]*Huv[i] + Huv[i+1]*Huv[i+1] );
*/
first2 = _mm_load_ps ( &Huv[i] );
next2 = _mm_load_ps ( &Huv[i+4] );
hr4 = _mm_shuffle_ps( first2, next2, 0x88 );
hi4 = _mm_shuffle_ps( first2, next2, 0xdd );
hsr4 = _mm_add_ps ( _mm_mul_ps (hr4 , hr4 ),
_mm_mul_ps (hi4 , hi4 ) );
/*
* Compute numerator term
*
* Complex Num = complex conj Huv * Guv
*
* numr = (float) ( Huv[i]*Guv[i] + Huv[i+1]*Guv[i+1] );
* numi = (float) ( Huv[i]*Guv[i+1] - Huv[i+1]*Guv[i ] );
*/
first2 = _mm_load_ps ( &Guv[i] );
next2 = _mm_load_ps ( &Guv[i+4] );
gr4 = _mm_shuffle_ps( first2, next2, 0x88 );
gi4 = _mm_shuffle_ps( first2, next2, 0xdd );
numr4 = _mm_add_ps ( _mm_mul_ps (hr4 , gr4),
_mm_mul_ps (hi4 , gi4) );
numi4 = _mm_sub_ps ( _mm_mul_ps (hr4 , gi4),
_mm_mul_ps (hi4 , gr4) );
/*
* Final computation
*
* Complex Image = Num / (Hs + D)
*
* if ( (hsr + dr) != 0.0 )
* {
* Img[i] = (float) (numr / (hsr + dr));
* Img[i+1] = (float) (numi / (hsr + dr));
* }
* else
* {
* Img[i] = (float) 0.0;
* Img[i+1] = (float) 0.0;
* }
*
* Do the reciprical division thing
*/
denom = _mm_add_ps( hsr4, dr4 );
rc = _mm_rcp_ps(denom);
rc = _mm_sub_ps( _mm_add_ps( rc, rc),
_mm_mul_ps( rc, _mm_mul_ps( rc, denom) ) );
inr4 = _mm_and_ps( _mm_mul_ps ( numr4 , rc ) ,
_mm_cmpneq_ps( denom, zero ) );
ini4 = _mm_and_ps( _mm_mul_ps ( numi4 , rc ) ,
_mm_cmpneq_ps( denom, zero ) );
_mm_store_ps( &Img[i ], _mm_unpacklo_ps ( inr4, ini4 ) );
_mm_store_ps( &Img[i+4], _mm_unpackhi_ps ( inr4, ini4 ) );
}
} /* intrin_wiener_rcp */
7.3 256 ビット組込み関数コード
void intrin_wiener_rcp_avx( float *Img,
float *Huv,
float *No,
float *Guv,
float gamma,
int rows,
int cols )
{
int i, sz;
__m256 first2, next2, nor4, noi4, nr4, inr4, ini4, dr4;
__m256 hr4, hi4, hsr4, gr4, gi4, numr4, numi4;
__m256 rc, denom;
__m256 zero = _mm256_setzero_ps();
sz = 2 * rows * cols;
assert( (sz > 3) & !(sz & 3) );
assert( !( ((int)Img) & 15 ) ); /* Assume alignment */
assert( !( ((int)Huv) & 15 ) );
assert( !( ((int)No) & 15 ) );
assert( !( ((int)Guv) & 15 ) );
for (i = 0; i < sz; i += 16)
{
* Compute (in place) the noise spectral density with Wiener gamma
*
* complex Noise = gamma * (Noise * complex conj Noise)
*
* No[i] = (float) ( gamma * ( No[i]*No[i] + No[i+1]*No[i+1] ) );
* No[i+1] = (float) 0.0;
*/
first2 = _mm256_load_ps ( &No[i] );
next2 = _mm256_load_ps ( &No[i+4*2] );
nor4 = _mm256_shuffle_ps( first2, next2, 0x88 );
noi4 = _mm256_shuffle_ps( first2, next2, 0xdd );
nr4 = _mm256_mul_ps ( _mm256_set1_ps( gamma ) ,
_mm256_add_ps ( _mm256_mul_ps( nor4 , nor4 ),
_mm256_mul_ps( noi4 , noi4 ) ) );
_mm256_store_ps( &No[i ], _mm256_unpacklo_ps ( nr4, zero ) );
_mm256_store_ps( &No[i+4*2], _mm256_unpackhi_ps ( nr4, zero ) );
/*
* Compute image spectral density
*
* Complex D = Image * complex conj Image
*
* dr = (float) ( Img[i]*Img[i] + Img[i+1]*Img[i+1] );
*/
first2 = _mm256_load_ps ( &Img[i] );
next2 = _mm256_load_ps ( &Img[i+4*2] );
inr4 = _mm256_shuffle_ps( first2, next2, 0x88 );
ini4 = _mm256_shuffle_ps( first2, next2, 0xdd );
dr4 = _mm256_add_ps ( _mm256_mul_ps( inr4 , inr4),
_mm256_mul_ps( ini4 , ini4) );
/*
* Compute denominator spectral density term
*
* Complex D = noise / D
*
* if (dr != 0.0)
* dr = (float) (No[i] / dr) ;
*
* Do that reciprical division thing!
*/
rc = _mm256_rcp_ps(dr4);
rc = _mm256_sub_ps( _mm256_add_ps( rc, rc),
_mm256_mul_ps( rc, _mm256_mul_ps( rc, dr4) ) );
dr4 = _mm256_and_ps ( _mm256_mul_ps ( nr4 , rc ),
_mm256_cmpneq_ps( dr4, zero ) );
/*
* Compute degradation power spectrum
*
* Complex Hs = Huv * complex conj Huv
*
* hsr = (float) ( Huv[i]*Huv[i] + Huv[i+1]*Huv[i+1] );
*/
first2 = _mm256_load_ps ( &Huv[i] );
next2 = _mm256_load_ps ( &Huv[i+4*2] );
hr4 = _mm256_shuffle_ps( first2, next2, 0x88 );
hi4 = _mm256_shuffle_ps( first2, next2, 0xdd );
hsr4 = _mm256_add_ps ( _mm256_mul_ps (hr4 , hr4 ),
_mm256_mul_ps (hi4 , hi4 ) );
/*
* Compute numerator term
*
* Complex Num = complex conj Huv * Guv
*
* numr = (float) ( Huv[i]*Guv[i] + Huv[i+1]*Guv[i+1] );
* numi = (float) ( Huv[i]*Guv[i+1] - Huv[i+1]*Guv[i ] );
*/
first2 = _mm256_load_ps ( &Guv[i] );
next2 = _mm256_load_ps ( &Guv[i+4*2] );
gr4 = _mm256_shuffle_ps( first2, next2, 0x88 );
gi4 = _mm256_shuffle_ps( first2, next2, 0xdd );
numr4 = _mm256_add_ps ( _mm256_mul_ps (hr4 , gr4),
_mm256_mul_ps (hi4 , gi4) );
numi4 = _mm256_sub_ps ( _mm256_mul_ps (hr4 , gi4),
_mm256_mul_ps (hi4 , gr4) );
/*
* Final computation
*
* Complex Image = Num / (Hs + D)
*
* if ( (hsr + dr) != 0.0 )
AP-807 Wiener Filtering Using Streaming SIMD Extensions
01/28/99 15
* {
* Img[i] = (float) (numr / (hsr + dr));
* Img[i+1] = (float) (numi / (hsr + dr));
* }
* else
* {
* Img[i] = (float) 0.0;
* Img[i+1] = (float) 0.0;
* }
*
* Do the reciprical division thing
*/
denom = _mm256_add_ps( hsr4, dr4 );
rc = _mm256_rcp_ps(denom);
rc = _mm256_sub_ps( _mm256_add_ps( rc, rc),
_mm256_mul_ps( rc, _mm256_mul_ps( rc, denom) ) );
inr4 = _mm256_and_ps( _mm256_mul_ps ( numr4 , rc ) ,
_mm256_cmpneq_ps( denom, zero ) );
ini4 = _mm256_and_ps( _mm256_mul_ps ( numi4 , rc ) ,
_mm256_cmpneq_ps( denom, zero ) );
_mm256_store_ps( &Img[i ], _mm256_unpacklo_ps ( inr4, ini4 ) );
_mm256_store_ps( &Img[i+4*2], _mm256_unpackhi_ps ( inr4, ini4 ) );
}
} /* intrin_wiener_rcp */
7.4 256 ビット組込み関数コードとグループ配列
blockHNG 構造
| Huv[0] | … | Huv[15] | No[0] | … | No[15] | Guv[0] | … | Guv[15] | Huv[16] | … | Huv[31] | No[16] | … |
void intrin_wiener_rcp_avx ( float *Img,
float *_blockHNG,
float gamma,
int rows,
int cols)
{
int sz;
__m256 first2, next2, nor4, noi4, nr4, inr4, ini4, dr4;
__m256 hr4, hi4, hsr4, gr4, gi4, numr4, numi4;
__m256 rc, denom;
__m256 zero = _mm256_setzero_ps();
sz = 2 * rows * cols;
assert( (sz > 3) & !(sz & 3) );
assert( !( ((int)Img) & 15 ) ); // Assume alignment
assert( !( ((int)_blockHNG) & 15 ) ); // Assume alignment
float *Huv;
float *No;
float *Guv;
int j = 0; // img index
for (int _blockHNG_tracker = 0; _blockHNG_tracker < 2 * rows * cols * 3; _blockHNG_tracker += 48)
{
Huv = &(_blockHNG[_blockHNG_tracker]);
No = &(_blockHNG[_blockHNG_tracker + 16]);
Guv = &(_blockHNG[_blockHNG_tracker + 32]);
/*
* Compute (in place) the noise spectral density with Wiener gamma
*
* complex Noise = gamma * (Noise * complex conj Noise)
*
AP-807 Wiener Filtering Using Streaming SIMD Extensions
01/28/99 13
* No[i] = (float) ( gamma * ( No[i]*No[i] + No[i+1]*No[i+1] ) );
* No[i+1] = (float) 0.0;
*/
first2 = _mm256_load_ps ( &No[0] );
next2 = _mm256_load_ps ( &No[8] );
nor4 = _mm256_shuffle_ps( first2, next2, 0x88 );
noi4 = _mm256_shuffle_ps( first2, next2, 0xdd );
nr4 = _mm256_mul_ps ( _mm256_set1_ps( gamma ) ,
_mm256_add_ps ( _mm256_mul_ps( nor4 , nor4 ),
_mm256_mul_ps( noi4 , noi4 ) ) );
_mm256_store_ps( &No[0], _mm256_unpacklo_ps ( nr4, zero ) );
_mm256_store_ps( &No[8], _mm256_unpackhi_ps ( nr4, zero ) );
/*
* Compute image spectral density
*
* Complex D = Image * complex conj Image
*
* dr = (float) ( Img[i]*Img[i] + Img[i+1]*Img[i+1] );
*/
first2 = _mm256_load_ps ( &Img[j] );
next2 = _mm256_load_ps ( &Img[j+8] );
inr4 = _mm256_shuffle_ps( first2, next2, 0x88 );
ini4 = _mm256_shuffle_ps( first2, next2, 0xdd );
dr4 = _mm256_add_ps ( _mm256_mul_ps( inr4 , inr4),
_mm256_mul_ps( ini4 , ini4) );
/*
* Compute denominator spectral density term
*
* Complex D = noise / D
*
* if (dr != 0.0)
* dr = (float) (No[i] / dr) ;
*
* Do that reciprical division thing!
*/
rc = _mm256_rcp_ps(dr4);
rc = _mm256_sub_ps( _mm256_add_ps( rc, rc),
_mm256_mul_ps( rc, _mm256_mul_ps( rc, dr4) ) );
dr4 = _mm256_and_ps ( _mm256_mul_ps ( nr4 , rc ),
_mm256_cmpneq_ps( dr4, zero ) );
/*
* Compute degradation power spectrum
*
* Complex Hs = Huv * complex conj Huv
*
* hsr = (float) ( Huv[i]*Huv[i] + Huv[i+1]*Huv[i+1] );
*/
first2 = _mm256_load_ps ( &Huv[0] );
next2 = _mm256_load_ps ( &Huv[8] );
hr4 = _mm256_shuffle_ps( first2, next2, 0x88 );
hi4 = _mm256_shuffle_ps( first2, next2, 0xdd );
hsr4 = _mm256_add_ps ( _mm256_mul_ps (hr4 , hr4 ),
_mm256_mul_ps (hi4 , hi4 ) );
/*
* Compute numerator term
*
* Complex Num = complex conj Huv * Guv
*
* numr = (float) ( Huv[i]*Guv[i] + Huv[i+1]*Guv[i+1] );
* numi = (float) ( Huv[i]*Guv[i+1] - Huv[i+1]*Guv[i ] );
*/
first2 = _mm256_load_ps ( &Guv[0] );
next2 = _mm256_load_ps ( &Guv[8] );
gr4 = _mm256_shuffle_ps( first2, next2, 0x88 );
gi4 = _mm256_shuffle_ps( first2, next2, 0xdd );
numr4 = _mm256_add_ps ( _mm256_mul_ps (hr4 , gr4),
_mm256_mul_ps (hi4 , gi4) );
numi4 = _mm256_sub_ps ( _mm256_mul_ps (hr4 , gi4),
_mm256_mul_ps (hi4 , gr4) );
/*
* Final computation
*
* Complex Image = Num / (Hs + D)
*
* if ( (hsr + dr) != 0.0 )
AP-807 Wiener Filtering Using Streaming SIMD Extensions
01/28/99 15
* {
* Img[i] = (float) (numr / (hsr + dr));
* Img[i+1] = (float) (numi / (hsr + dr));
* }
* else
* {
* Img[i] = (float) 0.0;
* Img[i+1] = (float) 0.0;
* }
*
* Do the reciprical division thing
*/
denom = _mm256_add_ps( hsr4, dr4 );
rc = _mm256_rcp_ps(denom);
rc = _mm256_sub_ps( _mm256_add_ps( rc, rc),
_mm256_mul_ps( rc, _mm256_mul_ps( rc, denom) ) );
inr4 = _mm256_and_ps( _mm256_mul_ps ( numr4 , rc ) ,
_mm256_cmpneq_ps( denom, zero ) );
ini4 = _mm256_and_ps( _mm256_mul_ps ( numi4 , rc ) ,
_mm256_cmpneq_ps( denom, zero ) );
_mm256_store_ps( &Img[j ], _mm256_unpacklo_ps ( inr4, ini4 ) );
_mm256_store_ps( &Img[j+8], _mm256_unpackhi_ps ( inr4, ini4 ) );
j+=16;
}
} /* Intrin_wiener_rcp_avx */
謝辞
Phil Kerly 氏、Raghu Muthyalampalli 氏、Justin Landon 氏には、コードのパフォーマンス評価、ホワイトペーパーのレビューでご協力いただきました。ここに感謝の意を表します。
参考文献
この記事では、以下の文献を参照しています。これらの文献には、この記事で説明している内容を理解するための背景情報が記述されています。
1. The Pocket Handbook of Image Processing Algorithms in C, by Harley R Myler and Aruthur R.
Weeks. ISBN 0-13-642240-3.
2. Increasing the Accuracy of the Results from the Reciprocal and Reciprocal Square Root
Instructions using the Newton-Raphson Method, Intel application note (AP-803, Order Number: 243637-001).
3. Split-Radix FFT, Intel application note (AP-808, Order Number: 243642-001).
4. Wiener Filtering Using Streaming SIMD Extensions, Intel application note (AP-807).
編集部追加
せっかく労力を掛けて最適化を行っても、全体として性能向上がそれに見合わなければ意味がありません。最適化の対象としてもっとも効果があるのは、アプリケーションでもっともホットな（処理時間が多い）関数や処理です。インテル® VTune™ Amplifier XE を用いることで、多数のソースコードやモジュールの中から Hotspot となる処理を素早く見つけ出すことができます。まずは１度お試しいただき、ご自分のアプリケーションの課題を見つけましょう。評価版のダウンロード、お問い合わせはエクセルソフト株式会社まで。
