%Parameters
FM = 1000;
AmpFM = 2.5;
Fc = 10000;
AmpFc = 2.5;
Kf = 80;
Kp = pi;
SampRate = 1000000;
Wm = 2*pi*FM;
Time = -1/FM:1/SampRate:1/FM;
betaValues = [1 2 5 10 20];

%At constant Wm
for i = 1:5
    SineSignal = AmpFM*sin(Wm*Time);
    Kfc = (Wm*betaValues(i))/AmpFM;
    OutSinePM = AmpFM*cos(2*pi*Fc*Time+Kp*SineSignal);
    IntegralSine = cumsum(OutSinePM)*(1/SampRate);
    OutFM = AmpFM*cos(2*pi*Fc*Time+Kfc*IntegralSine);
    OutFM_S = fftshift(fft(OutFM));
    figure(i);
    plot(Time,abs(OutFM_S));
    title(['Constant Wm / Beta = ' num2str(betaValues(i))]);
end


%At constant Delta W
for i = 1:5
    WmVar = (Kf*AmpFM)/betaValues(i);
    SineSignalVar = AmpFM*sin(WmVar*Time);
    OutSinePMVar = AmpFM*cos(2*pi*Fc*Time+Kp*SineSignalVar);
    IntegralSineVar = cumsum(OutSinePMVar)*(1/SampRate);
    OutFMVar = AmpFM*cos(2*pi*Fc*Time+Kf*IntegralSineVar);
    OutFM_SVar = fftshift(fft(OutFMVar));
    figure(i+5);
    plot(Time,abs(OutFM_SVar));
    title(['Constant Delta W / Beta = ' num2str(betaValues(i))]);
end