• 生体信号処理研究室
• 資料
• 信号処理
• フィルタ

LMS filter sample codes

最終更新

2011-07-05

著者

和田良太

LMS filter

Cによる実装

• lms.2011-06-24.zip

/*
 * main.c
 *
 *  - LMS algorithm sample code.
 *
 *   Author: Ryota Wada
 *     Date: 2011-06-24
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MU           0.01   // LMS アルゴリズムの mu に対応する値
#define M            31     // LMS アルゴリズムの M に対応する値
#define LENGTH       1000   // x, y, v, e の信号長（データ点数）
#define H_INIT_VALUE 0.1    // LMS アルゴリズムの h の全要素の初期値に対応する値

void lms (double *x, double *v, double *y, double *h, double *e);

/*
 * main: プログラムエントリポイント
 */

// 簡単のため，IO関連のエラーチェックはしていません。

int main (void) {
    int i;       // 汎用ループカウンタ
    FILE *fp;    // 汎用ファイルポインタ
    double *x;   // 行ベクトル ``x'' に対応する値格納変数（ユーザから与えられる）
    double *v;   // 行ベクトル ``v'' に対応する値格納変数（ユーザから与えられる）
    double *y;   // 行ベクトル ``y'' に対応する値格納変数（推定した出力信号）
    double *h;   // 列ベクトル ``h'' に対応する値格納変数（修正してゆくシステム）
    double *e;   // 行ベクトル ``e'' に対応する値格納変数（誤差）

    x  = (double *) malloc(LENGTH * sizeof(double));
    v  = (double *) malloc(LENGTH * sizeof(double)); 
    y  = (double *) malloc(LENGTH * sizeof(double));
    h  = (double *) malloc(     M * sizeof(double));
    e  = (double *) malloc(LENGTH * sizeof(double));

    fp = fopen("x.txt", "r");
    for (i = 0; i < LENGTH; i++) {
        fscanf(fp, "%lf", &x[i]);
    }
    fclose(fp);
    fp = fopen("v.txt", "r");
    for (i = 0; i < LENGTH; i++) {
        fscanf(fp, "%lf", &v[i]);
    }
    fclose(fp);
    for (i = 0; i < LENGTH; i++) {
        y[i] = 0.0;
        e[i] = 0.0;
    }
    for (i = 0; i < M; i++) {
        h[i] = H_INIT_VALUE;
    }

    // lms routine
    lms(x, v, y, h, e);
    
    // output ``y'' stream.
    fp = fopen("y.txt", "w");
    for (i = 0; i < LENGTH; i++) {
        fprintf(fp, "%.12lf\n", y[i]);
    }
    fclose(fp);
    // output ``e'' stream.
    fp = fopen("e.txt", "w");
    for (i = 0; i < LENGTH; i++) {
        fprintf(fp, "%.12lf\n", e[i]);//sqrt(e[i] * e[i]));
    }
    fclose(fp);

    free(x);
    free(v);
    free(y);
    free(h);
    free(e);

    return EXIT_SUCCESS;
}
/*
 * LMSアルゴリズムを用いた演算のルーチン
 */
void lms (double *x, double *v, double *y, double *h, double *e) {
    int i, k;

    for (i = M; i < LENGTH; i++) {
        for (k = 0; k < M; k++) {
            y[i] += h[k] * v[i - k];
        }
        e[i] = x[i] - y[i];
        for (k = 0; k < M; k++) {
            h[k] += MU * e[i] * v[i - k];
        }
    }
    return;
}

// EOF

Schemeによる実装

• lms.scm

;;;
;;; lms.scm
;;;   - LMS filter sample code
;;;
;;;         Author: Ryota Wada
;;;   Organization: BSP lab, Kinki University.
;;;           Date: 2011-07-05
;;; 
(define make-lms-result-object list)
(define lms-result-object-ref list-ref)
(define index-of-y 0)
(define index-of-e 1)
(define index-of-h 2)

(define lms
  (lambda (x v m mu)
    (let ([l (vector-length x)])
      (let update-lms ([y (make-vector l 0)]
                       [e (make-vector l 0)]
                       [h (make-vector m 0)]
                       [i m])
        (if (>= i l)
            (make-lms-result-object y e h i)
            (let* ([x-i (vector-ref x i)]
                   [y-i (let calc-sum ([k 0]
                                       [sum 0])
                          (if (>= k m)
                              sum
                              (calc-sum (+ k 1)
                                        (+ (* (vector-ref h k)
                                              (vector-ref v (- i k)))))))]
                   [e-i (- x-i y-i)])
              (update-lms (begin (vector-set! y i y-i) y)
                          (begin (vector-set! e i e-i) e)
                          (let update-h ([k 0])
                            (if (>= k m)
                                h
                                (begin
                                  (vector-set! h k (+ (vector-ref h k)
                                                      (* mu e-i (vector-ref v (- i k)))))
                                  (update-h (+ k 1)))))
                          (+ i 1))))))))
(define *length* 1000)
(define main
  (lambda args
    (print (lms-result-object-ref (lms (make-vector *length* 0.9)
                                       (make-vector *length* 1)
                                       31
                                       0.01)
                                  index-of-y))))
;; EOF