%RISSET_HELPER      CONTAINS THE MAIN PARTS OF THE RISSET ILLUSION m-FILE. 
%                   CALLED BY "start demo" button ON RISSET.M GUI.   

if sp~=0;       
    delete(sp);     % to delete old specgram and animation when we start the next demo. look  
end                 % at the beginning and end of risset.m, and at the end of this code where 
sp=0;                    % animation is generated.in this way when graphic objects exist, sp and
if anim~=0;         %anim are handles for them, when they don't exist, they are both 0.
    delete(anim);   %
end                 %
anim=0;

t_planned=12;

for j=100:2:5000;                             % time to reach the next octav
    t_0 = j*2*log(2)/f;                        % Since sin starts from 0, our signal should  
    if t_0>t_planned; break, end                      % also be 0 at t=t_0 forends to match . 
end                                            % By generating t_0 this way, we have the
                                               % parameters of all sine componentsf our vector 
                                               % multiples of 2*pi, which makes our vector 0 at
                                               %t_0.t_0 is very near to planned t, t_planned
                                                    
repeat_num = 3;                             % #of times we repeat the pattern

sampling_rate = 44100;                      % sampling rate  

t=0:1/sampling_rate:(t_0-1/sampling_rate);  % generates the dummy variable vector t which  
                                            % keeps the time, and built so that when we add
                                            % the parts together start and end exactly match
                                            


%f_central = f*2^(n/2); the frequency at the peak of the gaussian.
exp_pow=2;         %power of theexpression in the exponantial 
x=0;
n=round(n_harmonics);
for i=0:n;                                                     
    x=x+ exp(-((abs((i-(t/t_0)-(n/2))/(n*width_coeff))).^exp_pow)).*...% -the envelope(gaussian)
          sin(2*pi*2^i*t_0*f*(1-2.^(-t/t_0))/log(2));          % -harmonics with exponentially
end                                                            % decreasing frequency.

% building the single vector which will be repeated to make the risset 
% scale. Important points:
% - I put the integral of wdt inside the sine
% - Every harmonic reaches the lower one in time t_0
% - Harmonics are saperated and also move in a speed 
%   s.t when their value decreases to their lower 
%   harmonics,their amplitude also becomes the 
%   value of this lower harmonic.
% Overall: my code builds a group of harmonics which 
% move under a "constant" envelope, all with the same 
% speed than I will replicate this and play it. Look 
% at the spectrum if you are confused and you can 
% see how repeating this pattern gives a continous 
% structure
%



temp = repmat(x,1,repeat_num);          % repeating the single vector 

soundsc(temp,sampling_rate)                           % playing the sound

%%%%%%%%%%%%% animation  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
f_start=f;                              % to protect the animation from changing if
width_coeff_start=width_coeff;          % we change parameters during itd playing
axes('position',[.30 .55 .60 .35])      % specifying the axes for the enimation

win_len=3350;                            % it is chosen to make it as near
                                        % as possible ( for generaly used parameters) in
                                        % the computer this code is written. you can find your 
                                        % optimum by trying
for i=1:repeat_num;                     % animation replays itself at the number
                                        % of repititions in risset
    if i                                % 
        for i=0:floor(length(x)/win_len)-1;          % 
            a=i*win_len+1; b=(i+1)*win_len;
            [Pxx,F]=pwelch(x(a:b),win_len,[],2^13,44100);
            plot(log2(F+0.1),...        
                sqrt(Pxx/max(Pxx)),...  
                log2(F+0.1),...         
                exp(-((log2(F+0.1)-log2(f_start+0.1)-(n/2))/(n*width_coeff_start)).^2));
            xlabel('log_2f');
            ylabel('amplitude');
            title('Time evolution of spectrum');
            grid
            drawnow
            % sqrt(Pxx) because Pxx is proportional to square of amplitude.
            %
            %I devided by max(Pxx) to normalize so that amplitudes fits in the gaussian  
            %
            %0.01 is added to F to preventlog2(0)
            %
        end
    end
end
anim=gca;
%%%%%%%%%%%%%%%%%%%%%  spectrogram  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%