// IIR Hilbert transform filter pair demo program (4th order)

clear;

INFILE='male_weather.wav';
OUTFILE='effect_'+INFILE;

[indata, label]=loadwave(INFILE);
fs=label(3);
datlen=length(indata);
outdata=zeros(1:datlen);

p1=-0.947997;
p2= 0.814444;
p3=-0.583754;
p4= 0.218006;

h1=p1*p3;
g1=p1+p3;
h2=p2*p4;
g2=p2+p4;

zx1i=0;  zx2i=0;
zy1i=0;  zy2i=0;
zzx1i=0; zzx2i=0;
zzy1i=0; zzy2i=0;

zx1q=0;  zx2q=0;
zy1q=0;  zy2q=0;
zzx1q=0; zzx2q=0;
zzy1q=0; zzy2q=0;

sumy=0;
for i=1:datlen
  x=indata(i);

  //----------------------------------------
  // IIR Hilbert transform filter (1/2)
  
  x1i=x;

  y1i=h1*x1i-g1*zx1i+zzx1i+g1*zy1i-h1*zzy1i;
  zzx1i=zx1i;
  zx1i=x1i;
  zzy1i=zy1i;
  zy1i=y1i;

  x2i=y1i;
  y2i=h2*x2i-g2*zx2i+zzx2i+g2*zy2i-h2*zzy2i;
  zzx2i=zx2i;
  zx2i=x2i;
  zzy2i=zy2i;
  zy2i=y2i;

  di=y2i;

  //----------------------------------------
  // IIR Hilbert transform filter (2/2)
  
  x1q=x;

  y1q=h1*x1q+g1*zx1q+zzx1q-g1*zy1q-h1*zzy1q;
  zzx1q=zx1q;
  zx1q=x1q;
  zzy1q=zy1q;
  zy1q=y1q;

  x2q=y1q;
  y2q=h2*x2q+g2*zx2q+zzx2q-g2*zy2q-h2*zzy2q;
  zzx2q=zx2q;
  zx2q=x2q;
  zzy2q=zy2q;
  zy2q=y2q;

  dq=y2q;

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

  y=sqrt(di^2+dq^2);       // calculate envelope of input signal
  sumy=0.9*sumy+0.1*y;     // calculate DC offset
  outdata(i)=y-sumy;       // remove DC offset
end

savewave(OUTFILE, 0.95*outdata/max(abs(outdata)), fs);