// adaptive filter : equalizer (reverse filter design)
//
//            +----------------+
//       +--->| modeling delay +---+ d
//       |    +----------------+   |
//       |                         v +
//  x ---+                  ^     (+)---+---> e
//       |                  |      ^ -  |
//       |     +---+      +---+    |    |
//       +---->| h +----->| w +----+    |
//             +---+  fx  +---+ y       |
//                          |           |
//                          +-----------+

clear;
stacksize('max');

DATLEN=100000;                 // simulation data length
FIRLEN=10;                     // unknown system firlter length
ADPLEN=100;                    // adaptive filter length
DLYLEN=50;                     // delay length
MU=0.0001;                     // step size parameter

rand('normal');                // Gaussian noise mode
indata=rand(1:DATLEN);         // generate noise

outdata=zeros(1:DATLEN);

//rand('seed', 1234);          // set random number seed
h=rand(1:FIRLEN);
hbuf=zeros(1:FIRLEN);
w=zeros(1:ADPLEN);
wbuf=zeros(1:ADPLEN);
delay=zeros(1:DLYLEN);

for i=1:DATLEN
    
  x=indata(i);
    
  d=delay(DLYLEN);
  delay=[x, delay(1:DLYLEN-1)];

  hbuf=[x, hbuf(1:FIRLEN-1)];
  fx=hbuf*h';
  wbuf=[fx, wbuf(1:ADPLEN-1)];
  y=wbuf*w';
  e=d-y;
  w=w+2*MU*e*wbuf;
  
  outdata(i)=e;
    
end

scf(0); clf; plot2d(h);
scf(1); clf; plot2d(w);
scf(2); clf; plot2d(outdata);
scf(3); clf; plot2d(10*log10(outdata^2+1e-9));

scf(4);
clf;
hspec=20*log10(abs(fft([h, zeros(1:2000-length(h))], -1))+1e-9);
wspec=20*log10(abs(fft([w, zeros(1:2000-length(w))], -1))+1e-9);
plot2d(hspec(1:1000), style=2);
plot2d(wspec(1:1000), style=5);

