%Sampling and playback rate of 22050 samples per second
Fs=22050;
t=0:1/Fs:1;
%w is the angular frequency of the fundamental (100 Hz)
w=2*pi*100;
y1=sin(w*t);
y2=sin(2*w*t);
y4=sin(4*w*t);
y8=sin(8*w*t);
y16=sin(16*w*t);
y32=sin(32*w*t);
y64=sin(64*w*t);
%the two tones ynote1 and ynote2 are so-called Shepard tone consisting
%of harmonics 1,2,4,8,16,32,64.  This harmonic sequence makes the pitch
%ambiguous. The two tones are half and octave apart (fundamentals
%differ by root(2)). Some people will hear the second tone as higher 
%than the first whereas others will hear it lower. 
ynote1=exp(-((log(w)-log(4600))^2)/2)*y1+exp(-((log(2*w)-log(4600))^2)/2)*y2+exp(-((log(4*w)-log(4600))^2)/2)*y4+exp(-((log(8*w)-log(4600))^2)/2)*y8+exp(-((log(16*w)-log(4600))^2)/2)*y16+exp(-((log(32*w)-log(4600))^2)/2)*y32+exp(-((log(64*w)-log(4600))^2)/2)*y64;
ynote1=(ynote1/max(ynote1))*0.5;
pause
sound(ynote1,Fs)
pause
w=2*pi*100*sqrt(2);
y1=sin(w*t);
y2=sin(2*w*t);
y4=sin(4*w*t);
y8=sin(8*w*t);
y16=sin(16*w*t);
y32=sin(32*w*t);
y64=sin(64*w*t);
ynote2=exp(-((log(w)-log(4600))^2)/2)*y1+exp(-((log(2*w)-log(4600))^2)/2)*y2+exp(-((log(4*w)-log(4600))^2)/2)*y4+exp(-((log(8*w)-log(4600))^2)/2)*y8+exp(-((log(16*w)-log(4600))^2)/2)*y16+exp(-((log(32*w)-log(4600))^2)/2)*y32+exp(-((log(64*w)-log(4600))^2)/2)*y64;
ynote2=(ynote2/max(ynote2))*0.5;
pause
sound(ynote2,Fs)
pause
fy1=fft(ynote1);
fy2=fft(ynote2);
plot(abs(fy1(1:10000)))
hold
pause
plot(abs(fy2(1:10000)),'r--')
hold
