// http://www.cepstrum.co.jp/
// howling canceller simulation program for Scilab
//
//  +-------------------------------------------------------------------+
//  |                                                                   |
//  |                       Gain             s (quasi-white noise)      |
//  |                        |               |                          |
//  |                        v     +---+ d   v                          |
//  |                  +--->(x)--->| h +--->(+)---+ ds                  |
//  |                  |           +---+          |                     |
//  |    +---------+   |                        + v         +-------+   |
//  +--->| limiter +---+ x           ^           (+)---+--->| delay +---+
//       +---------+   |             |          - ^    | e  +-------+ ze
//                     |           +---+          |    |
//                     +---------->| w +----------+    |
//                                 +---+ y             |
//                                   |                 |
//                                   +-----------------+

clear;                         // clear all variables
  
//Gain=2.0;                      // maximum gain of acoustic system
//MU  =0.003;                    // step size parameter of adaptive filter

Gain=1000.0;                 // maximum gain of acoustic system
MU  =0.000000003;            // step size parameter of adaptive filter

DATLEN=10000;

FIRLEN=100;                    // adaptive filter legnth
DLYLEN=100;                    // delay length

// generate amplitude limited quasi-white noise
rand('uniform');
indata=zeros(1:DATLEN);
for i=1:DATLEN
  indata(i)=0.7*(sum(rand(1:20))-10.0)/10.0;   
end

outdata=zeros(1:DATLEN);
missalignment=zeros(1:DATLEN);
err=zeros(1:DATLEN);

// generate impulse response of acoustice system (h)
rand('normal');
rand('seed', 5678);
h=[zeros(1:5), rand(1:FIRLEN-5)];
h=h.*(0.96^(1:FIRLEN)); 
h=h/max(abs(fft([h, zeros(1:1000)], -1)));

buf=zeros(1:FIRLEN);
w=zeros(1:FIRLEN);
delay=zeros(1:DLYLEN);

ze=0;
for i=1:DATLEN 
  x=ze;

  // limiter
  if Gain<abs(x)
    x=Gain*sign(x);
  end

  buf=[x, buf(1:FIRLEN-1)];
  d=Gain*(h*buf');
  s=indata(i);
  ds=d+s;
  y=w*buf';
  e=ds-y;
  
  ze=delay(DLYLEN);
  delay=[e, delay(1:DLYLEN-1)];

  w=w+2*MU*e*buf;              // LMS algorithm

  outdata(i)=x;
  missalignment(i)=10.0*log10(sum(abs(w-Gain*h).^2)/sum((Gain*h).^2)+1e-10);
  err(i)=e-s;
end

// plot missalignment of adaptive filter coefficient [dB]
scf(1);
clf;
title('miss alignment [dB]', 'fontsize', 3);
plot2d(1:DATLEN, missalignment, style=2, rect=[0, -80, DATLEN, 0], axesflag=1);

// plot impulse response of h
scf(2);
clf;
title('imuplse response of h', 'fontsize', 3);
plot2d(1:FIRLEN, h, style=2, axesflag=1);       

// plot amplitude response of h
scf(3);
ampres=[h, zeros(1:1000)];
ampres=ampres(1:1000);
ampres=20*log10(abs(fft(ampres, -1)));
ampres=ampres(1:500);
clf;
title('amplitude response of h [dB]', 'fontsize', 3);
plot2d((1:500)/1000, ampres, style=2, axesflag=1, rect=[0, -30, 0.5, 0]);   

// file output
fd=mopen('out.txt', 'wt');
for i=1:DATLEN
  mfprintf(fd, '%5i %12.9f\n', i-1, missalignment(i));
end
mclose(fd);