lms.c

dsp AD公司ADSP21的代码,里面有FFT FIR IIR EQULIZER G722_21F 等可以在项目中直接应用的代码.此代码的来源是ADI公司自己出版的书籍,此书在美国购得

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "rtdspc.h"

#define N 351
#define L 20            /* filter order, (length L+1) */

/* set convergence parameter */
    float mu = 0.01;

void main()
{
    float lms(float,float,float *,int,float,float);
    static float d[N],b[21];
    float signal_amp,noise_amp,arg,x,y;
    int k;

/* create signal plus noise */
    signal_amp = sqrt(2.0);
    noise_amp = 0.2*sqrt(12.0);
    arg = 2.0*PI/20.0;
    for(k = 0 ; k < N ; k++) {
        d[k] = signal_amp*sin(arg*k) + noise_amp*gaussian();
    }

/* scale based on L */
    mu = 2.0*mu/(L+1);

    x = 0.0;
    for(k = 0 ; k < N ; k++) {
        sendout(lms(x,d[k],b,L,mu,0.01));
/* delay x one sample */
        x = d[k];
    }
}

/*
      function lms(x,d,b,l,mu,alpha)

Implements NLMS Algorithm b(k+1)=b(k)+2*mu*e*x(k)/((l+1)*sig)

x      = input data
d      = desired signal
b[0:l] = Adaptive coefficients of lth order fir filter
l      = order of filter (> 1)
mu     = Convergence parameter (0.0 to 1.0)
alpha  = Forgetting factor   sig(k)=alpha*(x(k)**2)+(1-alpha)*sig(k-1)
         (>= 0.0 and < 1.0)

returns the filter output
*/

float lms(float x,float d,float *b,int l,
                  float mu,float alpha)
{
    int ll;
    float e,mu_e,lms_const,y;
    static float px[51];      /* max L = 50 */
    static float sigma = 2.0; /* start at 2 and update internally */

    px[0]=x;

/* calculate filter output */
    y=b[0]*px[0];
    for(ll = 1 ; ll <= l ; ll++)
        y=y+b[ll]*px[ll];

/* error signal */
    e=d-y;

/* update sigma */
    sigma=alpha*(px[0]*px[0])+(1-alpha)*sigma;
    mu_e=mu*e/sigma;

/* update coefficients */
    for(ll = 0 ; ll <= l ; ll++)
        b[ll]=b[ll]+mu_e*px[ll];
/* update history */
    for(ll = l ; ll >= 1 ; ll--)
        px[ll]=px[ll-1];

    return(y);
}