clear all; close all; clc
Fs = 68e3;                    % Sampling frequency
Ts = 1/Fs;                    % Sample time
L = 1e3;                      % Length of signal
t = 0:Ts:L*Ts-Ts;             % Time vector
f = 500;

y = (2^11)/2+round((2^11-10)/2*sin(2*pi*f*t));     % Sinusoid waveform + offset
yy = round((2^11-10)/2*sin(2*pi*f*t));             % Sinusoid waveform without offset

figure(1); plot(1000*t,yy)    % time in miliseconds
title('Sinusový signal 1000Hz')   
xlabel('t [ms]')

Fs2 = Fs*2^11;               % f sampling for PWM has to be resolution times higher
nPWM=1;                      % create PWM signal from sinusoidal waveform
PWM = zeros(1,L*2^11);      
tPWM = 0:1/Fs2:L*Ts - 1/(Fs2);  % time vector

for j = 1:L
    s = y(j);
    for i = 1 : 2^11-1
        if i < s
            PWM(nPWM + i) = 1;
        else
            PWM(nPWM + i) = -1;
        end
    end
    nPWM = nPWM + 2^11-1;
end

figure(2)                           %draw PWM in time base
plot(1000*tPWM,PWM)
title('Signal PWM')
xlabel('time (milliseconds)')
%calculate spectrum
NFFT = 2^nextpow2(L);               % Next power of 2 from length of y
Y = fft(yy,NFFT)/L;
Y_n = 20*log10( 2*abs(Y(1:NFFT/2+1)) );
Y_n = Y_n - max(Y_n);               %the highest value will be 0 dB

f = Fs/2*linspace(0,1,NFFT/2+1);    %frequency axis


figure(3); subplot(211); % Plot single-sided amplitude spectrum.
plot(f*1e-3, Y_n ,'LineWidth',2 ) 
title('Jednostranné amplitudové spektrum 500Hz sinuvého signálu, Fs = 68kHz')
% title('Jednostranné amplitudové spektrum 500Hz sinuvého signálu, Fs = 68kHz')
xlabel('Frekvence (kHz)')
ylabel('20log|Y(f)| - 20log|Y_M(f)| [dB]')
grid on
axis([0 35 -80 2])

LPWM = L*2^11;
NFFT = 2^nextpow2(LPWM);              % Next power of 2 from length of y
Y2 = fft(PWM,NFFT)/LPWM;
f = Fs2/2*linspace(0,1,NFFT/2+1);
Y2_n = 20*log10( 2*abs(Y2(1:NFFT/2+1)) );
Y2_n = Y2_n - max(Y2_n);

subplot(212); plot(f*1e-3,Y2_n,'LineWidth',2) 
title('Jednostranné amplitudové spektrum PWM signálu, Fs = 140MHz')
xlabel('Frekvence (kHz)')
ylabel('20log|Y(f)| - 20log|Y_M(f)| [dB]')
grid on
% axis([0 204 -80 2])
axis([0 80 -80 2])

