%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Multimedia Signal Synthesis at CTU Prague 
% Copyright © 2014 by Roman Cmejla
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Suggestions for 8th laboratory: 
% Modulation synthesis
% November 13, 2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Recommended links: 
% http://music.columbia.edu/cmc/musicandcomputers/chapter4/04_05.php
% http://music.columbia.edu/cmc/musicandcomputers/chapter4/04_07.php
% http://amber.feld.cvut.cz/vyu/zzs/zzs12/cv12.htm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Example 1:  Amplitude and ring modulation
fs = 8000;                    
Ts = 1/fs;                    
Tmax = 1.0;                   
t = [0:Ts:Tmax-Ts];           
x = 0.5*cos(2*pi*750*t) + 0.5*cos(2*pi*450*t);
x = sin(2*pi*600*t).*sin(2*pi*150*t);
%%
x = (1 + 1*sin(2*pi*150*t)).*sin(2*pi*600*t);
%%
mi    = t/Tmax;            % linear increasing of the modulation depth od 0 do 1
x = (1 + mi.*sin(2*pi*150*t)).*sin(2*pi*600*t);
%%
mi    = t/Tmax              
x = (1 + mi.*sin(2*pi*4*t)).*sin(2*pi*400*t);
%%
mi=0.5;
x = (1 + mi.*sin(2*pi*4*t)).*sin(2*pi*400*t).*exp(-t);

% Example 2: Frequency modulation
% Parameters
doba = 2;	               % duration [s]
fs   = 8000;	               % sampling frequency [Hz]
fc   = 500;	               % carrier frequency (Hz)
fm   = 50;                     % modulation frequency (Hz)
                               % 2.5 5 10 25 50 100 250
                               % 4.05 8.09 16.2 40.5 80.9 162 405
mi   = 1;                      % modulation index
t    = 0 : 1/fs : doba-1/fs;   % time vector

% Example 3: Synthetic bells
% BELS
fc   = [110,220,110,110,250,250];             % carrier frequency [Hz]
fm   = [220,440,220,220,350,350];             % modulation frequency [Hz]  
                                              % ratio 1:2
Io   = [10,5,10,10,5,3];                      % modulation index
tau  = [2,2,12,0.3,2,1];                      % time constant
doba = [6,6,3,3,5,5];                         % duration [s]
fs   = 4.*[11025,11025,11025,11025,11025,11025]; % sampling frequency [Hz]

% Example 4:  Synthetic wind 
% WIND %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% frome tone to noise
% Greater "amplitude" (Melusina) reaches
% 1. increasing the modulation index
% 2. increasing the number of frequency components of signals, eg. 500
fmin = 261.6;                 % c1
fmax = 523.3;                 % c2
fc   = linspace((fmin),(fmax),100);
fm   = 0.8;                   % modulation frequency [Hz]  
doba = 4;                     % duration [s]
fs   = 8000;                  % sampling frequency [Hz]
t    = 0:1/fs:doba-1/fs;      % time axis
mi   = 1;                     % modulation index

% Example 5: Simple synthetic FM music instruments
%                        duration fc   H =fm/fc Imax Imin
%          Bell          15       200  1.4      10   0
%          Drum          0.2      80   0.688    25   0
%          Clarinet      0.5      900  0.667     4   2
%          Bassoon       0.5      1100 0.2       1.5 0
%          Brass         0.6      440  1.0       5    0

% Example 6: Warning - alerts
fc   = [1500,1500,1500,1500,1500];          % carrier frequency [Hz]
fm   = [1.25,2,5,5,10,25];                  % modulation frequency [Hz]  
Io   = [400,200,100,50,20];                 % modulation index
doba = [2,2,2,2,2];                         % duration [s]
fs   = [44100,44100,44100,44100,44100];     % sampling frequency [Hz]

% Example 7: Synthetic water drops
fs    = 44100;
doba1 = 0.02;         % duration of first part
doba2 = 0.18;         % duration of second part
t1    = 0:1/fs:doba1-1/fs; 
t2    = 0:1/fs:doba2-1/fs;
t     = 0:1/fs:doba1+doba2-1/fs;
fm    = 7;
mi    = 20;
tau   = 0.02;
pocet_kapek = 5;

% Example 7: Vibrato
fs = 22050;
f0 = 440; 
fm = 5; 
doba = 5;
t = 0:1/fs:doba-1/fs;
alfa = [1000 300 .1 .1 .1 3 .5 1 4];
% vibrato s aditivni syntezou
x = zeros(1,length(t));
for n = 1:9
    x = x + alfa(n)*sin(2*pi*f0*n.*t + 7*sqrt(n)*sin(2*pi*fm.*t));
end
% envelope
X = [0 .2 .8 .95 1].*t(end); Y = [0  1 .3 .1  0]; O=interp1(X,Y,t);