%% BPSK модель
% Количество бит данных
N = 4;
% Данные
data = [1 0 1 0];
%data = randi(1,N);
pData = data*2 - 1;
% Несущая частота (100 кГц)
fn= 100000;
% частота выборок
fns = fn*10;
% Период выборки
Tns = 1/fns;
% Период несущей
Tn = 1/fn;
% Количество периодов несущей на бит-символ
M = 4;
% Длина пакета данных
n = M*length(data);
% Текущее время
tn = 0:Tns:n*Tn;
% Несущая в пакете
car = sin(2*pi*fn*tn);

%% Преобразование данных в прямоугольные импульсы
tpn = 0:Tns:Tn*M;
exdata = [];
for (i = 1:length(data))
for(j = 1:length(tpn) - 1)
exdata = [exdata pData(i)];
end
end
exdata = [exdata 0];

%% Модуляция
% Перенос exdata на несущую
mSig = exdata.*car;

%% BPSK модель
% Количество бит данных
N = 4;
% Данные
data = [1 0 1 0];
%data = randi(1,N);
pData = data*2 - 1;
% Несущая частота (100 кГц)
fn = 100000;
% Доплер
fd = 100;
% частота выборок
fns = fn*8;
% Период выборки
Tns = 1/fns;
% Период несущей
Tn = 1/fn;
% Количество периодов несущей на бит-символ
M = 4;
% Длина пакета данных
n = M*length(data);
% Текущее время
tn = 0:Tns:n*Tn;
% Несущая в пакете
car_sin = sin(2*pi*fn*tn);
car_cos = cos(2*pi*fn*tn);

%% Преобразование данных в прямоугольные импульсы
tpn_i = 0:Tns:Tn*M;
exdata_i = [];
for (i = 1:length(data))
for(j = 1:length(tpn_i) - 1)
exdata_i = [exdata_i pData(i)];
end
end
exdata_i = [exdata_i 0];

tpn_q = 0:Tns:Tn*M;
exdata_q = [];
for (i = 1:length(data))
for(j = 1:length(tpn_q) - 1)
% exdata_q = [exdata_q pData(i)];
if (length(tpn_q) == 17)
exdata_q = [exdata_q 1];
else
exdata_q = [exdata_q -1];
end
end

end
exdata_q = [exdata_q 0];

%% Модуляция
% Перенос exdata на несущую
mSig_i = exdata_i.*car_cos;
mSig_q = exdata_q.*car_sin;

%%

figure;
plot(exdata_i, 'r-', 'LineWidth', 4);
hold on;
grid on;
plot(exdata_q, 'b-', 'LineWidth', 4);
hold on;
plot(car_sin, 'g-');
hold on
hold off;

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

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

%% === TRANSMITTER ===
M = pngen(5);
%M = 0.5;
D = rand(1, 200*8) > 0.5;
%D = [1 0 1 0];
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);

%% === GRAPHICS ===

figure;
plot(pkts, 'g-');
grid on;
hold on
plot(pktc, 'r-');
hold on

figure;
plot(Ir, 'g-');
grid on;
hold on
plot(Qr, 'r-');
hold on
plot(pktfu, 'r-');
hold on

figure;
plot(pkt, 'g-');
grid on;
hold on

clc;
% Переданные данные
%pkt

clear all;
Br = 400e3;
Fd = Br*8;

Fdopl = 100;

D = (rand(1, 10000) > 0.5) * 2 - 1;
% redescr
Ndkr = Fd/Br;
pkt = D(ceil(0.01:1/Ndkr:length(D)));
l = length(pkt);

s = sin(2*pi*Fdopl*(0:1/Fd:l/Fd));
c = cos(2*pi*Fdopl*(0:1/Fd:l/Fd));

I = s(1:l) .* pkt;
Q = c(1:l) .* pkt;

figure; hold on; plot(I); plot(Q, 'r');
figure; hold on; plot(I, Q, '*');

%% Петля Костаса на базе VCO

%% Модуляция
clear;

% Битовая картика профиля pacman
image = [ 1, 1, 1, 0, 0, 1, 1, 1;
1, 1, 0, 0, 0, 0, 1, 1;
1, 0, 0, 0, 1, 0, 0, 1;
0, 0, 0, 0, 0, 0, 1, 1;
0, 0, 0, 0, 0, 1, 1, 1;
1, 0, 0, 0, 0, 0, 0, 1;
1, 1, 0, 0, 0, 0, 1, 1;
1, 1, 1, 0, 0, 1, 1, 1 ];

% Начальный символ приема
sync_symbol = [ 1 ];
% Битовый поток с начальным символом
bit_stream = [ sync_symbol ];

% Генерация битового потока из матрицы профиля pacman
for j = 1:8
bit_stream = [bit_stream image(j,];
end

%bit_stream = (rand(1, 6) > 0.5);
%bit_stream = [1 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0];
bit_stream = [1 1 0 1 0];
N = length(bit_stream);
fprintf('bit_stream %d\n', bit_stream);

% Частота семплирования 288 кГц
fs = 288000;
% Несущая частота относительно частоты выборок 36 кГц
f = 0.125;
% Частота доплера 300 Гц
fdop = f/1200;
% Символьная частота 100 Гц
fsim = 100;
% Количество выборок на период несущей 8
ns = 1/f;
% Количество выборок необх. для одного символа (360*8 периодов нес. f)
% Символьная скорость 100 Гц
tsymbol = fs*f*ns/fsim;
fprintf('tsymbol %d\n', tsymbol);
% Счетчик выборок в 0
t = 0;
% Количество символов переданных с нулевой фазой
nDummySymbols = 4;
% Общее количество выборок
%ttotal = N*tsymbol + nDummySymbols * tsymbol + tsymbol*3;
ttotal = N*tsymbol + nDummySymbols * tsymbol;
fprintf('ttotal %d\n', ttotal);
% Создать пустой сигнал
BPSK_signal = zeros(ttotal,1);

% Символы с нулевой фазой
for i = 1:nDummySymbols
for j = 1:tsymbol
BPSK_signal(t+1) = cos( 2*pi*(f+fdop)*t );
BPSK_signal_nomod(t+1) = cos(2*pi*(f+fdop)*t);
t = t + 1;
end
end

fprintf('t %d\n', t);

% BPSK сигнал
for i = 1:length(bit_stream)
phase = pi * bit_stream(i);
for j = 1:tsymbol
BPSK_signal(t+1) = cos(2*pi*(f+fdop)*t + phase);
BPSK_signal_mod(t+1) = cos(2*pi*(f+fdop)*t + phase);
BPSK_signal_nomod(t+1) = 0;
t = t + 1;
end
end

fprintf('length BPSK_signal %d\n', length(BPSK_signal));
fprintf('t %d\n', t);

%% AWGN (аддитивный белый гауссов шум) канал передачи данных
% отношение сигнал/шум в децибелах
SNR = 20;
% мощность сигнала в децибелах
SIGPOWER = 10;
% смесь мод.сигнала с шумом
BPSK_signal = awgn( BPSK_signal, SNR, SIGPOWER);

%% Демодуляция

% полученная картинка в 0
rec_image = zeros(8,8);

% Координаты картинки
x = 1;
y = 1;

% начальный шаг в 0
t = 0;

% Все матрицы в 0
carrier_inph = zeros(ttotal,1);
mixer_inph = zeros(ttotal,1);
demod_thr_inph = zeros(ttotal,1);

% Фильтр демодулятора
b = fir1(100, 0.1);
% Память состояния демодулятора
stateLowpassInph = zeros(numel(-1,1);

% переменные для PLL осциллятора
% cosine output
c = 1;
% delayed cosine by one timestep
c_delay = 0;
% sine output
s = 0;
% delayed sine by one timestep
s_delay = 0;

% Данные из VCO для целей отладки
sine = zeros(ttotal,1);
cosine =zeros(ttotal,1);
vco = zeros(ttotal,1);

% PLL loop фильтр
bPLL = fir1(40, 0.1);
% and the state memory
zfPLLa = zeros(numel(bPLL)-1,1);
zfPLLc = zeros(numel(bPLL)-1,1);
zfPLLs = zeros(numel(bPLL)-1,1);

% Обратный отсчет. При нулевом значении отбирается символ.
% Начальный обратный отсчет учитывает задержку
% fir lowpass и пропускает начальный символ
sampleMoment = tsymbol/2+tsymbol;

% Это время необходимое PLL для стабилизации на несущей
nPLLsettle = (nDummySymbols/2) * tsymbol;

% Игнорировать данные, пока не получен символ "1"
ignoreData = 1;

% VCO шаг
VCOgain = 0.005;

% VCO, которое управляет частотой. при v = 0 он точно равен f.
v = 0;

% Основной цикл
for step = 1:ttotal

% Получить значение BPSK_signal
sample = BPSK_signal(step);

% PLL - "voltage" управляемый осциллятор с +/- v входом
% w0 частота VCO
w0 = 2*pi*(f+v*VCOgain);
% Cохранить предыдущее состояние
c_delay = c;
s_delay = s;
% Подсчет нового состояния
c = c_delay * cos(w0) - s_delay * sin(w0);
s = s_delay * cos(w0) + c_delay * sin(w0);
% Cохраним все в удобных векторах для графики
sine(step) = s;
cosine(step) = c;
vco(step) = v;
% конец VCO

% Сохраняем sine/cosine. Обратите внимание, что инфазная волна
% теперь является синусоидальной волной, потому что PLL с мультипликацией
% в качестве фазового детектора имеет сдвиг фазы на 90 градусов между
% входом и выходом
carrier_inph(step) = c;

% Это решает, находимся ли мы в режиме синхронизации PLL или режима замораживания
% Примечание: это должно быть заменено блокировкой PLL в реальном приеме,
% потому что мы не знаем, когда передача фактически
% начинается. Здесь мы просто предполагаем, что она начинается с самого начала

% PLL in action: изменение частоты через v
% Фазовый детектор
pc = sample*c;
ps = sample*-s;

% filtering out the 2f term
% Фильтрация pc
[vc,zfPLLc] = filter(bPLL,1,pc,zfPLLc);

% filtering out the 2f term
% Фильтрация ps
[vs,zfPLLs] = filter(bPLL,1,ps,zfPLLs);

[v,zfPLLa] = filter(bPLL,1,vc*vs,zfPLLa);

% we sneak here a -1 in because this is ambigous
% шаг на -1
mixer_inph(step) = -sample*carrier_inph(step);

[lp_demod_inph(step),stateLowpassInph] = filter(b,1,mixer_inph(step),stateLowpassInph);

demod_thr_inph(step) = lp_demod_inph(step);

if (ignoreData)
% Ждем стартового символа 1
if ( lp_demod_inph(step) > 0.25 )
% got it!
% Стартовый символ прошел
ignoreData = 0;
end
else
% Ожидание середины символа для выборки демодулятора
sampleMoment = sampleMoment - 1;
if (sampleMoment == 0)
sampleMoment = tsymbol;
if ( demod_thr_inph(step) > 0.2 )
lsb = 1;
else
lsb = 0;
end

fprintf('lsb %d\n', lsb);

% debug samples: add small spikes
% отладка: добавка небольших пиков
lp_demod_inph(step) = lp_demod_inph(step) + 0.1;

% Сохранить результат в картинку
if (((x<9) && (y<9)) && (ignoreData == 0))
rec_image(y,x) = lsb;
x = x + 1;
if (x > 8)
x = 1;
y = y + 1;
end % if (x > 8)
end % x>9 && y>9
end % sampleMoment == 0
end % ignore data
% Следующий шаг в выборке
t = t + 1;
end

%% Графика

% Немодулированный сигнал
figure;
plot(BPSK_signal_nomod, 'b-', 'LineWidth', 2);
hold on;
grid on;
title('Немодулированный сигнал');

% Модулированный сигнал
figure;
plot(BPSK_signal_mod, 'b-', 'LineWidth', 2);
hold on;
grid on;
title('Модулированный сигнал');

% Полный BPSK_signal
figure;
plot(BPSK_signal, 'b-', 'LineWidth', 2);
hold on;
grid on;
title('BPSK сигнал');

% Демодулированный сигнал
figure;
subplot(2,1,1);
plot(lp_demod_inph,'LineWidth',2);
xlabel('Samples');
ylabel('Amplitude');
title('demod / lowpass filtered');
grid on;

% Картинка
subplot(2,1,2);
imshow(rec_image);

% Напряжение VCO
figure;
subplot(1,1,1);
plot(vco);

%% Петля Костаса на базе VCO

%% Модуляция
clear;

% Битовая картика профиля pacman
% image = [ 1, 1, 1, 0, 0, 1, 1, 1;
% 1, 1, 0, 0, 0, 0, 1, 1;
% 1, 0, 0, 0, 1, 0, 0, 1;
% 0, 0, 0, 0, 0, 0, 1, 1;
% 0, 0, 0, 0, 0, 1, 1, 1;
% 1, 0, 0, 0, 0, 0, 0, 1;
% 1, 1, 0, 0, 0, 0, 1, 1;
% 1, 1, 1, 0, 0, 1, 1, 1 ];

image = [ 1, 0, 1, 0, 1, 0, 1, 0;
0, 1, 0, 1, 0, 1, 0, 1;
1, 0, 1, 0, 1, 0, 1, 0;
0, 1, 0, 1, 0, 1, 0, 1;
1, 1, 1, 1, 1, 1, 1, 1;
0, 0, 0, 0, 0, 0, 0, 0;
1, 1, 1, 1, 1, 1, 1, 1;
0, 0, 0, 0, 0, 0, 0, 0 ];

% Начальный символ приема
sync_symbol = [ 1 ];
% Битовый поток с начальным символом
bit_stream = [ sync_symbol ];

% Генерация битового потока из матрицы профиля pacman
for j = 1:8
bit_stream = [bit_stream image(j,];
end

%bit_stream = (rand(1, 20) > 0.5);
%bit_stream = [1 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 1 0 1 0 1 0 1 0];
%bit_stream = [1 1 0 1 0 1 0 1 0];
N = length(bit_stream);
fprintf('bit_stream %d\n', bit_stream);

% Частота семплирования 288 кГц (8 выборок на период)
fs = 288000;
% Несущая частота относительно частоты выборок 36 кГц
f = 0.125;
% Частота доплера 60 Гц (~2.5 м.сек)
fdop = f/600;

% % Частота семплирования 144 кГц (4 выборки на период)
% fs = 144000;
% % Несущая частота относительно частоты выборок 36 кГц
% f = 0.25;
% % Частота доплера 60 Гц (~2.5 м/сек)
% fdop = -f/600;

% Частота среза НЧ фильтра в демодуляторе
fc = 100;
% Частота Найквиста
fn = fs/2;
% Порядок НЧ фильтра демодулятора
lowpfiltord = 50;
% Порядок loop фильтра демодулятора
loopfiltord = 50;

% Амплитуда
A = 1.0;
% Символьная скорость 100 Гц
fsim = 100;
% Количество выборок на период несущей
ns = 1/f;
% Количество выборок необх. для одного символа (360*8 периодов нес. f)
% при символьной скорости fsim
tsymbol = fs*f*ns/fsim;
fprintf('tsymbol %d\n', tsymbol);
% Счетчик выборок в 0
t = 0;
% Количество символов переданных с нулевой фазой (пилот)
nDummySymbols = 4;
% Общее количество выборок
%ttotal = N*tsymbol + nDummySymbols * tsymbol + tsymbol*3;
ttotal = N*tsymbol + nDummySymbols * tsymbol;
fprintf('ttotal %d\n', ttotal);
% Создать пустой BPSK сигнал
BPSK_signal_t = zeros(ttotal,1);

% Пилот сигнал
for i = 1:nDummySymbols
for j = 1:tsymbol
BPSK_signal_t(t+1) = A*sin( 2*pi*(f+fdop)*t );
BPSK_signal_nomod(t+1) = A*sin(2*pi*(f+fdop)*t);
t = t + 1;
end
end

fprintf('t %d\n', t);

% BPSK сигнал
for i = 1:length(bit_stream)
phase = pi * bit_stream(i);
for j = 1:tsymbol
BPSK_signal_t(t+1) = A*sin(2*pi*(f+fdop)*t + phase);
BPSK_signal_mod(t+1) = A*sin(2*pi*(f+fdop)*t + phase);
BPSK_signal_nomod(t+1) = 0;
t = t + 1;
end
end

fprintf('length BPSK_signal %d\n', length(BPSK_signal_t));
fprintf('t %d\n', t);

%% AWGN (аддитивный белый гауссов шум) + полосовая фильтрация

% отношение сигнал/шум в децибелах
SNR = 9;
% мощность сигнала в децибелах
SIGPOWER = 12;
% смесь мод.сигнала с шумом
BPSK_signal_awgn = awgn( BPSK_signal_t, SNR, SIGPOWER );

% Порядок полосового фильтра
bandpfiltord = 50;
% нижняя частота полосового фильтра 35КГц
W1=35e3/(0.5*fs);
% верхняя частота полосового фильтра 37КГц
W2=37e3/(0.5*fs);
% расчет окна Ханна
k=hann(bandpfiltord+1);
% расчет коэффициентов входного полосового фильтра
bp=fir1(bandpfiltord, [W1 W2], k);
% Память состояния фильтра
stateBandpassBpsk = zeros(numel(bp)-1,1);
% Полосовая фильтрация BPSK сигнала + AWGN шум
[BPSK_signal_awgn_filtr,stateBandpassBpsk] = filter(bp,1,BPSK_signal_awgn,stateBandpassBpsk);
% BPSK_signal_awgn_filtr = filter(bp, 1, BPSK_signal_awgn);

% BPSK_signal
%BPSK_signal = BPSK_signal_awgn;
BPSK_signal = BPSK_signal_awgn_filtr;


%% Демодуляция

% полученная картинка в 0
rec_image = zeros(8,8);

% Координаты картинки
x = 1;
y = 1;

% начальный шаг в 0
t = 0;

% Все матрицы в 0
carrier_inph = zeros(ttotal,1);
mixer_inph = zeros(ttotal,1);
demod_thr_inph = zeros(ttotal,1);

%расчет окна Кайзера
kb=kaiser(lowpfiltord+1);
% Коэффициенты НЧ фильтра демодулятора порядок lowpfiltord
b = fir1(lowpfiltord, fc/fn, kb);
% Память состояния для НЧ фильтра демодул.
stateLowpassInph = zeros(numel(-1,1);

% переменные для PLL осциллятора
% cosine output
c = 1;
% delayed cosine by one timestep
c_delay = 0;
% sine output
s = 0;
% delayed sine by one timestep
s_delay = 0;

% Данные из VCO для целей отладки
sine = zeros(ttotal,1);
cosine =zeros(ttotal,1);
vco = zeros(ttotal,1);

%расчет окна Ханна
kp=hann(loopfiltord+1);
% Коэффициенты PLL loop фильтра
bPLL = fir1(loopfiltord, fc/fn, kp);
% Обнулить память
zfPLLa = zeros(numel(bPLL)-1,1);
zfPLLc = zeros(numel(bPLL)-1,1);
zfPLLs = zeros(numel(bPLL)-1,1);

% Обратный отсчет. При нулевом значении отбирается символ.
% Начальный обратный отсчет учитывает задержку
% fir lowpass и пропускает начальный символ
sampleMoment = tsymbol/2+tsymbol;

% Это время необходимое PLL для стабилизации на несущей
nPLLsettle = (nDummySymbols/2) * tsymbol;

% Игнорировать данные, пока не получен символ "1" (синхробит)
ignoreData = 1;

% VCO шаг
VCOgain = 0.005;

% VCO, которое управляет частотой. при v = 0 он точно равен f.
v = 0;

% Основной цикл
for step = 1:ttotal

% Полосовая фильтрация BPSK сигнала + AWGN шум
% [BPSK_signal(step),stateBandpassBpsk] = filter(bp,1,BPSK_signal_awgn(step),stateBandpassBpsk);

% Получить значение BPSK_signal в sample
sample = BPSK_signal(step);

% PLL - управляемый осциллятор с +/- v входом
% w0 частота VCO
w0 = 2*pi*(f+v*VCOgain);
% Cохранить предыдущее состояние
c_delay = c;
s_delay = s;
% Подсчет нового состояния
c = c_delay * cos(w0) - s_delay * sin(w0);
s = s_delay * cos(w0) + c_delay * sin(w0);
% Cохраним все в удобных векторах для графики
sine(step) = s;
cosine(step) = c;
vco(step) = v;
% конец VCO

% Сохраняем sine/cosine. Обратите внимание, что инфазная волна
% теперь является синусоидальной волной, потому что PLL с мультипликацией
% в качестве фазового детектора имеет сдвиг фазы на 90 градусов между
% входом и выходом
carrier_inph(step) = c;

% Это решает, находимся ли мы в режиме синхронизации PLL или режиме замораживания
% Примечание: это должно быть заменено блокировкой PLL в реальном приеме,
% потому что мы не знаем, когда передача фактически
% начинается. Здесь мы просто предполагаем, что она начинается с самого начала

% PLL изменение частоты через v
% Фазовый детектор
pc = sample*c;
ps = sample*-s;

% Фильтрация pc
[vc,zfPLLc] = filter(bPLL,1,pc,zfPLLc);

% Фильтрация ps
[vs,zfPLLs] = filter(bPLL,1,ps,zfPLLs);

[v,zfPLLa] = filter(bPLL,1,vc*vs,zfPLLa);

% шаг на -1
mixer_inph(step) = sample*carrier_inph(step);

[lp_demod_inph(step),stateLowpassInph] = filter(b,1,mixer_inph(step),stateLowpassInph);

demod_thr_inph(step) = lp_demod_inph(step);

if (ignoreData)
% Ждем стартового символа 1
if ( lp_demod_inph(step) > 0.25 )
% Стартовый символ прошел
ignoreData = 0;
end
else
% Ожидание середины символа для выборки демодулятора
sampleMoment = sampleMoment - 1;
if (sampleMoment == 0)
sampleMoment = tsymbol;
if ( demod_thr_inph(step) > 0.1 )
lsb = 1;
else
lsb = 0;
end

fprintf('lsb %d\n', lsb);

% отладка: добавка небольших пиков
lp_demod_inph(step) = lp_demod_inph(step) + 0.1;

% Сохранить результат в картинку 8:8
if (((x<9) && (y<9)) && (ignoreData == 0))
rec_image(y,x) = lsb;
x = x + 1;
if (x > 8)
x = 1;
y = y + 1;
end % if (x > 8)
end % x>9 && y>9
end % sampleMoment == 0
end % ignoreData
% Следующий шаг в выборке
t = t + 1;
end

%% Графика

% % Немодулированный сигнал
% figure;
% subplot(2,1,1);
% plot(BPSK_signal_nomod, 'b-', 'LineWidth', 2);
% hold on;
% grid on;
% title('Немодулированный сигнал (пилот)');
%
% % Модулированный сигнал
% %figure;
% subplot(2,1,2);
% plot(BPSK_signal_mod, 'r-', 'LineWidth', 2);
% hold on;
% grid on;
% title('Модулированный сигнал (данные)');

% Полный BPSK_signal
figure;
subplot(2,1,1);
plot(BPSK_signal_awgn, 'b-', 'LineWidth', 2);
xlabel('Выборки');
hold on;
grid on;
title('Полный BPSK сигнал + шум AWGN');

subplot(2,1,2);
plot(BPSK_signal, 'r-', 'LineWidth', 2);
xlabel('Выборки');
hold on;
grid on;
title('Полный BPSK сигнал + шум AWGN + Фильтр');

% Демодулированный сигнал
figure;
subplot(2,1,1);
plot(lp_demod_inph,'LineWidth',2);
xlabel('Выборки');
ylabel('Амплитуда');
title('Демодулированный сигнал / НЧ фильтрация');
grid on;

% Картинка
subplot(2,1,2);
imshow(rec_image);
title('Битовая картинка 8:8 квадратиков 1-белый, 0-черный кв. ' );

% Напряжение VCO
% figure
% plot(vco);
% title('Напряжение VCO');
%
% % Фильтр демодулятора
% figure
% % АЧХ НЧ фильтра
% [h,f]=freqz(b,1);
% % Параметры являются соответственно частотой и амплитудой
% plot(f*fs/(2*pi),20*log10(abs(h)))
% xlabel('частота/Hz');ylabel('амплитуда/dB');title('АЧХ НЧ фильтра демодулятора');
%
% % Фильтр полосовой на 36 кГц
% figure
% [h,fp]=freqz(bp,1);%amplitude-frequency characteristic graph
% plot(fp*fs/(2*pi),20*log10(abs(h)))%parameters are respectively frequency and amplitude
% xlabel('частота/Hz');ylabel('амплитуда/dB');title('АЧХ полосового фильтра');