Smile Engineering Blog

ジェイエスピーからTipsや技術特集、プロジェクト物語を発信します

適応フィルタを作ってみる

音響処理

信号処理の分野では、適応フィルタが色々な目的で使用されています。ちょっとした理由で適応フィルタのプログラムを作ってみました。

ADF【Adaptive Filter】

適応フィルタ【Adaptive Filter】とは、出力が目的とする信号に近づく（収束する）ように係数を自動で更新していく適応信号処理で、応用例として代表的なものはエコーキャンセラやノイズキャンセラだと思います。

適応フィルタの構成

f:id:jspnet:20200306002903p:plain:right:w400 適応フィルタの構成は、適応アルゴリズムとフィルタ部です。 フィルタ部でd(n)に似たy(n)を合成、

{ 　y(n) = \sum_{k=0}^{N}h_k x(n-k) }

入力x(n)と誤差e(n)に基づいて、フィルタ係数h(n)を最適係数に近づけます。

{ 　e(n)=d(n)-y(n) }

誤差e(n)のパワーを最小にする。

適応アルゴリズム

代表的な適応アルゴリズムとしては、LMS【Least Mean Square】や、NLMS【Normalized Least Mean Square】が思いつきます。その理由は、処理量（演算量）やメモリ構成が現実的なものなので組み込みでも可能になる事だと思います。

LMS【Least Mean Square】

{ 　h_m(n+1)=h_m(n)+𝜇e(n)x(n-m) }


NLMS【Normalized Least Mean Square】

{ 　\displaystyle h_m(n+1)=h_m(n)+𝜇\frac{e(n)x(n-m)}{\sum_{l=0}^{N-1} x^{2}(n-l)} }

　※𝜇：ステップサイズ

サンプルコード：NLMSで作ってみる

エコーキャンセラやノイズキャンセラを想定した音響信号処理向けで、16ビットのPCMを処理するサンプルコードです。シンプルに書きたかったので組み込み向けの最適化は考えずに、フレーム構成ではなくサンプル処理です。適応アルゴリズムにはNLMS、フィルタにはFIR【Finite Impulse Response】を使用しました。

FIRフィルタ

f:id:jspnet:20200306005824p:plain:right:w600 まずはフィルタ部のFIRです。

{ 　y(n) = \sum_{k=0}^{N}h_k x(n-k) }


サンプルコード

#define FIR_TAP          512     // 目的に合わせて調整する

// clip
static short double_saturat(double in)
{
    if( in >32767.0) in = 32767.0;
    if( in < -32768.0) in = - 32768.0;
    return (short)in;
}

// FIR
short fir_sample(float* coef, short *mem, short in )
{
    int i;
    double sum;
    short out;
    short *buff;

    for(i=0; i<(FIR_TAP-1); i++){
        mem[i] = mem[i+1];
    }
    mem[i]=in;

    sum=0.;
    buff = &mem[FIR_TAP - 1];

    // convolution
    for(i=0; i< FIR_TAP; i++){
        sum += coef[i] * (double)*buff--;
    }
    out = double_saturat(sum);
    return out;
}

ADF

適応アルゴリズムはNLMSです。目的信号d(n)からフィルタ出力y(n)を減算した誤差e(n)に基づいて、フィルタ係数h(n)を最適係数に近づけます。

{ 　e(n)=d(n)-y(n) }

誤差e(n)のパワーを最小にする。

{ 　\displaystyle h_m(n+1)=h_m(n)+𝜇\frac{e(n)x(n-m)}{\sum_{l=0}^{N-1} x^{2}(n-l)} }

　※𝜇：ステップサイズ

サンプルコード

#define ADF_STEP_SIZE    0.025   // 目的に合わせて調整する

// float coef[FIR_TAP];  // ゼロ初期化
// short mem[FIR_TAP];   // ゼロ初期化
// short x(n)
// short d(n)
// short e(n)

// Normalised least mean squares
void adf_sample(float *coef, short *mem, short x, short d, short *e)
{
    int i;
    double sum;
    short *buff;
    short y;

    // FIR
    y = fir_sample(coef, mem, x);
    *e = d - y; // error

    // NLMS
    sum = 1.0e-6;
    for (i = 0; i < FIR_TAP; i++) {
        sum += mem[i] * mem[i];
    }

    // update coef
    buff = &mem[FIR_TAP - 1];
    for (i = 0; i < FIR_TAP; i++) {
        coef[i] += (float)(ADF_STEP_SIZE * *e * *buff-- / sum);
    }
}