clc;clear;
n = 50:70; % rozsah harmonickych, melo by byt symetricke okolo nosne frekvence [kHz]
o = 2*pi*1000; % frekvence obalky
w = 2*pi*60000; % frekvence nosne
pm = 6.4; % amplituda pm
p0 = 2e-005; % prah slysitelnosti
t = linspace(0,4*pi/o,10000); % normalni cas --- pro vykresleni prubehu buzeni
f = 0; % vysledna funkce
harmonics = 1:120; % pole harmonickych
peak = zeros(1,120);

for harm=n
    deviation = small_a(w,o,harm)*sin(harm*o*t);
    peak(harm) = pm*small_a(w,o,harm);
    %speed = n(i)*o*pm*small_a(w,o,n(i))*cos(n(i)*o*t);
    f = f + deviation;

end

Lp = 20*log10(peak/p0); % logaritmus Lp = 20*log(p/p0)

g = sqrt(1-cos(o*t));
h = -g;

%% ----------PLOT------------
figure
hold on
grid on
plot(t,f,t,g,t,h)
xlabel('t [s]')
ylabel('w [-]')
figure
hold on
grid on
bar(harmonics, Lp, 'red')
xlabel('frekvence [kHz]')
ylabel('L_p [dB]')

%% ----------------FUNCTIONS-------------------

function an = big_a(w,o,n)
    an = sqrt(2)*o^2*(-16*n*w*o-(2*w+o+2*n*o)*(2*w-o+2*n*o)+cos(2*pi*(n-w/o))*(2*w-o-2*n*o)*(2*w+o-2*n*o)*cos((2*pi*(w+n*o))/o));
end

function a = small_a(w,o,n)
    a = big_a(w,o,n)/(pi*(-8*(4*n^2+1)*w^2*o^2+(1-4*n^2)^2*o^4+16*w^4));
end