Shepard Tone

Roger Shepard created an amusing auditory paradox, which is after him named: Shepard tone. It is based on a self-similar sequence of notes that consist on a superposition of 12 tones, each an octave higher then the lower neighbor. The 12 tones extend from the lower frequency limit of auditory perception to the upper frequency limit of hearing. In the present approach it goes from 10 Hz to 20,480 Hz, with all octave tones interleaved. The set of frequencies composing our Shepard tone is 10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10,240, and 20,480 Hz. In order to create a paradoxal sound of a continouos ingreasing frequency sound, we create a sweep, an exponential chirp for each composing tone, going from its initial frequency to its following octave.

The instantaneous frequency of each tone is given by
$f(t) = f_0 k^t$.
As we want the final frequency after $T$ seconds to be $2 f_0$,
$f(T) = 2 f_0 = f_0 k^T$
That leads to
$k = e^{\frac{\ln 2}{T}}$
and the instantaneous frequency may be written as
$f(t) = f_0 e^{\frac{\ln 2}{T} t}$
An exponential chirp tone $x(t)$ is given by
$x(t) = \sin \left( 2 \pi \int_{0}^{t} f(\tau) d\tau \right)$
$x(t) = \sin \left( 2 \pi f_0 \frac{e^{\frac{\ln 2}{T}t}-1}{\ln 2 / T} \right)$

To acchieve a better result I have weighted the tones using a gaussian window, and I created the chirp tone with the double sampling frequency and in the end I made a resample to the desired sampling frequency, to attenuate the aliasing caused by the higher chirping tones.

fs = 44100;
fs2 = 2*fs;
ts = 1/fs2;
F = [10 20 40 80 160 320 640 1280 2560 5120 10240 20480];
w = gausswin(length(F));
x = [];
T = 10;
t = [0:1/fs2:T-1/fs2]';
x = zeros(size(t));
w(end+1)=0;
for i = 1 : length(F),
    x += linspace(w(i),w(i+1),length(t))' .* sin(2*pi*F(i)*(exp(log(2)/T * t)-1)/(log(2)/T));
end;
x = resample(x,1,2);

I used the function bellow to concatenate several chirping tones:

function y = concatenate(varargin)
y = varargin{1};
for i = 2 : nargin,
    x = varargin{i};
    o1 = linspace(1,0,11025)';
    o2 = linspace(0,1,11025)';
    y = [y(1:end-11025); o1.*y(end-11025+1:end)+o2.*x(1:11025); x(11025+1:end)];
end;

And the result is plotted in a Spectogram (see figure bellow) and might also be heard.

y = concatenate(x,x,x);
specgram(y,[],fs,[],[]);
soundsc(y,fs);