% function alpha_M = eflslbr_tv1(lambda,M,delta) 

% 

% Implementation of the Recursive Least Squares Lattice (Recursive LSL)  

% adaption algorithm using a priori estimation errors and error 

% feedback. (Carezia version using a modified EF-LSL algorithm) 

% 

% use int2log and log2int routines 

% 

% lambda- Forgetting factor  

% delta- Small positive initialization constant 

% 

% alpha_M- A priori estimation error 

% 

% Felix Albu ---- Dublin ---- 10 October 2001 

% F. Albu, J. Kadlec, N. Coleman, A. Fagan, "Pipelined Implementations of the A Priori Error-Feedback LSL Algorithm Using Logarithmic Number System", 

% Proceedings of ICASSP 2002, pp.2681-2684, Orlando, Florida, U.S.A

% felix_albu@ieee.org

%

% ------------------------------------------------------------------------------ 

% Set default values 

  

clear 

load ca2 

  

%load datagen 

%d=yy4; 

%u=yy6; 

%L = length(d); 

  

prec=23; 

  

L= 10; 

M = 5; 

  

d=mtlb_c(1:L); 

u=mtlb_d(1:L);

 

ufpga= (2^(23-prec))*fix((2^prec)* u);

u=fix((2^prec)* u)/((2^prec));

for k=1:L

lu(k)=int2log2(ufpga(k));

end

 

dfpga= (2^(23-prec))*fix((2^prec)* d);

d=fix((2^prec)* d)/((2^prec));

for k=1:L

ld(k)=int2log2(dfpga(k));

end

 

 

fd=d2sp(d);fu=d2sp(u);

 

lambda = 0.999;

flambda=d2sp(lambda);

llambda=d2logs(lambda);

 

teta=2^(-7);

fteta=d2sp(teta);

lteta=d2logs(teta);

 

delta = lambda*teta*(1-lambda);

fdelta=d2sp(delta);

ldelta=d2logs(delta);

 

sirm=d2logs(zeros(M,1));

 

% ------------------------------------------------------------------------------

% The initialisation. (variables at time n=0)

 

gamma_old = ones(M,1);

psi_old = zeros(M,1);

F_old = ones(M,1) .* delta;

B_old = ones(M,1) .* delta;

Gamma_f_old = zeros(M,1);

Gamma_b_old = zeros(M,1);

kappa_old = zeros(M,1);

bn_old = zeros(M,1);

bn = zeros(M,1);

 

fgamma_old = ones(M,1);

fpsi_old = zeros(M,1);

fF_old = d2sp(F_old); % fF for m=M-1 is not calculated

fB_old = d2sp(B_old);

fGamma_f_old = zeros(M,1);

fGamma_b_old = zeros(M,1);

fkappa_old = zeros(M,1);

fbn_old = zeros(M,1);

fbn = zeros(M,1);

 

lgamma_old = d2logs(gamma_old);

lpsi_old = sirm;

lF_old = d2logs(F_old);

lB_old = d2logs(B_old);

lGamma_f_old = sirm;

lGamma_b_old = sirm;

lkappa_old = sirm;

lbn_old = sirm;

lbn = sirm;

 

% ------------------------------------------------------------------------------

% Declare variables

 

gamma= ones(M,1);

psi= zeros(M,1);

F= F_old;

B= B_old;

Gamma_b= zeros(M,1);

Gamma_b= zeros(M,1);

kappa= zeros(M,1);

eta= zeros(M,1);

 

fgamma= ones(M,1);

fpsi= zeros(M,1);

fF= fF_old;

fB= fB_old;

fGamma_b= zeros(M,1);

fGamma_b= zeros(M,1);

fkappa= zeros(M,1);

feta= zeros(M,1);

 

lgamma= d2logs(gamma);

lpsi= sirm;

lF= lF_old;

lB= lB_old;

lGamma_b= sirm;

lGamma_b= sirm;

lkappa= sirm;

leta= sirm;

 

alpha_M = zeros(1,L);

falpha_M = zeros(1,L);

lalpha_M = d2logs(zeros(1,L));

 

% ------------------------------------------------------------------------------

% After the initialisation (for n>=1)

%

 

for n = 1:L

 

% Order initialization for each n>=1

 

gamma(1)= 1;fgamma(1) = 1;lgamma(1)= d2log(1);

 

eta (1)= u(n);psi(1) = u(n);

feta(1)= fu(n);fpsi(1) = fu(n);

leta (1)= lu(n);lpsi(1) = lu(n);

alpha= d(n);falpha = fd(n);lalpha= ld(n);

 

% Order update of the lattice filter

for mm=2:M

 

% Line (1)

eta(mm)= eta(mm-1) - Gamma_f_old(mm) * psi_old(mm-1);

feta (mm)= fsub(feta(mm-1) , fmul(fGamma_f_old(mm) , fpsi_old(mm-1)));

leta(mm)= lsub(leta(mm-1) , lmul(lGamma_f_old(mm) , lpsi_old(mm-1)));

 

 

% Line (2)

f = gamma_old(mm-1) * eta(mm-1);

ff = fmul(fgamma_old(mm-1), feta(mm-1));

lf = lmul(lgamma_old(mm-1), leta(mm-1));

 

 

% Line (3)

psi(mm)= psi_old(mm-1) - Gamma_b_old(mm) * eta(mm-1);

fpsi(mm)= fsub(fpsi_old(mm-1) , fmul(fGamma_b_old(mm) , feta(mm-1)));

lpsi(mm)= lsub(lpsi_old(mm-1) , lmul(lGamma_b_old(mm) , leta(mm-1)));

 

 

% Line (4)

b = gamma(mm-1) * psi(mm-1);

fb = fmul(fgamma(mm-1), fpsi(mm-1));

lb = lmul(lgamma(mm-1), lpsi(mm-1));

 

 

% Line (5)

alpha = alpha - kappa_old(mm) * psi(mm-1);

falpha = fsub(falpha , fmul(fkappa_old(mm) , fpsi(mm-1)));

lalpha = lsub(lalpha , lmul(lkappa_old(mm) , lpsi(mm-1)));

 

 

% Line (6)

Gamma_f(mm) = Gamma_f_old(mm) + bn_old(mm-1)*eta(mm);

fGamma_f(mm) = fadd(fGamma_f_old(mm) , fmul(fbn_old(mm-1), feta(mm)));

lGamma_f(mm) = ladd(lGamma_f_old(mm) , lmul(lbn_old(mm-1), leta(mm)));

 

 

% Line (7)

F(mm-1)= lambda * F_old(mm-1) + f*eta(mm-1) + teta;

fF(mm-1)= fadd(fadd(fmul(flambda , fF_old(mm-1)) , fmul(ff, feta(mm-1))), fteta);

lF(mm-1)= ladd(ladd(lmul(llambda , lF_old(mm-1)) , lmul(lf, leta(mm-1))), lteta);

 

% Line (8)

B(mm-1)= lambda * B_old(mm-1) + b*psi(mm-1) + teta;

fB(mm-1)= fadd(fadd(fmul(flambda , fB_old(mm-1)) , fmul(fb, fpsi(mm-1))), fteta);

lB(mm-1)= ladd(ladd(lmul(llambda , lB_old(mm-1)) , lmul(lb, lpsi(mm-1))), lteta);

 

 % Line (9)

fn = f/F(mm-1);

ffn = fdiv(ff, fF(mm-1));

lfn = ldiv(lf, lF(mm-1));

 

 

% Line (10)

bn(mm-1) = b/B(mm-1);

fbn(mm-1) = fdiv(fb, fB(mm-1));

lbn(mm-1) = ldiv(lb, lB(mm-1));

 

 

% Line (11)

Gamma_b(mm) = Gamma_b_old(mm) + fn*psi(mm);

fGamma_b(mm) = fadd(fGamma_b_old(mm) , fmul(ffn, fpsi(mm)));

lGamma_b(mm) = ladd(lGamma_b_old(mm) , lmul(lfn, lpsi(mm)));

 

 

% Line (12)

kappa(mm)= kappa_old(mm) + bn(mm-1)*alpha;

fkappa(mm)= fadd(fkappa_old(mm) , fmul(fbn(mm-1), falpha));

lkappa(mm)= ladd(lkappa_old(mm) , lmul(lbn(mm-1), lalpha));

 

 

% Line (13)

gamma(mm) = gamma(mm-1) - bn(mm-1)*b;

fgamma(mm)= fsub(fgamma(mm-1), fmul(fbn(mm-1), fb));

lgamma(mm)= lsub(lgamma(mm-1), lmul(lbn(mm-1), lb));

 

end

 

 

% Update output

 

alpha_M(n) = alpha;

falpha_M(n) = falpha;

lalpha_M(n)=log2int2(lalpha)/2^prec;

 

erlns(n)=alpha_M(n)-lalpha_M(n);

ersp(n)=alpha_M(n)-falpha_M(n);

 

% Update old values

gamma_old= gamma;

psi_old= psi;

F_old= F;

B_old= B;

Gamma_f_old = Gamma_f;

Gamma_b_old = Gamma_b;

kappa_old= kappa;

bn_old= bn;

 

fgamma_old= fgamma;

fpsi_old= fpsi;

fF_old= fF;

fB_old= fB;

fGamma_f_old = fGamma_f;

fGamma_b_old = fGamma_b;

fkappa_old= fkappa;

fbn_old= fbn;

 

lgamma_old= lgamma;

lpsi_old= lpsi;

lF_old= lF;

lB_old= lB;

lGamma_f_old = lGamma_f;

lGamma_b_old = lGamma_b;

lkappa_old= lkappa;

lbn_old= lbn;

 

end

 

ssl=0;ssf=0;

sl=zeros(L,1);sf=zeros(L,1);

 

for i=1:L,

ssl=ssl+abs(erlns(i)); sl(i)=ssl;

ssf=ssf+abs(ersp(i)); sf(i)=ssf;

end;

 

 

fid=fopen('.\Debug\input.dat','w'); % The Handel C input file

fprintf(fid, '%d\n', L);

fprintf(fid, '%d\n', M);

fprintf(fid, '%d\n', ldelta);

fprintf(fid, '%d\n', llambda);

fprintf(fid, '%d\n', lteta);

for i=1:L

fprintf(fid,'%d\n', ufpga(i));

fprintf(fid,'%d\n', dfpga(i));

end

fclose(fid);

 

%

% ======================== end of file writing ====================

HC_JOB='handelc -s -ss 1000 eflslbr_tv1.c';

!del outfpga.dat

dos('del input.dat');

dos('copy .\Debug\input.dat input.dat');

 

dos(HC_JOB);

 

load outfpga.dat

 

r=outfpga-lalpha_M'*2^23;

if sum(r)==0

disp(['test is OK']);

end;

 

 

 

 

/*------------------------------------------------------

eflslbr01.c: 10.10.2001(on chip RAM) 

LOG LATTICE RECURSIVE LEAST-SQUARES 

 A PRIORI ERROR FEEDBACK ALGORITHM (CAREZIA's version) 

it takes the initial varibles values from the file input.dat

returns the output errors (LOG) in outfpga.dat

 Contact information: 

 felix_albu@ieee.org 

*/

#define PP1000_32BIT_RAMS

#define PP1000_DIVIDE1

#define PP1000_CLOCK PP1000_MCLK

#include "Pp1000.h"

#include <stdlib.h>

#include "stdlib.c"

set part = "V2000BG560-6";

#include "lalu.h"

int 32 idk, iuk, ialpha, luk, ldk, ldelta, llambda, lf, lb, lfn, leta, lalpha;

int 32ltemp3, ltemp1, ltemp2, ltemp4, ltemp5, ltemp6, ltemp7, ltemp8, ltemp9, ltemp10, ltemp11;

ram int 32 lgamma_old[16], lpsi_old[16], lF_old[16], lB_old[16], lGamma_f_old[16], lGamma_b_old[16], lkappa_old[16];

ram int 32 lgamma[16], lbn[16], lpsi[16], lF[16], lB[16], lGamma_b[16], lGamma_f[16], lkappa[16], lbn_old[16];

unsigned 4 lM1, lj;

unsigned 16 lindex;

unsigned 16 lL1;

int 32 lM, lL;

macro proc eflatice()

{

 par

 {

 ltemp1=lm(lGamma_f_old[lj], lpsi_old[lj-1]); // ltemp1=Gamma_f_old(mm) * psi_old(mm-1);

 lf = lm(lgamma_old[lj-1], leta[lj-1]); // (2) !!!!! f = gamma_old(mm-1) * eta(mm-1);

 ltemp2= lm(lGamma_b_old[lj],leta[lj-1]);// ltemp2=Gamma_b_old(mm) * eta(mm-1);

 lb = lm(lgamma[lj-1], lpsi[lj-1]); // (4) !!!!! lb=gamma(mm-1) * psi(mm-1)

 ltemp3= lm(lkappa_old[lj],lpsi[lj-1]);// ltemp3=kappa_old(mm) *psi(mm-1);

 }

 par

 {

 lsub(leta[lj-1], ltemp1, leta[lj]); // (1)!!! eta(mm)= eta(mm-1) + Gamma_f_old(mm) * psi_old(mm-1);

 {

 ltemp7= lm(lb, lpsi[lj-1]); // ltemp7= b*psi(mm-1)

 }

 }

 par

 {

 lsub(lpsi_old[lj-1], ltemp2, lpsi[lj]);// (3) !!!!!!! psi(mm) = psi_old(mm-1) - Gamma_b_old(mm) * eta(mm-1);

 {

 ltemp6= lm(lf, leta[lj-1]); // ltemp6= f*eta(mm-1)

 ltemp8= lm(lbn_old[lj-1], leta[lj]);// ltemp8= bn_old(mm-1)*eta(mm);

 ltemp4= lm(llambda, lF_old[lj-1]); // ltemp4=lambda * F_old(mm-1); 

 ltemp5= lm(llambda, lB_old[lj-1]); // ltemp5=lambda * B_old(mm-1);

 }

 }

 lsub(lalpha, ltemp3, lalpha); //(5) !!!! alpha = alpha - kappa_old(mm) *psi(mm-1);

 ladd(lGamma_f_old[lj], ltemp8, lGamma_f[lj]); // (6) !!! Gamma_f(mm) = Gamma_f_old(mm) + bn_old(mm-1)*eta(mm);

 ladd(ltemp4, ltemp6, ltemp10); 

  ladd(ltemp10, lteta, lF[lj-1]); // (7) !!F(mm-1)= lambda * F_old(mm-1) + f*eta(mm-1) + teta;

 ladd(ltemp5, ltemp7, ltemp11); 

 ladd(ltemp11, lteta, lB[lj-1]); // (8) !!B(mm-1)= lambda * B_old(mm-1) + b*psi(mm-1) + teta;

 par

 {

 lfn= ld(lf, lF[lj-1]); // (9) !!!! fn = f/F(mm-1);

 lbn[lj-1]= ld(lb, lB[lj-1]); // bn(mm-1) = b/B(mm-1);

 }

 ltemp9= lm(lbn[lj-1], lb); // ltemp9=bn(mm-1)*b;

 par

 {

 lsub(lgamma[lj-1], ltemp9, lgamma[lj]); // (13) !! gamma(mm) = gamma(mm-1) - bn(mm-1)*b;

 {

 ltemp12= lm(lfn, lpsi[lj]); // ltemp4= fn(mm-1)*psi(mm);

 ltemp13= lm(lbn[lj-1], lalpha); // ltemp5= bn(mm-1)*alpha;

 }

 }

 ladd(lGamma_b_old[lj], ltemp12, lGamma_b[lj]);// (11) !! Gamma_b(mm) = Gamma_b_old(mm) + fn*psi(mm);

 ladd(lkappa_old[lj], ltemp13, lkappa[lj]); // (12) !!kappa(mm)= kappa_old(mm) + bn(mm-1)*alpha;

}

void main(void)

{  

 lL=pc2alg();

 lM=pc2alg();

 ldelta=pc2alg();

 llambda=pc2alg();

 lteta=pc2alg();

 par

 {

 lM1=(unsigned 4)lM[3:0];

 lL1=(unsigned 16)lL[15:0];

 }

 for (lj=0; lj<lM1; lj=lj+1)

 {

 par

 {

 lgamma_old[lj]=(int 32)0;

 lpsi_old[lj]=ZERO_LNSALU;

 lB_old[lj]=ldelta;

 lF_old[lj]=ldelta;

 lGamma_b_old[lj]=ZERO_LNSALU;

 lGamma_f_old[lj]=ZERO_LNSALU;

 lkappa_old[lj]=ZERO_LNSALU;

 lgamma[lj]=(int 32)0;

 lbn[lj]=ZERO_LNSALU;

 lpsi[lj]=ZERO_LNSALU;

 lB[lj]=ldelta;

 lF[lj]=ldelta;

 lGamma_b[lj]=ZERO_LNSALU;

 lGamma_f[lj]=ZERO_LNSALU;

 lkappa[lj]=ZERO_LNSALU;

 lbn_old[lj]=ZERO_LNSALU;

 }

 }

 for (lindex=0; lindex<lL1; lindex=lindex+1)

 {

iuk=pc2alg();

 luk=int2log(iuk);

idk=pc2alg();

 ldk=int2log(idk);

 par

 {

 lgamma[0]=(int 32)0;

 leta=luk;

 lpsi[0]=luk; 

 lalpha=ldk;

 }

 for (lj=1; lj<lM1; lj=lj+1)// ciclu lj 

 eflatice();

 par

 {

 ialpha=log2int(lalpha);

 {

 for (lj=0; lj<lM1; lj=lj+1)

 {

 par

 {

 lgamma_old[lj]=lgamma[lj];

 lpsi_old[lj]=lpsi[lj];

 lF_old[lj]=lF[lj];

 lB_old[lj]=lB[lj];

 lGamma_f_old[lj]=lGamma_f[lj];

 lGamma_b_old[lj]=lGamma_b[lj];

 lkappa_old[lj]=lkappa[lj];

 lbn_old[lj]=lbn[lj];

 }

 }

 }

 }

 alg2pc(ialpha);

 } // end program EFLSL (simplified version)