clc;clear;
n = 40:60; % rozsah harmonickych, melo by byt symetricke okolo nosne frekvence [kHz]
o = 2*pi*1000; % frekvence obalky
w = 2*pi*50000; % frekvence nosne
pm = 6.4; % amplituda pm
global m;
m = 1;
p0 = 2e-005; % prah slysitelnosti
t = linspace(-2*pi/o,2*pi/o,10000);
tau = linspace(-1,1,10000);
f = 0; % vysledna funkce
harmonics = 1:100; % pole harmonickych
peak = zeros(1,100);

for harm=n
    coeff_a = an(o,w,harm)
    coeff_b = bn(o,w,harm)
    deviation = coeff_a*cos(harm*o*t)+coeff_b*sin(harm*o*t);
    peak(harm) = pm*sqrt(coeff_a^2+coeff_b^2);
    f = f + deviation;
end

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

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

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

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

function sum = an(o,w,n)
    global m;
    dt = 0.000001;
    t = 0:dt:2*pi/o;
    sum = 0;
    for i=1:length(t)
        sum = sum + dt*sqrt(1+sin(o*t(i)))*sin(w*t(i))*cos(n*o*t(i));
    end
    sum = sum*o/pi;
end

function sum = bn(o,w,n)
    global m;
    dt = 0.000001;
    t = 0:dt:2*pi/o;
    sum = 0;
    for i=1:length(t)
        sum = sum + dt*sqrt(1+sin(o*t(i)))*sin(w*t(i))*sin(n*o*t(i));
    end
    sum = sum*o/pi;
end