// howling simulation program for Scilab
//
//  +----------------------------------------------+
//  |                                              |
//  |        x  +-----+ r  +------+ y  +-------+ m |
//  +--->(+)--->| FIR +--->| GAIN +--->| LIMIT +---+
//        ^     +-----+    +------+    +-------+    
//        |
//        |
//        S

clear;
rand('uniform');
rand('seed', 0);

SIMLEN=10000;
FIRLEN=20;      
FC    =7;              
DLYLEN=20;           
GAIN  =1.1;
LIMIT =1.0;
TRIGLVL=0.0001;

out=zeros(1:SIMLEN);
h=zeros(1:FIRLEN);
firbuf=zeros(1:FIRLEN+DLYLEN-1);

t=[0, 0.7, 0.7, 0.7, 0.85, 1.0, 0.85, 0.6, 0.7, 0.37, 0.0];
t=[t, mtlb_fliplr(t(2:10))];
h=real(fft(t, 1));
h=[zeros(1:DLYLEN), h(FIRLEN/2+2:FIRLEN), h(1:FIRLEN/2)];
hspec=fft([h, zeros(1:2000-length(h))], -1);
hspec=20*log10(abs(hspec(1:1000))+1e-100);

m=0.0;
for i=1:SIMLEN
  // generate trigger signal
  if i<=FIRLEN
    s=(2.0*rand(1)-1.0)*TRIGLVL;
  else
    s=0.0;
  end
  // add trigger signal
  x=m+s;
  // FIR filter
  firbuf=[x, firbuf(1:FIRLEN+DLYLEN-2)];
  r=h*firbuf';
  // apply gain
  y=r*GAIN;
  // limiter
  if LIMIT<=abs(y)
    m=LIMIT*sign(y);
  else
    m=y;
  end
  // save output data
  out(i)=y;
end

scf(1); clf; plot(out);     // plot howling signal waveform
scf(2); clf; plot(h);       // plot impulse response
scf(3); clf; plot(hspec);   // plot amplitude response

savewave('howling001out.wav', 0.95*out/max(abs(out)));