%%%%%%%%
% Granular synthesis example
% FOF (fonction d'onde formantique)
%%%%%%%%

%%% base config
fs = 48000;
Ts = 1/fs;
dur = 1;

%%% fundamental frequency -------
f0 = 150;
T0 = 1/f0;

%%% formant frequencies  -------
f1 = 400;
f2 = 2300;

%f1 = 700;
%f2 = 1200;


%%% formant shape -------- 
%%% (as a simplification, we make them common for both formants in this example)

% attack (will be implemented as half the period of one cosine)
a = 0.1;
attack_time = a*T0;

% decay
d = 0.5;
decay_time = d*T0;

%%% grain duration -------
% ( we make it to be enough to hold one RT60 of decay, i.e. 6.91 * decay time )
grain_dur = 6.91*decay_time;


%%% construct grain envelope -------
L_grain = round(grain_dur * fs);
time_axis = (0:L_grain-1)*Ts;

figure;
hold on;

% decay
grain_envelope = exp(-time_axis/decay_time);

plot(time_axis,grain_envelope,'k:');

% attack portion
L_attack = round(attack_time * fs);
cos_half_period = attack_time;
cos_period = 2 * cos_half_period ;
cos_freq = 1/cos_period;
for n = 1:L_attack
    grain_envelope(n) = grain_envelope(n) * 0.5*(1-cos(2*pi*cos_freq*time_axis(n)));
end

plot(time_axis,grain_envelope,'g','linewidth',1);


% now construct a grain ----------

formant1_noEnv = sin(2*pi*f1*time_axis);
formant2_noEnv = sin(2*pi*f2*time_axis);

formant_1_env = grain_envelope.*formant1_noEnv;
formant_2_env = grain_envelope.*formant2_noEnv;


grain = formant_1_env + formant_2_env;

plot(time_axis,formant_1_env,'b--','linewidth',1);
plot(time_axis,formant_2_env,'r--','linewidth',1);
plot(time_axis,grain,'k','linewidth',2);

xlabel('seconds');




% prepare FOF grains

L = round(dur*fs);
N0 = round(fs*T0);
N_grains = floor(L/N0);

x = zeros(L,1);
c = 1;
for i = 1:N_grains
    
    % determine start and end positions in buffer
    grain_start = c;
    grain_end = min(L,c+L_grain);
    
    % add to buffer at current grain start
    k = 1;
    for n = grain_start:grain_end
        if k>L_grain % if our grain is smaller than the fundamental period, stop
            break;
        end
        x(n) = x(n) + grain(k);
        k = k+1;
    end
    
    % advance to next grain
    c = c+N0;
    
end



% plot time-domain
time_axis_tot = (0:L-1)/fs;

figure;
subplot(2,1,1);
plot(time_axis_tot(1:1000),x(1:1000));
xlabel('seconds');


% plot freq-domain
freq_axis_tot = (0:L-1)*fs/L;

xf = fft(x);
subplot(2,1,2);
plot(freq_axis_tot(1:floor(end/2)),db(xf(1:floor(end/2))));
xlabel('Hz');
ylabel('dB');

% play sound
soundsc(x,fs);

