Implementation in C of Least Mean Square (LMS) algorithm

[radustefan]

I'm a student. In a project for my Bachelor of Science Degree i have to implement in C a LMS algorithm. The algorithm is put in an IIR noth filter. The error signal for the adaptive filter is e(n)=-y(n). I implemented the algorithm but it doesen't work. I applied a signal at the input of the filter. The signal is a sinusoid limited at half of the amplitude. The filter should increase/decrease it's notch frequency until it matches the frequency of the signal. Instead the notch frequency becomes 0. Here is the C program, it is something wrong in it. I also put the matlab program with which I generate the signal for the C program. I also put a function the matlab program uses.

-----------------------------iir.c-------------------------------------

#include <prototype.h>

#include <stdio.h>

#define DataBlockSize 128    /* the size of a data block */

#define BlockLength 50        /* the number of the data blocks */

Word16 x[DataBlockSize],y[DataBlockSize],e[DataBlockSize];

Word16 w,w1,w2,_miu1=WORD16(0.05),_miu2=WORD16(0.0666),rho=WORD16(0.8),aux,aux2;

//Word16 b[]={WORD16(b0/2), WORD16(b1/2), WORD16(b2/2)};

//Word16 a[]={WORD16(a1/2), WORD16(a2/2)};

Word16 b[]={WORD16(0.5), WORD16(-1.4142/2), WORD16(0.5)};

Word16 a[]={WORD16(-1.1314/2), WORD16(0.64/2)};

Word32 sum;

int main()

{

short n,i;

FILE *fpx,*fpy;

fpx=fopen("x.dat","r+b");

if (!fpx)

    printf("\nIt didn't opened");

fpy=fopen("y.dat","w+b");

if (!fpy)

    printf("\nIt didn't opened");

for (i=0; i<BlockLength; i++)

{

    fread(x,sizeof(Word16),DataBlockSize,fpx);

    for (n=0; n<DataBlockSize; n++)

    {

        sum = L_deposit_h(shr(x[n],1));

        sum = L_msu(sum,a[0],w1);

        sum = L_msu(sum,a[1],w2);

        w = shl(round(sum),1);

        sum = L_mult(b[0],w);

        sum = L_mac(sum,b[1],w1);

        y[n]= shl(mac_r(sum,b[2],w2),1);

        w2=w1;

        w1=w;

        e[n]=-y[n];

        if(e[n]<0)

            aux=-div_s(-e[n],_miu1);

        else

            aux=div_s(e[n],_miu1);

        if(x[n]<0)

            aux2=-div_s(-x[n],_miu2);

        else

            aux2=div_s(x[n],_miu2);

        aux=mult(aux,aux2);

        b[1]=add(b[1],aux);

        //b[1]=b[1]+miu*e[n]*x[n];

        a[0]=mult(rho,b[1]);

        //a[0]=rho*b[1];

    }

    fwrite(y,sizeof(Word16),DataBlockSize,fpy);

}

fclose(fpx);

fclose(fpy);

return(0);

}

---------------------------matlab.m--------------------------------

clear all

close all

clc

Fs = 6400;

w0=2*pi/8;

miu=300;

b=[1 -2*cos(w0) 1]

a=[1 -2*0.8*cos(w0) 0.64]

figure(1),freqz(b,a)

figure(2),zplane(b,a)

% Simulation with a sinusoidal signal with the amplitude limited at half of the sinuoid's

% amplitude

t=0:1/Fs:1-1/Fs;              % 1 secs @ 6,4kHz sample rate

W=2*pi*Fs/12;

  for i=1:6400

      x(i)=min(cos(W*t(i)),0.5);

      x(i)=max(x(i),-0.5);

      %x(i)=cos(W*t(i));

  end

x=0.01*x;

k=fix(length(x)/128);

L=k*128;

fid=fopen('..\x.dat','w','b');

fwrite(fid,x.*2^15,'int16');

fclose(fid);

display('PRESS A KEY');

pause

fid=fopen('..\y.dat','r','b');

y=fread(fid,L,'int16');

fclose(fid);

y=y'/(2^15);

for i=1:L

   %adaptive filter

    if i<=2

        if i==1

            [y_ref(i),b,a,e]=filtru_adaptiv_notch([0 0 x(i)],[0 0],a,b,miu);

        else    

            [y_ref(i),b,a,e]=filtru_adaptiv_notch([0 x(i-1:i)],[0 y_ref(i-1)],a,b,miu);

        end

    else

        [y_ref(i),b,a,e]=filtru_adaptiv_notch(x(i-2:i),y_ref(i-2:i-1),a,b,miu);

    end

end

n=0:L-1;

w=-pi:2*pi/Fs:pi-2*pi/Fs;

figure(3),plot(w,abs(fftshift(fft(x))));%plot(n,x)

xlabel('Normalized frequencies');

ylabel('Amplitude');

w=-pi:2*pi/512:pi-2*pi/512;

figure(4),plot(w,abs(fftshift(fft(y_ref(6400-512+1:6400)))));

xlabel('Normalized frequencies');

ylabel('Amplitude');

figure(5),plot(w,abs(fftshift(fft(y(6400-512+1:6400)))));%plot(n,y),grid

xlabel('Normalized frequencies');

ylabel('Amplitude');

-----------------------adaptive_notch_filter.m---------------------------

function [y,b,a,e]=adaptive_notch_filter(x,y,a,b,miu)

r=0.8;

y=2*(x(3-0)*b(1+0)/2+(y(3-1)*(-a(1))/2+x(3-1)*b(1+1)/2+x(3-2)*b(1+2)/2+y(3-2)...

            *(-a(2))/2))

y=y/4;

e=-y;

b=[b(1) b(1+1)+miu*e*x(3) b(3)];

a=[a(1) r*b(2) a(3)];

[radustefan · Accepted solution]

I resolved the problem meanwhile, you don't have to find the problem.

[sauloqueiroz]

Hi radustefan,

Does your code solve the nonlinear least-square problem or have you any idea in this sense?

thank you in advance!