clear;
rand('uniform');
  
FS=8000;
LEN=2000;
ENVELOPE_SHAPE=0.2;
C_FREQ=261.62557;

toneheight=[0, 2, 4, 5, 7, 9, 11];

for scale_no=1:7
  x=zeros(1:LEN);
  f0=(C_FREQ/32)*2^(toneheight(scale_no)/12);
  for i=1:16
    if (100<f0) && (f0<(FS/2-100))
      x=x+sin(2.0*%pi*(1:LEN)*f0/FS+%pi*rand());
    end
    f0=2*f0;
  end
  
  // add harmonics
  f0=16.0*(C_FREQ/32)*2^(toneheight(scale_no)/12);  
  f1=3*f0;
  if f1<(FS/2-100) 
    x=x+sin(2.0*%pi*(1:LEN)*f1/FS+0.5*sin(2.0*%pi*(1:LEN)/733)+%pi*rand());
  end
  f1=5*f0;
  if f1<(FS/2-100) 
    x=x+sin(2.0*%pi*(1:LEN)*f1/FS+0.3*sin(2.0*%pi*(1:LEN)/661)+%pi*rand());
  end
  f1=6*f0;
  if f1<(FS/2-100) 
    x=x+sin(2.0*%pi*(1:LEN)*f1/FS+0.2*sin(2.0*%pi*(1:LEN)/503)+%pi*rand()); 
  end
  f1=10*f0;
  if f1<(FS/2-100) 
    x=x+sin(2.0*%pi*(1:LEN)*f1/FS+0.4*sin(2.0*%pi*(1:LEN)/467)+%pi*rand());     
  end
  f1=12*f0;
  if f1<(FS/2-100) 
    x=x+sin(2.0*%pi*(1:LEN)*f1/FS+0.7*sin(2.0*%pi*(1:LEN)/593)+%pi*rand());      
  end
  f1=14*f0;
  if f1<(FS/2-100) 
    x=x+sin(2.0*%pi*(1:LEN)*f1/FS+0.7*sin(2.0*%pi*(1:LEN)/859)+%pi*rand());
  end
  f1=18*f0;
  if f1<(FS/2-100) 
    x=x+sin(2.0*%pi*(1:LEN)*f1/FS+%pi*rand());      
  end
  f1=20*f0;
  if f1<(FS/2-100) 
    x=x+sin(2.0*%pi*(1:LEN)*f1/FS++%pi*rand());        
  end
  f1=26*f0;
  if f1<(FS/2-100) 
    x=x+sin(2.0*%pi*(1:LEN)*f1/FS+%pi*rand());        
  end                
   
  x=x.*(sin(%pi*(1:LEN)/LEN)^ENVELOPE_SHAPE);   
  savewave('tone'+msprintf('%i', toneheight(scale_no))+'.wav', 0.1*x, FS);
end