Untitled

%Data: 2023.10.30
%Indeks: 187590
%Imię i nazwisko:   Łukasz Kowalka
%Zadanie 3. Modulator/demodulator sygnałów GMSK
%Komputer nr  226

%Dane wejsciowe
N = 60;
M = 20;
Z = 10; %częstotliwość wartości -1 (im wieksza wartosc tym wieksza częstotliwosc)
% Z z przedzialu (10 do 40) najlepsza wartosc 10

%wygenerować wektor wejsciowy

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

%-----------------------------------------------------------------------
%MODULATOR
%-----------------------------------------------------------------------

%konwesja 0/1 na -1/1

for i=1:N
    if(input(i) == 0)
        input2(i) = -1;
    end
end

figure();
subplot(2,1,1);
stem(input);
title('Wektor wejsciowy 1/0');  
subplot(2,1,2);
stem(input2);
title('Wektor wejsciowy 1/-1');

%nadpróbkowanie

input_probki = zeros(1, N*M);

for i=1:N
    for k=1:M
        input_probki(1,i*M-M+k) = input2(i);
    end
end

figure();
stem(input_probki);
title('Wektor wejsciowy - nadprobkowany M razy');

%odpowiedź impulsowa filtru

% Parametry filtru Gaussa
BT = 0.5;
alpha = sqrt(log10(2)/2/pi/BT);
%alpha = 0.3;

% Czas próbkowania
fs = 10;  % Przykładowa częstotliwość próbkowania (Hz)
T = 1/fs;   % Okres próbkowania

% Zakres czasu
t = -2.55 * alpha : T : 2.7 * alpha;

% Generowanie filtru Gaussa
h = gauspuls(t, alpha, BT);
h = h/Z*2*pi/M;

figure()
stem(h);
title('Odpowiedź impulsowa filtru');

%filtrowanie sygnału
gauss = conv(input_probki, h);

figure();
plot(gauss);
title('Sygnał po filtracji');

%całkowanie

suma_cal = zeros(1,N*M);

for i=2:N*M
    suma_cal(i) = suma_cal(i-1)+gauss(i);
end

figure();
plot(suma_cal);
title('Po calkowaniu');

%Generowanie SIN i COS

% Parametry
numSamples = M*N; % Liczba próbek
numPeriods = 2*N;    % Ilość okresów

% Tworzenie wektora czasu od 0 do 2*pi (1 okres) z określoną liczbą próbek
t = linspace(0, 2*pi*numPeriods, numSamples);

% Generowanie sinusoidy i cosinusoidy

SIN = sin(t);
COS = cos(t);

figure();
subplot(2,1,1);
plot(SIN)
title('SIN');

subplot(2,1,2);
plot(COS);
title('COS');

%generowanie nosnej (sinus i cosinus)

phaseShift = pi/2;

nosna_sin = sin(t + phaseShift);
nosna_cos = cos(t + phaseShift);

%mnozenie przez nosnne

pmod_i = sin(suma_cal);
pmod_q = cos(suma_cal);

GMSK_sin = pmod_i.* nosna_sin;
GMSK_cos = pmod_q.* nosna_cos;

%modulacja

GMSK = GMSK_sin + GMSK_cos;

figure();
subplot(3,1,1);
plot(GMSK_sin)
title('GMSK sin');

subplot(3,1,2);
plot(GMSK_cos);
title('GMSK cos');

subplot(3,1,3);
plot(GMSK);
title('GMSK');

%-----------------------------------------------------------------------
%DEMODULATOR
%-----------------------------------------------------------------------

%mnożenie przez nośne

GMSK_sinD = GMSK.* nosna_sin;
GMSK_cosD = GMSK.* nosna_cos;

figure();
subplot(2,1,1);
plot(GMSK_sinD)
title('GMSK sinD');

subplot(2,1,2);
plot(GMSK_cosD);
title('GMSK cosD');

%generownie odpowiedzi impulsowej filtr dolnoprzepustowego

fc_f_dol = 4; % Częstotliwość odcięcia (0.0 - 1.0, gdzie 1.0 to Nyquist)
N_f_dol = 80; % Długość filtra (nieparzysta liczba)

% Tworzenie wektora próbek czasu
t_f_dol = -N_f_dol/2:N_f_dol/2;

% Tworzenie filtra sinc
FDP = sinc(2 * fc_f_dol * (t_f_dol/N_f_dol));

% Normalizacja filtra
FDP = FDP / 1;

%for i=1:(4*M+1)
%   filtrSinc(i) = sinc((i-2*M-1)/(M/2));
%end

figure();
plot(FDP);
title('filtr DP');

%filtrowanie

I = conv(GMSK_sinD, FDP);
Q = conv(GMSK_cosD, FDP);

figure()
subplot(2,2,1);
plot(I);
title('I');
subplot(2,2,2);
plot(Q);
title('Q');

%różniczkowanie

dIdt = diff(I);
dIdt(end+1) = 0;
dQdt = diff(Q);
dQdt(end+1) = 0;

subplot(2,2,3);
plot(dIdt);
title('dIdt');
subplot(2,2,4);
plot(dQdt);
title('dQdt');

%sygnał różnicowy

skladowa_I = dQdt.*I;
skladowa_Q = dIdt.*Q;
suma_skladowych = skladowa_Q - skladowa_I;

figure()
subplot(3,1,1);
plot(skladowa_I);
title('skladowa I');
subplot(3,1,2);
plot(skladowa_Q);
title('skladowa Q');
subplot(3,1,3);
plot(suma_skladowych);
title('I - Q');

%całkowanie

probkiScalkowaneWy = cumsum(suma_skladowych);
% Liczba próbek, o które chcemy przesunąć sygnał
przesuniecie = 2*M;
% o 40 próbek w prawo (zniwelowanie opóźnienia)
probkiScalkowaneWyy = circshift(probkiScalkowaneWy, -przesuniecie);
probkiScalkowaneWyyy = probkiScalkowaneWyy(1:end-((N_f_dol)-1));

figure();
plot(probkiScalkowaneWyyy);
title('całka sumaD');

%detekcja

sygnWy = zeros(1, N);
stala = 0;
z = 1;
% Detekcja (-1 albo 1)

for i=1:N*M
    if probkiScalkowaneWyyy(i) >= probkiScalkowaneWyyy(i+1)
        stala = stala + 1;
    else
        stala = stala - 1;
    end
    if mod(i, M) == 0 && stala <= 0 
        sygnWy(z) = 1;
        z = z + 1;
        stala = 0;
    end
    if mod(i, M) == 0 && stala > 0
        sygnWy(z) = -1;
        z = z + 1;
        stala = 0;
    end
end

output = zeros(1, N);
% Zamiana -1, 1 na 0/1

for i=1:N
    output(i) = sygnWy(i);
    if output(i) == -1
        output(i) = 0;
    end
end

%porównanie wektorów wejściowego z wyjściowym
figure();
subplot(2,1,1);
stem(input);
title('Wektor wejsciowy');
subplot(2,1,2);
stem(output);
title('Wektor wyjsciowy');