function y = risset_beats(fund,incr,harm,tones,dur,fs,bwr)
%RISSET_BEATS Generates Risset's musical beats
%   y = risset_beats(fund,incr,harm,tones,dur,fsmbwr)
%   fund: the fundamental of the first tone in Hz
%   incr: the increment in Hz for the subsequent tones
%   harm: the number of harmonincs in each tone
%   dur:  the duration of the resulting sound
%   fs:   the sampling frequency
%   bwr:  write / not write jpeg and wav
%
%   Example: risset_beats(100,0.1,7,7,10,8000);
%
%   This function generates a sum of tones each of
%   which has fund+incr frequency and contains
%   harm number of harmonics with 1/f amplitude.
%   It writes WAV file and plots and saves specgram
%
%   First done by Jean-Claude Risset at Bell Labs

% ------- risset_beats.m -----------------------------------
% Spring 2002
% Marios Athineos, marios@ee.columbia.edu
% http://www.ee.columbia.edu/~marios/
% Copyright (c) 2002 by Columbia University.
% All rights reserved.
% ----------------------------------------------------------

% Defaults
if nargin < 7; bwr = 0; end

nr      = fs/2;
nfft    = 1024;
% nWindow is the 70th prime number !
% Think about it ;)
nWindow = 349;
twopi   = 2*pi;

% The time vector
t=[0:1/fs:dur].';
% Initialize the output
y=zeros(size(t));

for j=1:tones
    for i=1:harm
        % Check for aliasing
        if(i*fund>nr)
            break;
        end        
        x=(1/i)*sin(twopi*i*fund*t);
        y=y+x;
    end
    fund=fund+incr;
end

y=scale(y);

figure;
specgram(y,nfft,fs,nWindow,floor(0.75*nWindow));
title('Risset Musical Beats', 'FontWeight', 'bold')
if bwr
    print -djpeg80 risset_beats;
    wavwrite(y,fs,16,'risset_beats.wav');
end
sound(y,fs);

% ----------------------------------------------------------
function y = scale(x)
% Subfunction taken from soundsc, used ofr scaling

% Determine scaling vector, SLIM:
xmin = min(x(:));
xmax = max(x(:));
slim = [xmin xmax];

% Scale the data so that the limits in
% SLIM are scaled to the range [-1 +1]
dx = diff(slim);
if dx==0,
    % Protect against divide-by-zero errors:
    y = zeros(size(x));
else
    y = (x-slim(1))/dx*2-1;
end