%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SYN - Náměty pro 3. cvičení:
%       Aditivní syntéza
%       8. října, 2025
%       Roman Čmejla
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Skript vyžaduje soubor: banjo.wav
%               

%% -----------------------------------------------------------
% PŘÍKLAD 3.1 – Základní periodické průběhy a vliv fáze
%   Postupně přidávejte harmonické složky a sledujte,
%   jak se mění tvar a zvuk výsledného signálu.
%   Porovnejte stejný průběh s náhodnými fázemi.
%   Cíl: ukázat, že lidské ucho je citlivější na složení
%   spektra než na fázi.
%
% a) Sled impulsů
% b) Obdélníkový průběh
% c) Pilový průběh

% a) Sled impulsů
figure
clear, close all
f = 440; fs = 16000; duration = 0.5; t = 0:1/fs:duration;

snd_1a(1,:) = cos(2*pi*f*t); 
snd_1b(1,:) = cos(2*pi*f*t + 2*pi*rand);
for k = 2:10
    snd_1a(k,:) = cos(k*2*pi*f*t);            
    snd_1b(k,:) = cos(k*2*pi*f*t + 2*pi*rand);
    subplot(321), plot(t(1:200), sum(snd_1a(:,1:200)));
    subplot(322), plot(t(1:200), sum(snd_1b(:,1:200)));
    soundsc(sum(snd_1a), fs); pause(1.5*duration)
    soundsc(sum(snd_1b), fs); pause(1.5*duration)
end;
pause
% b) Obdélníkový průběh
snd_2a(1,:) = sin(2*pi*f*t);
snd_2b(1,:) = sin(2*pi*f*t + 2*pi*rand);
for k = 3:2:18
    snd_2a(k,:) = (1/k) * sin(k*2*pi*f*t);     
    snd_2b(k,:) = (1/k) * sin(k*2*pi*f*t + 2*pi*rand);
    subplot(323), plot(t(1:200), sum(snd_2a(:,1:200)));
    subplot(324), plot(t(1:200), sum(snd_2b(:,1:200)));
    soundsc(sum(snd_2a), fs); pause(1.2*duration)
    soundsc(sum(snd_2b), fs); pause(1.2*duration)
end;
pause
% c) Pilový průběh
snd_3a(1,:) = sin(2*pi*f*t); 
snd_3b(1,:) = sin(2*pi*f*t + 2*pi*rand);
for k = 2:18
    snd_3a(k,:) = (1/k) * sin(k*2*pi*f*t);  
    snd_3b(k,:) = (1/k) * sin(k*2*pi*f*t + 2*pi*rand);
    subplot(325), plot(t(1:200), sum(snd_3a(:,1:200)))
    subplot(326), plot(t(1:200), sum(snd_3b(:,1:200)))
    soundsc(sum(snd_3a), fs), pause(1.2*duration)
    soundsc(sum(snd_3b), fs), pause(1.2*duration)
end;

%% -----------------------------------------------------------
% PŘÍKLAD 3.2 – Modelování barvy zvuku
%   Aditivní syntéza umožňuje vytvořit barvu zvuku nástroje
%   podle rozložení jeho harmonických.
%   Flétna / dřevěné nástroje → převaha lichých harmonických.
%   Žestě → výraznější sudé harmonické.
%   Úkol: vygenerujte krátkou stupnici pro oba případy.
%

figure
c1 = 261.63; d1 = 293.66; e1 = 329.63; 
f1 = 349.23; g1 = 392; a1 = 440;  
h1 = 493.88; c2 = 523.25;

notes = [c1 d1 e1 f1 g1 a1 h1 c2]; % melody
duration = [.5 .5 .5 .5 .5 .5 .5 .75]; % note duration [s]

signal = []; % variable initialization
brass = []; flute = [];

fs = 11025; % sampling frequency

for i = 1:length(notes)
    brass = [brass 0.2*sin(1*2*pi*notes(i)*(0:1/fs:duration(i))) + ... 
        0.6*sin(2*2*pi*notes(i)*(0:1/fs:duration(i))) + ...
        0.4*sin(3*2*pi*notes(i)*(0:1/fs:duration(i))) + ...
        0.6*sin(4*2*pi*notes(i)*(0:1/fs:duration(i))) + ...
        0.4*sin(5*2*pi*notes(i)*(0:1/fs:duration(i)))];
    
    flute = [flute 0.8*sin(1*2*pi*notes(i)*(0:1/fs:duration(i))) + ...
        0.4*sin(3*2*pi*notes(i)*(0:1/fs:duration(i))) + ...
        0.2*sin(5*2*pi*notes(i)*(0:1/fs:duration(i)))];
end

sound([brass flute], fs)

subplot(211), plot(flute(1:499))
subplot(212), plot(brass(1:499))

%% -----------------------------------------------------------
% PŘÍKLAD 3.3 – Harmonická analýza v MATLABu
%   Analyzujte reálný zvuk (např. „banjo.wav“).
%   Pomocí Fourierovy transformace zobrazte amplitudové spektrum.
%   Identifikujte dominantní harmonické složky.
% 
% >> analyza('banjo')
% >> [X,Y]=ginput(10)

% function analyza(soubor)
% % funkce analyza(soubor) vykresli amplitudove
% % spektrum *.wav souboru.
% [signal,fs] = wavread(soubor);    
[signal,fs] = audioread('banjo.wav');
N = length(signal);               
c = fft(signal)/N;                    
A = 2*abs(c(2:floor(N/2)));         
f = (1:floor(N/2)-1)*fs/N;         
plot(f,A,'r')
figure
semilogy(f,A,'r')

%% -----------------------------------------------------------
% PŘÍKLAD 3.4 – Aditivní syntéza banja
%   Vytvořte součet několika harmonických složek
%   s amplitudami odpovídajícími reálnému spektru banja.
%   Použijte exponenciální obálku pro doznívání tónu.
%   Porovnejte s originálním zvukem (banjo.wav).
% viz přednáška

figure
% fs   = 16000;
fs   = 11127;
doba = 5;
tau  = .1;
f0   = 392;
nT   = 0:1/fs:doba-1/fs;
ampl = [0.0051 0.0144 0.0246 0.0290 0.0422 0.0101 0.0052 0.0096 0.0022 0.0015];
o=exp(-nT./tau);
x=o.*[ampl*sin(2*pi*[1:length(ampl)]'*f0*nT)];
soundsc(x,fs)
plot(nT,x), title('banjo'), 
axis tight, xlabel('---> cas [s]')


%% -----------------------------------------------------------
% PŘÍKLAD 3.5 – Aditivní syntéza automobilových zvuků
%   a) Volnoběh motoru – harmonické brumění.
%   b) Hluk pneumatik – součet nízkofrekvenčních složek.
%   Cíl: ukázat využití aditivní syntézy i mimo hudbu.

%   a) Volnoběh motoru – harmonické brumění
fs = 16000;     % vzorkovací frekvence
T  = 2;         % délka signálu [s]
t  = 0:1/fs:T-1/fs;
f0 = 55; A = 1./(1:10); y=zeros(size(t));
for k=1:10, y=y+A(k)*sin(2*pi*k*f0*t); end
y=y.*(1+0.05*sin(2*pi*0.7*t)); 
figure, subplot(211), plot(t,y), title('Motor na volnoběh (harmonické)'), xlim([0 0.2])
subplot(212), spectrogram(y,256,200,1024,fs,'yaxis'); 
soundsc(y,fs); pause(T+0.5);

pause(2)

%   b) Hluk pneumatik – "brumění od silnice"
fs = 16000;     % vzorkovací frekvence
T  = 2;         % délka signálu [s]
t  = 0:1/fs:T-1/fs;
f0=70; A=[.5 .25 .12 .06]; y=zeros(size(t));
for k=1:numel(A), y=y+A(k)*sin(2*pi*k*f0*t); end
figure, subplot(211), plot(y(1:2000)), title('Hluk pneumatik')
subplot(212), spectrogram(y,256,200,1024,fs,'yaxis'); ylim([0 1])
soundsc(0.6*y,fs); pause(T+0.5);


%% -----------------------------------------------------------
% PŘÍKLAD 3.6 – Aditivní syntéza ptačího zpěvu
%   a) Výr velký („hu-hůů“) – dvě slabiky s různou délkou a výškou.
%   b) Červenka / drozd – melodické pískání 2–5 kHz.
%   c) Kukačka – dvě slabiky (interval malé tercie), periodické opakování.

% Výr velký („hu-hůů“)
% první slabika: krátký hluboký tón s rychlým náběhem a útlumem
% druhá slabika: delší, doznívající, často o trochu nižší frekvencí

fs = 16000; 
t1 = 0:1/fs:0.4;   % 0.4 s první slabika
t2 = 0:1/fs:0.7;   % 0.7 s druhá slabika

f1 = 250;  % první tón (Hz)
% f2 = 220;  % druhý tón (trochu nižší)
f2 = 245;  % druhý tón (trochu nižší)

env1 = (1-exp(-t1/0.05)).*exp(-t1/0.3);   % hu
env2 = (1-exp(-t2/0.05)).*exp(-t2/0.7);   % hůůů

y1 = sin(2*pi*f1*t1).*env1;
y2 = sin(2*pi*f2*t2).*env2;

owl = [y1 zeros(1,fs*0.1) y2];
soundsc(owl,fs);

figure, 
subplot(211), plot(owl), title('výr'),
subplot(212), spectrogram(owl,256,200,1024,fs,'yaxis'); ylim([0 1])

 
pause(3)
%% Červenka / drozd – krátké melodické flétnové „pískání“
% několik sinusů v rozsahu 2–5 kHz
%  krátké motivy

fs=16000; t=0:1/fs:0.15; 
freqs = [3000 3500 4200]; 
amps = [1 0.5 0.3];
env = (1-exp(-t/0.01)).*exp(-t/0.05);
y = sum(amps'.*sin(2*pi*freqs'*t),1).*env;
song = repmat([y zeros(size(y))],1,5);
soundsc(song,fs);

figure, 
subplot(211), plot(song), title('drozd'),
subplot(212), spectrogram(song,256,200,1024,fs,'yaxis'); ylim([0 8])


pause(3)

%% Kukání kukačky
% dvě slabiky (interval malé tercie: f a d)
% pravidelně se opakují

fs=16000; t=0:1/fs:0.25; 
env = (1-exp(-t/0.02)).*exp(-t/0.3);
f_low = 600; f_high=480; 

low = sin(2*pi*f_low*t).*env; 
high= sin(2*pi*f_high*t).*env; 

% cuckoo = repmat([low zeros(1,fs*0.05) high zeros(1,fs*0.3)],1,5);
cuckoo = repmat([low zeros(1,fs*0.1) high zeros(1,fs*0.8)],1,5);
soundsc(cuckoo,fs);

figure, 
subplot(211), plot(cuckoo), title('kukačka'),
subplot(212), spectrogram(cuckoo,256,200,1024,fs,'yaxis'); ylim([0 2])

%% -----------------------------------------------------------
% PŘÍKLAD 3.7 – Analýza a tvarování obálky
%   Obálka určuje dynamiku zvuku:
%     obálka(t) = A * t^n * exp(-t/τ)
%   Experimentujte s různými hodnotami n (náběh) a τ (útlum).
%   Sledujte vliv tvaru obálky na výsledný zvuk.

% Envelope Type: Attack - Decay = A * t^n * exp(-t/tau)
% n   = 0.01; 0.1; 0.7;
% tau = 0.4;  0.1; 0.7;

figure
% Envelope Analysis  
fs = 8000;   % Sampling frequency
tmax = 3;    % Maximum time duration
t = 0:1/fs:tmax-1/fs;  % Time vector

% Parameters for different envelopes
% n  = [0.75 .75 .75 .75 .75]';
% tau= [.7  .3  .7  .7  .7]';

n  = 0.01;
tau = 0.2;

% Plot the first envelope (n = 0.01, tau = 0.2)
subplot(331)
    plot(t, t.^n)
    ylabel('t^n')
    title('Attack (n) = 0.01')

subplot(334)
    plot(t, exp(-t/tau))
    ylabel('exp(-t/\tau)')
    title('Decay (tau) = 0.2')

subplot(337)
    plot(t, (t.^n) .* exp(-t./tau))
    ylabel('t^n * exp(-t/\tau)')
    xlabel('Time [s]')
    title('Envelope')

n  = 0.5;
tau = 0.5;

% Plot the second envelope (n = 0.5, tau = 0.5)
subplot(332)
    plot(t, t.^n)
    ylabel('t^n')
    title('Attack (n) = 0.5')

subplot(335)
    plot(t, exp(-t/tau))
    ylabel('exp(-t/\tau)')
    title('Decay (tau) = 0.5')

subplot(338)
    plot(t, (t.^n) .* exp(-t./tau))
    ylabel('t^n * exp(-t/\tau)')
    xlabel('Time [s]')
    title('Envelope')

n  = 1;
tau = 1;

% Plot the third envelope (n = 1, tau = 1)
subplot(333)
    plot(t, t.^n)
    ylabel('t^n')
    title('Attack (n) = 1')

subplot(336)
    plot(t, exp(-t/tau))
    ylabel('exp(-t/\tau)')
    title('Decay (tau) = 1')

subplot(339)
    plot(t, (t.^n) .* exp(-t./tau))
    ylabel('t^n * exp(-t/\tau)')
    xlabel('Time [s]')
    title('Envelope')

%% -----------------------------------------------------------
% PŘÍKLAD 3.8 – Aditivní syntéza hudebních nástrojů (harmonické)
%   Vytvořte syntézu pro různé nástroje:
%   a) Klarinet – liché harmonické, pomalý dozvuk
%   b) Drnknutí struny – rychlý náběh, krátký útlum
%   c) Flétna 1 / Flétna 2 – různé poměry harmonických
%   d) Trubka – silný základ, výrazné vyšší harmonické
%   Porovnejte časové průběhy a spektra.
% 
%   a) Klarinet 
%     Frequency component = k * f0
%     Amplitude envelope   = A * t^n * exp(-t/tau)
%     f0 = 262 Hz
%     ------------------------------------
%     k       | A       | n       | tau
%     ------------------------------------
%     1.0000  | 1.0000  | 0.7500  | 0.7000
%     2.0000  |      0  | 0.7500  | 0.3000
%     3.0000  | 0.5000  | 0.7500  | 0.7000
%     4.0000  |      0  | 0.7500  | 0.7000
%     5.0000  | 0.2000  | 0.7500  | 0.7000

clear all, 
fs = 16000;
tmax = 3;
t = 0:1/fs:tmax-1/fs;
f0 = 262 ;
k  = [1   2   3   4   5]';
A  = [1   0  .5   0  .2]';
n  = [.75 .75 .75 .75 .75]';
tau= [.7  .3  .7  .7  .7]';
%  obal = (A*t).^(n*t.^0).*exp((-1./tau)*t);
obal = (A*t).^n.*exp((-1./tau)*t);
% obal = A.*t.^(n*t.^0).*exp((-1./tau)*t);
x=sum([obal.*sin(2*pi*k*f0*t)]);
soundsc(x,fs)
figure
subplot(221)
    plot(t,x), title('klarinet')
subplot(222); 
    plot(t,obal)
    legend('1.','2.','3.','4.','5.')
    ylabel('amplituda frekvencni slozky')
    xlabel('---> cas [s]')
subplot(223);
    N = length(x); c = fft(x)/N; H = 2*abs(c(2:floor(N/2)));         
    f = (1:floor(N/2)-1)*fs/N; plot(f,H,'r')
subplot(224);
    spectrogram(x,[1024],[],[],fs, 'yaxis')


%% b) Drnknutí struny
%     f0 = 262 Hz
%     Amplitude envelope   = A * t^n * exp(-t/tau)
%     ------------------------------------
%     k       | A       | n       | tau
%     ------------------------------------
%     1.0000  | 1.0000  | 0.0100  | 0.7000
%     2.0000  | 1.0000  | 0.0100  | 0.3000
%     3.0000  | 0.5000  | 0.0100  | 0.1000
%     4.0000  | 0.3000  | 0.0100  | 0.1000
%     5.0000  | 0.1000  | 0.0100  | 0.1000

fs = 8000;
tmax = 3;
t = 0:1/fs:tmax-1/fs;

f0 = 262 ;
k  = [1   2   3   4   5]';
A  = [1   1  .5  .3  .1]';
n  = [.01 .01 .01 .01 .01]';
tau= [.7  .3  .1  .1  .1]';

obal = (A*t).^(n*t.^0).*exp((-1./tau)*t);
x=sum([obal.*sin(2*pi*k*f0*t)]);
soundsc(x,fs)

figure
subplot(221)
    plot(t,x), title('drnknutí struny')
subplot(222); 
    plot(t,obal)
    legend('1.','2.','3.','4.','5.')
    ylabel('amplituda frekvencni slozky')
    xlabel('---> cas [s]')
subplot(223);
    N = length(x); c = fft(x)/N; H = 2*abs(c(2:floor(N/2)));         
    f = (1:floor(N/2)-1)*fs/N; plot(f,H,'r')
subplot(224);
    spectrogram(x,[1024],[],[],fs, 'yaxis')  

%% c1) Flétna 1
%     f0 = 440 Hz         T = 1 s
%     f  = f0*2^(3/12) Hz T = 1 s
%     Frequency components = k * f0
%     Envelopes        = A * interp1(X,Y,t)
%                      X = [0 .2 .9 1]
%                      Y = [0 1  .9 0]
%     -------------------
%     k       | A       |
%     -------------------
%     1.0000  | 1.0000  |
%     2.0000  | 0.2500  |
%     3.0000  | 0.0625  |
%     4.0000  | 0.0156  |
%     5.0000  | 0.0625  |
%     6.0000  | 0.0156  |

fs = 16000;
tmax = 3;
t = 0:1/fs:tmax-1/fs;

f0 = 440;
% f0 = 440*2^(3/12);
% f0 = 440*2^(7/12);

k  = [1   2   3     4     5     6]';
%
A  = [1 .25 .0625 .0156 .0625 .0156]';   % 1. flétna

X =  [0 .2 .9 1];
Y =  [0 1  .9 0];
obal = interp1(X,Y,t./tmax);
x = obal.*(A'*sin(2*pi*k*f0*t));
soundsc(x,fs)
subplot(221)
    plot(t,x), title('flétna A')
subplot(222); 
    plot(t,obal)
    legend('1.','2.','3.','4.','5.')
    ylabel('amplituda frekvencni slozky')
    xlabel('---> cas [s]')
subplot(223);
    N = length(x); c = fft(x)/N; H = 2*abs(c(2:floor(N/2)));         
    f = (1:floor(N/2)-1)*fs/N; plot(f,H,'r')
subplot(224);
    spectrogram(x,[1024],[],[],fs, 'yaxis')      

%% c2) Flétna 2
%     f0 = 440 Hz         T = 1 s
%     f  = f0*2^(7/12) Hz T = 1 s
%     Frequency components = k * f0
%     Envelopes        = A * interp1(X,Y,t)
%                      X = [0 .2 .9 1]
%                      Y = [0 1  .9 0]
%     ------------------
%     k       | A      |
%     ----------------------------
%     1.0000  | 1.0000 |
%     2.0000  | 0.0500 |
%     3.0000  | 0.0250 |
%     4.0000  | 0.0005 |
%     5.0000  | 0.0010 |
%     6.0000  | 0.0005 |

fs = 16000;
tmax = 3;
t = 0:1/fs:tmax-1/fs;

f0 = 440;
% f0 = 440*2^(3/12);
% f0 = 440*2^(7/12);

k  = [1   2   3     4     5     6]';
%
A  = [1 .05 .025 .0005 0.001 0.0005]';  % 2. flétna

X =  [0 .2 .9 1];
Y =  [0 1  .9 0];
obal = interp1(X,Y,t./tmax);
x = obal.*(A'*sin(2*pi*k*f0*t));
soundsc(x,fs)

subplot(221)
    plot(t,x), title('flétna B')
subplot(222); 
    plot(t,obal)
    legend('1.','2.','3.','4.','5.')
    ylabel('amplituda frekvencni slozky')
    xlabel('---> cas [s]')
subplot(223);
    N = length(x); c = fft(x)/N; H = 2*abs(c(2:floor(N/2)));         
    f = (1:floor(N/2)-1)*fs/N; plot(f,H,'r')
subplot(224);
    spectrogram(x,[1024],[],[],fs, 'yaxis')      

%% d) Trubka
%     f0 = 440 Hz         T = 3 s
%     Frequency components = k * f0
%     Envelopes            = A * interp1(X,Y,t)
%                            Envelope Y  = [0 1 .8 0]
%     ------------------------------------------------
%     k       | A       |     Envelope X             |
%     ------------------------------------------------
%     1.0000  | 1.0000  | 0    0.0071    0.9875    1 |
%     2.0000  | 1.8222  | 0    0.0286    0.9750    1 |
%     3.0000  | 1.4444  | 0    0.0643    0.9625    1 |
%     4.0000  | 0.5556  | 0    0.1143    0.9500    1 |
%     5.0000  | 0.2222  | 0    0.1786    0.9375    1 |
%     6.0000  | 0.0889  | 0    0.2571    0.9250    1 |
%     7.0000  | 0.0667  | 0    0.3500    0.9125    1 |
%     8.0000  | 0.0444  | 0    0.4571    0.9000    1 |
%     9.0000  | 0.0222  | 0    0.5786    0.8875    1 |
%     ------------------------------------------------

clear
% TRUBKA
fs = 16000;
tmax = 2;
t = 0:1/fs:tmax-1/fs;
f0 = 440 ;
A  = [ 1.0000    1.8222    1.4444    0.5556    0.2222 ...
       0.0889    0.0667    0.0444    0.0222]';
for n=1:length(A);
    X  = [0 .05*n^2/7 1-.05*n/4 1];
    Y  = [0 1 .8 0];
    obal(n,:) = A(n)*interp1(X,Y,t/t(end));
end;
sig=obal.*sin(2*pi*[1:9]'*f0*t);
x=sum([obal.*sig]);
x=x./max(abs(x));
sound(.9*x,fs)

figure
subplot(221)
    plot(t,x), title('trubka')
subplot(222); 
    plot(t,obal)
    legend('1.','2.','3.','4.','5.','6.','7.','8.','9.')
    ylabel('amplituda frekvencni slozky')
    xlabel('---> cas [s]')
subplot(223);
    N = length(x); c = fft(x)/N; H = 2*abs(c(2:floor(N/2)));         
    f = (1:floor(N/2)-1)*fs/N; plot(f,H,'r')
subplot(224);
    spectrogram(x,[1024],[],[],fs, 'yaxis')  

%%
% Obálka Rissetovy trumpety
figure
for n=1:length(A);
    X  = [0 .05*n^2/7 1-.05*n/4 1];
    Y  = [0 1 .8 0];
    obal(n,:) = interp1(X,Y,t/t(end));
    plot(obal(n,:)), hold on
end;
legend('1','2','3','4','5','6','7','8','9')

%% -----------------------------------------------------------
% PŘÍKLAD 3.9 – Aditivní syntéza s neharmonickými parciály
%   Využijte složky, které nejsou přesnými násobky f0:
%   a) Perkuse – krátké inharmonické složky
%   b) Zvon – dlouhé doznívání, inharmonické parciály
%   c) Tympány – kombinace více frekvenčních poměrů a obálek

% a) Bicí
%     f0 = 262 Hz
%     Amplitude envelope   = A * t^n * exp(-t/tau)
%     ------------------------------------
%     k       | A       | n       | tau
%     ------------------------------------
%     1.0000  | 1.0000  | 0.0100  | 0.7000
%     1.5800  | 1.0000  | 0.0100  | 0.3000
%     3.0000  | 1.0000  | 0.0100  | 0.1000
%     2.2400  | 0.3000  | 0.0100  | 0.4000
%     2.5500  | 0.3000  | 0.0100  | 0.1000

fs = 8000;
tmax = 3;
t = 0:1/fs:tmax-1/fs;
f0 = 262;
k  = [1   1.58   3    2.24   2.55]';
A  = [1   1      1    0.3    0.3 ]';
n  = [.01 .01   .01  .01  .01]';
tau = [.7  .3    .1   .4   .1 ]';

envelope = (A * t).^n .* exp(-t ./ tau);
x = sum(envelope .* sin(2 * pi * k * f0 * t));
soundsc(x, fs);

figure;
subplot(221);
    plot(t, x), title('BICÍ');
subplot(222); 
    plot(t, envelope);
    legend('1', '2', '3', '4', '5');
    ylabel('Amplitude of frequency components');
    xlabel('Time [s]');
subplot(223);
    N = length(x); 
    spectrum = fft(x) / N; 
    H = 2 * abs(spectrum(2:floor(N/2)));
    f = (1:floor(N/2) - 1) * fs / N; 
    plot(f, H, 'r');
subplot(224);
    spectrogram(x, [1024], [], [], fs, 'yaxis');


%% b) ZVON (Inharmonic additive synthesis)
%     f0 = 116 Hz         T = 7 s
%     Frequency components = k * f0
%     Envelope            = A * exp(-0.8*t) 
%     ------------------
%     k       | A      |
%     ------------------
%     1.0000  | 1.0000 |
%     2.0000  | 2.7143 |
%     2.3760  | 1.4286 |
%     3.0600  | 0.4286 |
%     4.0000  | 2.0000 |
%     4.9400  | 0.2857 |
%     5.2140  | 0.7143 |
%     5.3000  | 1.0571 |
%     5.9820  | 2.8571 |
%     6.7340  | 0.5143 |
%     8.2740  | 0.8571 |
%     8.9740  | 0.2857 |
%     9.6580  | 0.4286 |
%    10.7700  | 0.8571 |
%    11.7260  | 0.2857 |
%    13.4180  | 0.2857 |
%    16.1540  | 0.1429 |
%    17.0940  | 0.0571 |
%    18.0340  | 0.0286 |
%    19.0600  | 0.1000 |
%    22.0520  | 0.0143 |
%    24.7860  | 0.0429 |
fs = 6000;
tmax = 2;
t = 0:1/fs:tmax-1/fs;
f0 = 110;
k_A =[
    1.0000   1.0000 
    2.0000   2.7143 
    2.3760   1.4286 
    3.0600   0.4286 
    4.0000   2.0000 
    4.9400   0.2857 
    5.2140   0.7143 
    5.3000   1.0571 
    5.9820   2.8571 
    6.7340   0.5143 
    8.2740   0.8571 
    8.9740   0.2857 
    9.6580   0.4286 
   10.7700   0.8571 
   11.7260   0.2857 
   13.4180   0.2857 
   16.1540   0.1429 
   17.0940   0.0571 
   18.0340   0.0286 
   19.0600   0.1000 
   22.0520   0.0143 
   24.7860   0.0429 
 ];

% X =  [0 .01 tmax*.99 tmax];
% Y =  [0 1  1   0];

obal = exp(-2.8*t); %.*interp1(X,Y,t);
x = obal.*(k_A(:,2)'*sin(2*pi*k_A(:,1)*f0*t));
soundsc(x,fs)

 figure
subplot(221)
    plot(t,x), title('zvon')
subplot(222); 
    plot(t,obal)
    legend('1.','2.','3.','4.','5.')
    ylabel('amplituda frekvencni slozky')
    xlabel('---> cas [s]')
subplot(223);
    N = length(x); c = fft(x)/N; H = 2*abs(c(2:floor(N/2)));         
    f = (1:floor(N/2)-1)*fs/N; plot(f,H,'r')
subplot(224);
    spectrogram(x,[1024],[],[],fs, 'yaxis')

    
%% c) Tympány
%     f0 = 55 Hz          T = 2 s
%     Frequency components = k * f0
%     Envelopes            = A * exp(-2.8*t) * interp1(X,Y,t)
%                          X = [0 .002 T*.99 T]
%                          Y = [0 1  .9   0]                          
%     ------------------
%     k       | A      |
%     ------------------
%     1.0000  | 1.0000 |
%     1.5958  | 0.4518 |
%     2.1417  | 0.2740 |
%     2.3000  | 0.2427 |
%     2.6583  | 0.1897 |
%     2.9250  | 0.1613 |
%     3.1625  | 0.1412 |
%     3.5083  | 0.1184 |
%     3.6042  | 0.1131 |
%     3.6542  | 0.1105 |
%     4.0667  | 0.0921 |
%     4.2375  | 0.0859 |
%     4.6083  | 0.0745 |
%     4.8417  | 0.0685 |
%     4.9125  | 0.0668 |
%     5.1417  | 0.0618 |
%     5.4250  | 0.0564 |
%     5.5500  | 0.0543 |
%     5.9875  | 0.0477 |
%     6.1667  | 0.0454 |
%     6.2208  | 0.0447 |
%     6.5417  | 0.0411 |
%     6.7583  | 0.0388 |
%     6.8625  | 0.0378 |
%     7.3417  | 0.0337 |
%     7.4833  | 0.0327 |
%     7.9083  | 0.0297 |
%     8.0875  | 0.0286 |
%     8.6792  | 0.0254 |
%     9.2583  | 0.0227 |

% TYMPÁNY
fs = 16000;
tmax = 4;
t = 0:1/fs:tmax-1/fs;
f0 = 132;
k_A =[
    1.0000   1.0000 
    1.5958   0.4518 
    2.1417   0.2740 
    2.3000   0.2427 
    2.6583   0.1897 
    2.9250   0.1613 
    3.1625   0.1412 
    3.5083   0.1184 
    3.6042   0.1131 
    3.6542   0.1105 
    4.0667   0.0921 
    4.2375   0.0859 
    4.6083   0.0745 
    4.8417   0.0685 
    4.9125   0.0668 
    5.1417   0.0618 
    5.4250   0.0564 
    5.5500   0.0543 
    5.9875   0.0477 
    6.1667   0.0454 
    6.2208   0.0447 
    6.5417   0.0411 
    6.7583   0.0388 
    6.8625   0.0378 
    7.3417   0.0337 
    7.4833   0.0327 
    7.9083   0.0297 
    8.0875   0.0286 
    8.6792   0.0254 
    9.2583   0.0227 ];

X =  [0 .002 tmax*.99 tmax];
Y =  [0 1  .9   0];
obal = exp(-1.8*t).*interp1(X,Y,t);
x = obal.*(k_A(:,2)'*sin(2*pi*k_A(:,1)*f0*t));
soundsc(x,fs)

 figure
subplot(221)
    plot(t,x), title('tympány')
subplot(222); 
    plot(t,obal)
    legend('1.','2.','3.','4.','5.')
    ylabel('amplituda frekvencni slozky')
    xlabel('---> cas [s]')
subplot(223);
    N = length(x); c = fft(x)/N; H = 2*abs(c(2:floor(N/2)));         
    f = (1:floor(N/2)-1)*fs/N; plot(f,H,'r')
subplot(224);
    spectrogram(x,[1024],[],[],fs, 'yaxis')    



%% -----------------------------------------------------------
% PŘÍKLAD 3.10 – Syntéza zvonku
%   Simulace malého zvonku nebo budíku.
%   Kombinace několika harmonických s různou amplitudou a útlumem.
%   Dva zvonky s mírně odlišnou základní frekvencí → vzniká interferenční „chvění“.
%   Zobrazte časový průběh a spektrogram.

clear
fs =44100;          % vzorkovaci frekvence
T1 =0.06;           % doba mezi udery  
T2 =0.480;          % delka posledniho uderu 
 f01=180;           % zakladni frekvence 1.zvonku
 f02=181;           % zakladni frekvence 2.zvonku

A=[.01 .1  1  1];   % amplitudy ctyr oscilatoru
K=[  5 10 20 40];   % nasobky zakl.frekvence                				
                    % jednotlivych oscilatoru 
M=2;                % pocet serii zvoneni
N=4;                % celkovy pocet uderu = 2*N+1

t=0:1/fs:T2-1/fs;
x=[];
for m=1:M            
    for n=1:N        
        x1=A*sin(2*pi*K'*(f01.*t)); % uder 1.zv.
        x1=x1.*exp(-t/T1);  % 1.zvonek s obalkou
        x2=A*sin(2*pi*K'*(f02.*t));% uder 2.zv.
        x2=x2.*exp(-t/T1); % 2.zvonek s obalkou
        x=[x x1(1:T1*fs) x2(1:T1*fs)];   
    end;                                
x=[x x1];  % pripojeni posl.uderu prvniho zvonku 
end;

subplot(211), plot(x)
subplot(212), specgram(x)
soundsc(x,fs)