% 23.10.2019
% 165678 Tomasz Kisielewski
% Zadanie 2 Modulator/demodulator sygnalu GMSK

N = 20; % symobli
M = 20; % sampli
BT = 0.5;

input=randi([0,1],[1,N]);

figure
plot(1:1:N, Wektor_Wejsciowy, 'k*');
title('Modulator - Wektor wejsciowy');
xlabel('Numer probki');
ylabel('Wartosc probki');
axis([1 N -0.5 1.5]);
ax = gca;
ax.YTick = [0 1];

samples = zeros(1,M*N);
for i = 1:1:N
    if input(1,i) == 0
        samples(1,i) = -1;
    end
end

samples = repelem(input,M);

impulse_response_gaussian = zeros(1,2*M+1);
delta = sqrt(log10(2)/(2*pi*BT));
j = 1;

for i = -M:1:M
    t = i/M;
    j = j+1;
    impulse_response_gaussian(1,j) = exp((-t^2)/(2*(delta^2)))/(sqrt(2*pi)*delta);
end

filtered_samples = conv(samples, impulse_response_gaussian, 'same');

figure
plot(1:1:M*N, filtered_samples);
title('Modulator - Filtr Gaussa');
xlabel('Numer probki');
ylabel('Wartosc probki');

integral_samples = zeros(1,N*M);

integral_samples(1,1) = filtered_samples(1,1);
for i = 2:1:N*M
    integral_samples(1,i) = integral_samples(1,i-1) + filtered_samples(1,i);
end

sin_integral_samples = zeros(1,N*M);
cos_integral_samples = zeros(1,N*M);
for i = 1:1:N*M
    sin_integral_samples(1,i) = sin(integral_samples(1,i)/M);
    cos_integral_samples(1,i) = cos(integral_samples(1,i)/M);
end

s = zeros(1, M*N);
c = zeros(1, M*N);
fi = 0;
for i=1:1:N*M
    s(1,i) = sin(fi);
    c(1,i) = cos(fi);
    fi = fi + 0.4*pi;
end

s_samples = sin_integral_samples.*s;
c_samples = cos_integral_samples.*c;
mod_out_samples = s_samples + c_samples;

figure
plot(1:1:M*N, mod_out_samples);
title('Modulator - Wektor Wyjsciowy');
xlabel('Numer probki');
ylabel('Wartosc probki');

% =========================================

d_in_samples = mod_out_samples;
d_s_samples = d_in_samples.*c;
d_c_samples = d_in_samples.*s;

f_lp = zeros(1, N*M);
fi = -20*pi;
for n = 1:1:N*M
    f_lp(1,n) = sin(fi)/fi;
    fi = fi + pi/10;
end

d_f_lp_I = conv(d_s_samples, f_lp, 'same');
d_f_lp_Q = conv(d_c_samples, f_lp, 'same');

figure
plot(1:1:M*N, d_f_lp_I, 'k-', 1:1:M*N, d_f_lp_Q, 'r-');
title('Demodulator - I i Q po filtrze dolnoprzepustowym');
xlabel('Numer probki');
ylabel('Wartosc probki');

d_f_lp_I_diff = zeros(1,N*M);
d_f_lp_Q_diff = zeros(1,N*M);
d_f_lp_I_diff(1,1) = d_f_lp_I(1,1);
d_f_lp_Q_diff(1,1) = d_f_lp_Q(1,1);
for i = 2:1:N*M
    d_f_lp_I_diff(1,i) = d_f_lp_I(1,i) - d_f_lp_I(1,i-1);
    d_f_lp_Q_diff(1,i) = d_f_lp_Q(1,i) - d_f_lp_Q(1,i-1);
end

d_f_lp_I_diff_Q = d_f_lp_I.*d_f_lp_Q_diff;
d_f_lp_Q_diff_I = d_f_lp_Q.*d_f_lp_I_diff;

d_diff_IQ = d_f_lp_I_diff_Q - d_f_lp_Q_diff_I;

figure
plot(1:1:M*N, d_diff_IQ);
title('Demodulator - I i Q zsumowane');
xlabel('Numer probki');
ylabel('Wartosc probki');

d_int_diff = zeros(1, N);

for i=1:1:N
    for j=1:1:M
        d_int_diff(1,i) = d_int_diff(1,i) + d_diff_IQ(1,(i-1)*M + j);
    end
end

d_out = zeros(1, N);
for i = 1:1:N
    if(d_int_diff(1,i) < 0)
        d_out(1,i) =  0;
    elseif(d_int_diff(1,i) > 0)
        d_out(1,i) =  1;
    end
end

figure
plot(1:1:N, d_out, 'k*');
title('Demodulator - Wektor wyjsciowy ');
xlabel('Numer probki');
ylabel('Wartosc probki');
axis([1 N -0.5 1.5]);
ax = gca;
ax.YTick = [0 1];

if(isequal(input, d_out))
    disp('Wektor wejsciowy i wyjsciowy sa rowne');
else
    disp('Odebrany sygnal nie jest nadanym!')
end