% 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);

fc = 100 ;
fm1 = fc * 3 ; 
fm2 = fc * 1.6 ;
tau_amp = .75 ;
attack_exp = 0.1 ;
tau_freq = 1 ;
freq_depth1 = 2.5;
freq_depth2 = 0.5;

a = (t.^attack_exp) .* exp(-t/tau_amp) ;
I1 = freq_depth1 * exp(-t/tau_freq) ;
I2 = freq_depth2 * exp(-t/tau_freq) ;

for i=1:Fs*tmax
    out(i) = a(i) * cos(2*pi*fc*t(i) ...
                        +I1(i)*sin(2*pi*fm1*t(i)) ...
                        +I2(i)*sin(2*pi*fm2*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;

