% see
% http://web.eecs.utk.edu/~qi/ece505/project/proj1.pdf
% info on FM synthesis

%   C2   65.4
%   C3   130.8
%   C4   262 Hz (middle C)   
%   D4   294
%	E4   330
%	F4   349
%	G4   392
%	A4   440  
%	B4   494
%	C5   523  
% see: http://www.phy.mtu.edu/~suits/notefreqs.html
clear all
% 
% sample freq and time to simulate
Fs = 8000;
tmax = 3;
t = 0:1/Fs:tmax ;
out = zeros(1,Fs*tmax);

f1 = 262 ;
f2 = 262 * 2 ; 
f3 = 262 * 3 ;
f4 = 131 * 2.24 ;
f5 = 131 * 2.55 ;

amp1 = 1; fall_time1 = .7 ; attack_time1 = 0.01 ;
amp2 = .3; fall_time2 = .3 ; attack_time2 = 0.01 ;
amp3 = .1; fall_time3 = .1 ; attack_time3 = 0.01;
amp4 = 0; fall_time4 = .4 ; attack_time4 = 0.01 ;
amp5 = 0; fall_time5 = .2 ; attack_time5 = 0.01 ;

a1 = amp1 * (t.^attack_time1) .* exp(-t/fall_time1) ;
a2 = amp2 * (t.^attack_time2) .* exp(-t/fall_time2) ;
a3 = amp3 * (t.^attack_time3) .* exp(-t/fall_time3) ;
a4 = amp4 * (t.^attack_time4) .* exp(-t/fall_time4) ;
a5 = amp5 * (t.^attack_time5) .* exp(-t/fall_time5) ;

for i=1:Fs*tmax
    out(i) = a1(i) * cos(2*pi*f1*t(i)) ...
             + a2(i) * cos(2*pi*f2*t(i)) ...         
             + a3(i) * cos(2*pi*f3*t(i)) ...  
             + a4(i) * cos(2*pi*f4*t(i)) ...  
             + a5(i) * cos(2*pi*f5*t(i)) ...  
         ;
end

soundsc(out,Fs)

figure(2)
plot(out)

figure(1); clf
pwelch(out,[],[],[],Fs)
%set(gca, 'xlim', [0 .2])

figure(3); clf
spectrogram(out,[1024],[],[],Fs, 'yaxis')

% chime
% fc = 131 ;
% fm1 = fc * 1.5 ; 
% fm2 = fc * 5 ;
% tau_amp = .75 ;
% attack_exp = .1 ;
% tau_freq = 1 ;
% freq_depth1 = 0.5;
% freq_depth2 = 1.5;

