clear all;
Br = 117.1875*10^3;
Fd = Br*8;

Fif = 20*10^6;
Fdif = 60*10^6;
Fdopl = 100;

%% === TRANSMITTER ===
M = pngen(5);
D = rand(1, 200*8) > 0.5;
CRC = zeros(1, 16);
for i=1:length(D),
CRC(16) = xor(CRC(16), D(i));
CRC = [CRC(2:16) CRC(1)];
end;
FILL = ones(1, 200)*0.5;
DATA = [dec2bin(length(D)/8, 11)=='1' D CRC];
packet = [FILL M M diff_enc([0 DATA]) FILL]*2-1;

% redescr
Ndkr = Fd/Br;
pkt = packet(ceil(0.01:1/Ndkr:length(packet)));

hch = firls(31, [0 .24 .3 .5]*2, [1 1 0 0]);
pktf = filter(hch, [1], pkt);

% upsample
t1 = upsample4(pktf, 31);
t2 = upsample4(t1, 31);
pktfu = upsample4(t2, 31)';

% -> IF
l = length(pktfu);
s = sin(2*pi*(Fif+Fdopl)*(0:1/Fdif:l/Fdif));
c = cos(2*pi*(Fif+Fdopl)*(0:1/Fdif:l/Fdif));
pkts = pktfu.*s(1:l);
pktc = pktfu.*c(1:l);

%% === TRANSMIT LINE ===
Ir = awgn(awgn(pkts, 0.05), 0.01);
Qr = awgn(awgn(pktc, 0.05), 0.01);

%% === RECEIVER ===
% -> 0
l = length(Ir);
s = sin(2*pi*Fif*(0:1/Fdif:l/Fdif));
c = cos(2*pi*Fif*(0:1/Fdif:l/Fdif));

Rreu = Ir.*s(1:l) + Qr.*c(1:l);
Rimu = Qr.*s(1:l) - Ir.*c(1:l);

tr00 = downsample4(Rreu);
tr01 = downsample4(tr00);
Rre_ = downsample4(tr01);
Rre = filter(hch, [1], Rre_);

tr10 = downsample4(Rimu);
tr11 = downsample4(tr10);
Rim_ = downsample4(tr11);
Rim = filter(hch, [1], Rim_);

%figure; hold on; plot(Rre); plot(Rim, 'r');

% detection
Ms = M(ceil(0.01:1/Ndkr:length(M)));
cfh = Ms(length(Ms):-1:1);
cre = filter(cfh, [1], Rre)/length(cfh);
cim = filter(cfh, [1], Rim)/length(cfh);
%cr = sqrt(cre.^2 + cim.^2);
cr = (7/8)*max(abs(cre), abs(cim)) + (1/8)*min(abs(cre), abs(cim));

%figure; hold on; plot(cre); plot(cim); plot(cr, 'r');

env(length(cr)) = 0;
for j=1:length(cr),
f = max(1, j-Ndkr*length(M)+2*Ndkr);
t = min(length(cr), max(2, j-Ndkr));
env(j) = sum(cr(f:t))/(t-f);
end;

%figure; hold on; plot(cr, 'r'); plot(env, 'g'); plot(env*5);

%
for i=length(M)*Ndkr+1:length(cr)-1,
c1 = cr(i) > env(i)*5;
% c2 = (cr(i) > cr(i-1)) && (cr(i) > cr(i+1));
c2 = (cr(i-1) > cr(i-2)) && (cr(i-1) > cr(i));
c3 = abs(cr(i)-cr(i-length(M)*Ndkr)) < 4*env(i);
c4 = cr(i-length(M)*Ndkr) > env(i-length(M)*Ndkr)*5;

if c1 && c2 && c3 && c4,
break;
end;
end;
i
re = Rre(i:length(Rre));
im = Rim(i:length(Rim));

%figure; hold on; plot(re); plot(im, 'r');
figure; hold on; plot(cr, 'r'); plot(env, 'g'); plot(env*5); plot((c1&c2&c3&c4)*1.3, 'k');

%
l = length(re);
lb = length(packet);
Yre(lb) = 0;
Yim(lb) = 0;
fre(lb) = 0;
fim(lb) = 0;
are(lb) = 0;
aim(lb) = 0;
mark(l) = 0;
bit(lb) = 0;

k = 1;
gamma = 0.1;
Yre(1) = re(1);
Yim(1) = im(1);

for i = 8:8:length(re)-10,
mark(i) = 1;
are(k) = sum(re(i-6:i-1));
aim(k) = sum(im(i-6:i-1));

% okre(k) = are(k)*Yim(k) - aim(k)*Yre(k);
okim(k) = aim(k)*Yim(k) + are(k)*Yre(k);

bit(k) = sign(okim(k));

fre(k) = are(k)*bit(k);
fim(k) = aim(k)*bit(k);

Yre(k+1) = (1-gamma)*Yre(k) + gamma*fre(k);
Yim(k+1) = (1-gamma)*Yim(k) + gamma*fim(k);

k = k + 1;
end;

bit_dec = diff_dec(sign(bit+1));
bits = bit_dec(2:length(bit_dec));

%figure; hold on; plot(okim); plot(okre, 'r');
%figure; hold on; plot(re); plot(Yre(ceil(0.01:1/Ndkr:length(bit))), 'r');
%figure; hold on; plot(im); plot(Yim(ceil(0.01:1/Ndkr:length(bit))), 'r');

figure; hold on; plot(im); plot(mark, 'k');

%
len = bin2dec(int2str(bits(1:11)));
data = bits(12:12+len*8-1);
rcrc = bits(12+len*8:12+len*8+15);
crc = zeros(1, 16);
for i=1:length(D),
crc(16) = xor(crc(16), data(i));
crc = [crc(2:16) crc(1)];
end;

disp('---');
disp(sprintf('Received packet length: %d bytes', len));
disp(sprintf('Data ok: %d', sum(abs(data-D))==0));
disp(sprintf('CRC ok: %d', sum(abs(crc-rcrc))==0));

clear
% несущаЯ
Fc = 2400;
% частота выборки бит из вектора Х
Fd = 2400;
% частота выборки выходного модулированного сигнала
Fs = 19200;
% количество сэмплов на 1 бит
Sb = Fs/Fd;
% количество 0 сэмплов перед сигналом
ZeroS = 10;%round(Sb/0.8);
% количество 0 сэмплов перед сигналом
NoiseLev = 0.2;
% ограничение тракта
Amax = 5;
% Количество наборов фаз на символ
M = 2;
% МодулирующаЯ последовательность
x = repmat([0 1 0 0 0 0 0 1 1 1 1 1 0 1 0 1 0 1 0 0 1 0 0 0],1,1); x=x(;
% Modulate, keeping track of time.
[y,t] = dmod(x,Fc,Fd,[Fs pi/2],'psk',M);
% Всего выборок в выходом сигнале
samples_total = length(y);

if ZeroS > 0
y(ZeroS + 1: samples_total) = y (1 : samples_total - ZeroS);
y(1 : ZeroS) = 0;
end

noise = NoiseLev*randn(size(y));

% yn сигнал на входе демодулЯтора
yn = y + noise;

% Смоделировать выход за динамический диапазон по входу
for sn = 1 : samples_total
if yn(sn) > Amax
yn(sn) = Amax;
elseif yn(sn) < -Amax
yn(sn) = -Amax;
end
end

% ДемодулЯциЯ
% Вектор номеров выборок
samples_vector = 1 : samples_total;
%ГенерациЯ векторов sin и cos локальной несущей длЯ коррелЯции c Y
dFs = 2*pi*Fc/Fs;
% sin и cos на длине одного периода
sinFc = sin(0:dFs:2*pi-dFs)';
cosFc = cos(0:dFs:2*pi-dFs)';
samples_sin = length(sinFc);
% sin и cos на длине вектора входного сигнала
sinY = repmat(sinFc, samples_total/samples_sin, 1);
cosY = repmat(cosFc, samples_total/samples_sin, 1);

% ДвигаЯсь по выборкам скользЯщим окном, равным длине битового интервала, вычислЯем коррелЯцию с несущей
demod_y(1:samples_total) = 0;
demod_I(1:samples_total) = yn.*cosY;
demod_Q(1:samples_total) = yn.*sinY;
demod_y = demod_I + demod_Q*j;

% Усреднение I и Q на длине битового интервала
demod_avr_I(1:samples_total) = 0;
demod_avr_Q(1:samples_total) = 0;
avr_interval = samples_sin;

for sn = avr_interval : samples_total
demod_avr_I(sn) = sum(demod_I( sn - avr_interval + 1: sn ))/avr_interval*2;
demod_avr_Q(sn) = sum(demod_Q( sn - avr_interval + 1: sn ))/avr_interval*2;
end

% Перемножение усредненных I и Q с задержанной на 1 бит версией
demod_dly_I(1:samples_total) = 0;
demod_dly_Q(1:samples_total) = 0;
demod_dly_I(1+Sb:samples_total) = demod_avr_I(1:samples_total - Sb).* demod_avr_I(1+Sb:samples_total);
demod_dly_Q(1+Sb:samples_total) = demod_avr_Q(1:samples_total - Sb).* demod_avr_Q(1+Sb:samples_total);
demod_dly = demod_dly_I + demod_dly_Q;

% ФильтрациЯ демодулированного сигнала
[demod_b,demod_a]=butter(2,(Fd/Fs));
demod_filt = filter(demod_b,demod_a,sign(demod_dly));


subplot(5,1,1), plot (samples_vector, y, 'b', samples_vector, y, 'r'), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,2), plot (samples_vector, yn), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,3), plot (samples_vector, demod_I, samples_vector, demod_Q, samples_vector, abs(demod_y)), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,4), plot (samples_vector, demod_dly_I, samples_vector, demod_dly_Q, samples_vector, demod_dly), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]');
subplot(5,1,5), plot (samples_vector, demod_dly, samples_vector, sign(demod_dly), samples_vector, sign(demod_filt)), grid on;
set (gca, 'XLimMode', 'manual',...
'XLim', [1,samples_total]',...
'XTickMode', 'manual',...
'XTick', [Sb:Sb:samples_total]',...
'YLimMode', 'manual',...
'YLim', [-1.5, 1.5]',...
'YTickMode', 'manual',...
'YTick', [-1.5:0.5:1.5]');

function x=pngen(p);
% x=pnegn(p) - Generate pseudo-random sequence
% p = either a scaler (2..12) to generate a PN sequence
% using a predefined polynomial of length p, or
% a binary generating polynomial.
% x = sequence, of length 2^length(p) - 1.

% Predefined seqence polynomials:
pl{2} = [1 1];
pl{3} = [1 1 0];
pl{4} = [1 0 0 1];
pl{5} = [1 0 1 1 1];
pl{6} = [1 1 1 0 0 1];
pl{7} = [1 1 0 0 1 0 1];
pl{8} = [1 0 1 0 1 1 1 1];
pl{9} = [1 0 1 0 1 1 0 1 1];
pl{10} = [1 0 1 0 1 0 1 0 1 1];
pl{11} = [1 0 1 0 1 0 1 0 1 0 1];
pl{12} = [1 1 0 1 0 0 0 0 1 1 0 1];
if length(p) == 1
p = pl{p};
end

L=length(p);
N=2^L-1;

s=ones(1,length(p));
x=zeros(1,N);
for i=1:N
snew = mod(sum(and(p,s)),2);
x(i) = snew;
s=[s(2:end),snew];
end

function enc=diff_enc(data, start);
    if ~(exist('start'))
        start = 0;

    de(1)=xor(data(1), start);
    for i=2:length(data), de(i)=xor(de(i-1), data(i)); end;
    enc = de;
end

%% BPSK модель
% Количество бит данных
N = 4;
% Данные
data = [1 0 1 0];
%data = randi(1,N);
pData = data*2 - 1;
% Несущая частота (100 кГц)
f= 100000;
% частота выборок
fs = f*100;
% Период выборки
Ts = 1/fs;
% Период несущей
T = 1/f;
% Количество периодов несущей на бит-символ
M = 4;
% Длина пакета данных
n = M*length(data);
% Текущее время
t = 0:Ts:n*T;
% Несущая в пакете
car = sin(2*pi*f*t);
% График данных
subplot(2, 1, 1);
stem(pData);
title('Данные');
% График несущей
subplot(2, 1, 2);
plot(car);
title('Несущая');

%% Преобразование данных в прямоугольные импульсы
tp = 0:Ts:T*M;
exdata = [];
for (i = 1:length(data))
for(j = 1:length(tp) - 1)
exdata = [exdata pData(i)];
end
end
exdata = [exdata 0];
figure;
plot(exdata, 'r-', 'LineWidth', 4);
hold on;
grid on;
plot(car, 'g-');
hold on

%% Модуляция
% Перенос exdata на несущую
mSig = exdata.*car;
plot(mSig, 'h-', 'LineWidth', 2);
hold off;
title('Прям. импульсы + несущая + промод.несущая');

%% AWGN (аддитивный белый гауссов шум) канал передачи данных
figure;
% отношение сигнал/шум в децибелах
SNR = 50;
% мощность сигнала в децибелах
SIGPOWER = 10;
% смесь мод.сигнала с шумом
rx = awgn( mSig, SNR, SIGPOWER);
plot(mSig, 'r-', 'LineWidth', 3);
hold on;
plot(rx, 'g-', 'LineWidth', 1);
grid on;
%% Демодуляция
% умножение мод.зашумленого сигн.на несущую
dem = rx.*car;
figure;
plot(dem);
grid on;
%% Декодирование
k = 1;
rcv = [];
for(i = 1:length(data))
sm = 0;
for(j = 1:length(tp) - 1)
sm = sm + dem(k);
k = k + 1;
end
if(sm > 0)
rcv = [rcv 1];
else
rcv = [rcv 0];
end
end
clc;
% Переданные данные
data
% Принятые данные
rcv