%
%
%   shepard_tone.m  version 1.0    December 7, 2012
%   By Tom Irvine   Email:  tomirvine@aol.com
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear a;
clear TT;
clear sine_sweep;
clear arg;
clear noc;
clear theta;
%
%
number_tones=4;           % number of tones
f1 = 27.5;                % starting freq (Hz)
f2 = (2^number_tones)*f1; % ending freq (Hz)
%
tmax=60;    % duration (sec)
sr=22050;   % sample rate (samples/sec)
%
%   log sweep rate
%
dt=1./sr;    % time step
ns=tmax*sr;  % number of sample/cycles
cycles=1;    % number of cycles
%  
oct=log(f2/f1)/log(2.);   % number of octaves
%
out1 = sprintf('\n Number of octaves = %8.2f  ',oct);    
disp(out1);       
%
ntimes = ns*cycles;
%
t2=tmax;
dur=tmax;
%
tpi=2.*pi;
%
rate=oct/dur;
%
out1 = sprintf('\n Sweep rate = %8.2f oct/sec ',rate);    
disp(out1);  
out1 = sprintf('\n            = %8.2f  oct/min \n',60.*rate);    
disp(out1);     
%
disp(' Calculating time vector... ')
TT=linspace(dt,(ntimes+1)*dt,ntimes); 
%
a = zeros(1,ntimes);
arg = zeros(1,ntimes);
freq = zeros(1,ntimes);
%
disp(' ')
disp(' Calculating log sweep....')
disp(' ')    
%  
%			fspectral = f1*pow(2.,rate*t);
%
a=zeros(1,ntimes);
theta=zeros(number_tones,1);
arg=zeros(number_tones,1);
key=zeros(number_tones,1);
noc=zeros(number_tones,1);
amp=zeros(number_tones,1);
%
for i=1:number_tones
    noc(i)=2^(i-1);
end
%
fspan=f2-f1;
theta_coefficient=tpi/(rate*log(2));
%
sd=150;
progressbar % Create figure and set starting time
for(i=1:ntimes)
%
    progressbar(i/ntimes); % Update figure  
%
    W=-1.+2^(rate*TT(i));
%
    arg(1)=f1*W;
%
    for j=1:number_tones
        if(key(j)==0)
            arg(j)=arg(1)*noc(j);       
        else
            W=-1.+2^(rate*TT(i-key(j)));
            arg(j)=f1*W;    
         end
    end
%
    for j=1:number_tones
        fspectral=noc(j)*f1*(2^(rate*TT(i-key(j))));
        delta=(fspectral-f1)/fspan;
        amp(j)=(1-cos(tpi*delta));
%
        if(fspectral>f2)
             key(j)=i;
             noc(j)=1;
        end
    end
%
    theta=theta_coefficient*arg;
%    
    a(i)=sum(amp.*sin(theta));
%
end
%
a(1)=0.;
TT=TT-dt;
%
% disp(' ')
% disp(' End of sweep ')
% disp(' ')
%
clear aa;
aa=detrend(a);
sine_sweep=[TT;aa]'; 
out5 = sprintf('\n\n The time history is renamed as "sine_sweep" ');
disp(out5); 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
np=ntimes;
%
disp(' ')
disp(' End of sweep ')
disp(' ')
clear signal;
signal=aa';
signal = signal/max(abs(signal));
disp(' ');
disp(' Enter output wav filename ');
out_file=input(' ','s');
wavwrite(signal,sr,out_file)
%