// adaptive filter : predictor (beat noise canceller)
//
//       +-------------------------+ d
//       |                         |
//       |                         v +
//  x ---+                   ^    (+)---+---> e
//       |                   |     ^ -  |
//       |    +-------+    +---+   |    |
//       +--->| delay +--->| w +---+    |
//            +-------+ zx +---+ y      |
//                           |          |
//                           +----------+

clear;
stacksize('max');

DATLEN=100000;                                       // simulation data length
FLTLEN=100;                                          // filter length
DLYLEN=5;                                            // delay length
MU=0.00003;                                          // step size parameter

indata=loadwave('male.wav');                         // male voice
//indata=loadwave('female.wav');                     // female voice
DATLEN=min([DATLEN, length(indata)]);
indata=indata(1:DATLEN);
for i=1:DATLEN
  if i<DATLEN/2
    indata(i)=indata(i)+0.7*sin(2*%pi/12.3*i);       // add sinusoidal signal
    indata(i)=indata(i)+0.3*sin(2*%pi/73.3*i);       // add sinusoidal signal    
  else
    indata(i)=indata(i)+0.4*sin(2*%pi/22.2*i);       // add sinusoidal signal
    indata(i)=indata(i)+0.6*sin(2*%pi/5.6*i);        // add sinusoidal signal
  end
end

outdata=zeros(1:DATLEN);
out_e=zeros(1:DATLEN);
out_y=zeros(1:DATLEN);

delay=zeros(1:DLYLEN);
xbuf=zeros(1:FLTLEN);
w=zeros(1:FLTLEN);

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

  xbuf=[zx, xbuf(1:FLTLEN-1)];
  y=xbuf*w';
  e=d-y;
  w=w+2*MU*e*xbuf;
  
  out_e(i)=e;
  out_y(i)=y;
    
end

scf(0); clf; plot2d(w);
scf(1); clf; plot2d(indata);
scf(2); clf; plot2d(out_e);
scf(3); clf; plot2d(out_y);

savewave('in_beat_cancel.wav', 0.99*indata/max(abs(indata)), 8000);
savewave('out_beat_cancel.wav', 0.99*out_e/max(abs(out_e)), 8000);

