// cepstrum.co.jp
//
// decorrelator effect simulation program of adaptive howling canceller
//
//-------------------------------------------------------------------------
// decorrelator off
//
//  +---------------------------------------------+
//  |                                             | 
//  |                                   s         |
//  |                                   |         |
//  |                                   |         |
//  |                       +-----+ d   v         |
//  |                  +--->| FIR +--->(+)---+    |
//  |                  |    +-----+          |    |
//  |                  |                   + |    |
//  |    +---------+   |                     v  e |
//  +--->| limiter +---+ x                  (+)---+
//       +---------+   |       ^             ^    |
//                     |       |           - |    |
//                     |    +-----+ y        |    |
//                     +--->| ADP +----------+    |
//                          +-----+               |
//                             |                  |
//                             |                  |
//                             +------------------+
//
//-------------------------------------------------------------------------
// decorrelator on (delay on)
//
//  +--------------------------------------------------------------+
//  |                                                              | 
//  |                                   s                          |
//  |                                   |                          |
//  |                                   |                          |
//  |                       +-----+ d   v                          |
//  |                  +--->| FIR +--->(+)---+                     |
//  |                  |    +-----+          |                     |
//  |                  |                   + |                     |
//  |    +---------+   |                     v  e      +-------+   |
//  +--->| limiter +---+ x                  (+)---+--->| delay +---+ ze
//       +---------+   |       ^             ^    |    +-------+
//                     |       |           - |    |
//                     |    +-----+ y        |    |
//                     +--->| ADP +----------+    |
//                          +-----+               |
//                             |                  |
//                             |                  |
//                             +------------------+
//
//-------------------------------------------------------------------------

clear;

rand('normal');
rand('seed', 234);

DATLEN=1000000;   // should be odd number
ADPLEN=100;        
DLYLEN=10;
MU    =0.001;
THLD  =2.0;

indata=rand(1:DATLEN);
indata=convol([2.0, 0.5, 0.5, 0.0, 1.0, 1.0, 1.0, 0.0, 0.5, 0.5 ], indata);
indata=indata(1:DATLEN);
indata=0.9*indata/max(abs(indata));

savewave('in.wav', indata, 8000);

outdata=zeros(1:DATLEN);

edly=zeros(1:DLYLEN);

h=rand(1:ADPLEN).*(0.97^(1:ADPLEN));
h=h-sum(h)/ADPLEN;
scf(0);
clf;
plot(h);

xbuf=zeros(1:ADPLEN);
w_decor_on=zeros(1:ADPLEN);
w_decor_off=zeros(1:ADPLEN);

out_decor_on=zeros(1:DATLEN);
out_decor_off=zeros(1:DATLEN);

//--------------------------------------------------------------
// decorrelator off

mprintf('running decorrelator off simulation\n');

e=0;
for i=1:DATLEN

  if THLD<abs(e) 
    x=sign(e)*THLD;       // hard limiter
  else
    x=e;
  end

  xbuf=[x, xbuf(1:ADPLEN-1)];
  d=h*xbuf';
  y=w_decor_off*xbuf';
  s=indata(i);
  e=(d+s)-y;
  w_decor_off=w_decor_off+2.0*MU*e*xbuf;

  out_decor_off(i)=x;

end
savewave('out_decor_off.wav', 0.9*out_decor_off/THLD, 8000);

//--------------------------------------------------------------
// decorrelator on

mprintf('running decorrelator on simulation\n');

ze=0;
for i=1:DATLEN

  if THLD<abs(ze) 
    x=sign(ze)*THLD;       // hard limiter
  else
    x=ze;
  end

  xbuf=[x, xbuf(1:ADPLEN-1)];
  d=h*xbuf';
  y=w_decor_on*xbuf';
  s=indata(i);
  e=(d+s)-y;
  w_decor_on=w_decor_on+2.0*MU*e*xbuf;
  ze=edly(DLYLEN);
  edly=[e, edly(1:DLYLEN-1)];

  out_decor_on(i)=x;

end
savewave('out_decor_on.wav', 0.9*out_decor_on/THLD, 8000);

//--------------------------------------------------------------

spec_in=indata(DATLEN-8192:DATLEN-1);
spec_in=20.0*log10(fft(spec_in, -1));
spec_in=spec_in(1:4096);
scf(1);
clf;
plot2d(spec_in, rect=[0, -40, 4096, 40], axesflag=1, style=1);  // black

spec_decor_off=out_decor_off(DATLEN-8192:DATLEN-1);
spec_decor_off=20.0*log10(fft(spec_decor_off, -1));
spec_decor_off=spec_decor_off(1:4096);
scf(2);
clf;
plot2d(spec_decor_off, rect=[0, -40, 4096, 40], axesflag=1, style=2);  // blue

spec_decor_on=out_decor_on(DATLEN-8192:DATLEN-1);
spec_decor_on=20.0*log10(fft(spec_decor_on, -1));
spec_decor_on=spec_decor_on(1:4096);
scf(3);
clf;
plot2d(spec_decor_on, rect=[0, -40, 4096, 40], axesflag=1, style=5);  // read