function out = hoth(sf, len);

%   out = hoth(sf, len);
%       sf: sampling freq [Hz]
%      len: length of the Hoth signal [sample]
%
% generate Hoth spectral noise defined by ITU-T P.800
%
% above 8000 Hz, -15 dB/oct. slope is applied.
% below 100 Hz, IEC60120 octave band filter

f = [0 25 50 71 100	125	160	200	250	315	400	500 ...
     630	800	1000	1250	1600	2000	2500	3150 ...
     4000	5000	6300	8000  10000  12600 16000];
db = [-42.6 -10.1 14.4 32.4 32.4 30.9 29.1 27.6	26.0 24.4 22.7	21.1 ...
      19.5	17.8	16.2	14.6	12.9	11.3	9.6	7.8 ...
      5.4	2.6	-1.3	-6.6   -11.6    -16.6   -21.6];

blen = 2048;
% wgn = randmtn(len+blen, 1);
wgn = randn(len+blen, 1);

nlen = length(f);

if max(f) < sf/2
  f(nlen+1) = sf/2;
  db(nlen+1) = db(nlen) - 15*(log(f(nlen+1)/f(nlen))/log(2));
elseif max(f) > sf/2
  nf = find(f > sf/2);
  f = f(1:nf-1);
  db = db(1:nf-1);
end

b = fir2(blen, f/(sf/2), 10.^((db-db(1))/20));
out = conv(wgn, b);
out = out(blen/2+1:len+blen/2)/10000;

