function fsk4_states = fsk4_init(fsk4_states,fsk4_info)

       Fs = 48000; fsk4_states.Fs = Fs;   %частота дискритизации

       Rs = fsk4_info.rs; fsk4_states.Rs = Rs;    %скорость передачи символов

   M = fsk4_states.Fs/fsk4_states.Rs; fsk4_states.M = M; %бит на символ

   

    % Настройка 4FSK приподнятый косинус.

 

    empty_filter = [zeros(1,99) 1];

 

    rf = (0:(Fs/2));

    if fsk4_info.no_filter

      fsk4_states.tx_filter = empty_filter;

      fsk4_states.rx_filter = empty_filter;

    else

      fsk4_states.tx_filter = fir2(400,rf/(Fs/2),fsk4_info.tx_filt_resp(rf));

      fsk4_states.rx_filter=fir2 (400, rf/(Fs/2),fsk4_info.rx_filt_resp(rf));

    end

 

    %fsk4_states.tx_filter = fsk4_states.rx_filter = [zeros(1,99) 1];

    %Настройка значений 4FSK

    fsk4_states.symmap = fsk4_info.syms / fsk4_info.max_dev;

   

    fm_states.Ts = M;

    fm_states.Fs = Fs;

    fm_states.fc = 0;

    fm_states.fm_max = fsk4_info.max_dev*2;

    fm_states.fd = fsk4_info.max_dev;

    fm_states.pre_emp  = 0; %= fm_states.de_emp

    fm_states.output_filter = 0;

    fm_states.ph_dont_limit = 1;

    fsk4_states.fm_states = analog_fm_init(fm_states);

    fsk4_states.modinfo = fsk4_info;

    fsk4_states.verbose = 0;

end

 

В функции fsk4_init происходит объявление параметров. Так мы видим, что частота дискретизации равна 48000 Гц, а также объявляется скорость передачи символов и значение бит, приходящихся на символ. Затем настраивается фильтр приподнятого косинуса и вводятся значения 4FSK модуляции такие как период дискретизации, частота дискретизации, частота среза, девиация, максимальное достигаемое значение частоты.

function d = idmp(data, M)

    d = zeros(1,length(data)/M);

    for i = 1:length(d)

      d(i) = sum(data(1+(i-1)*M:i*M));

    end

end

 

В функции idmp настраивается интеграция по каналу передачи данных и  обнуления каждые М выборок. Создается накопленная сумма входного сигнала дискретного по времени, в то время как нужно сбросить сумму до нуля в соответствии с установленным графиком. Когда начинается моделирование, блок отбрасывается количество образцов, указанных в параметрах. После этого начального периода, блок суммирует входной сигнал вдоль столбцов и сбрасывает сумму в ноль каждые M входных выборок.

% DMR модулятор ----------------------------------------------------------

 

function [tx, tx_filt, tx_stream] = fsk4_mod(fsk4_states, tx_bits)

  verbose = fsk4_states.verbose;

 

  hbits = tx_bits(1:2:length(tx_bits));

  lbits = tx_bits(2:2:length(tx_bits));

  %Последовательность нечётной длины бит

  if(length(hbits)~=length(lbits))

    lbits = [lbits 0];

  end

  tx_symbols = lbits + hbits*2 + 1;

  M = fsk4_states.M;

  nsym = length(tx_symbols);

  nsam = nsym*M;

 

  tx_stream = zeros(1,nsam);

  for i=1:nsym

    tx_stream(1+(i-1)*M:i*M) = fsk4_states.symmap(tx_symbols(i));

  end

  tx_filt = filter(fsk4_states.tx_filter, 1, tx_stream);

  tx = analog_fm_mod(fsk4_states.fm_states, tx_filt);

 

  if verbose

    figure(10);

    plot(20*log10(abs(fft(tx))))

    title('Spectrum of modulated 4FSK')

  end

 

end

 

Функция fsk4_mod непосредственно говорит о процессе модуляции DMR. Задаются битовые последовательности и накладываются определённые условия. В завершение функции строится спектр на выходе 4FSK модулятора.

% Интецрация в 4FSK демодулятора ----------------------------------------------------

 

function bits = fsk4_demod_thing(fsk4_states, rx)

 

  M = fsk4_states.M;

  Fs = fsk4_states.Fs;

  verbose = fsk4_states.verbose;

  t = (0:length(rx)-1);

  symup = fsk4_states.modinfo.syms;

 

  % Интегратор подобен КИХ-фильтру с прямоугольным окном.

 

  Fs = fsk4_states.Fs;

  rf = (0:(Fs/2));

  rx_filter_a = fir1(100 ,.2);

  rx_filter_b = fsk4_states.rx_filter;

  rx_filter_n = [zeros(1,99) 1];

 

  rx = filter(rx_filter_b, 1, rx);

 

  sym1m = exp(-1i*2*pi*(symup(1)/Fs)*t).*rx;

  sym2m = exp(-1i*2*pi*(symup(2)/Fs)*t).*rx;

  sym3m = exp(-1i*2*pi*(symup(3)/Fs)*t).*rx;

  sym4m = exp(-1i*2*pi*(symup(4)/Fs)*t).*rx;

 

  %Будет ли меняться импульсная хар-ка фильтра в зависимости от изменения задержки.

 

  % добавлении М к нему грубо синхронизация будет регулировать на 1.

 

  fine_timing = 54;   % точная синхронизация

 

  sym1m = idmp(sym1m(fine_timing:length(sym1m)),M);

sym1m =(real(sym1m).^2+imag(sym1m).^2);

  sym2m = idmp(sym2m(fine_timing:length(sym2m)),M);

sym2m = (real(sym2m).^2+imag(sym2m).^2);

  sym3m = idmp(sym3m(fine_timing:length(sym3m)),M);

sym3m = (real(sym3m).^2+imag(sym3m).^2);

  sym4m = idmp(sym4m(fine_timing:length(sym4m)),M);

sym4m = (real(sym4m).^2+imag(sym4m).^2);

 

 

  figure(2);

  nsym = 500;

  %subplot(411); plot(sym1m(1:nsym))

  %subplot(412); plot(sym2m(1:nsym))

  %subplot(413); plot(sym3m(1:nsym))

  %subplot(414); plot(sym4m(1:nsym))

  plot((1:nsym),sym1m(1:nsym),(1:nsym),sym2m(1:nsym),(1:nsym),sym3m(1:nsym),(1:nsym),sym4m(1:nsym))

 

  [x iv] = max ([sym1m; sym2m; sym3m; sym4m;]);

  bits = zeros(1,length(iv*2));

  figure(3);

  hist(iv);

  for i=1:length(iv)

    bits(1+(i-1)*2:i*2) = [[0 0];[0 1];[1 0];[1 1]]*(iv((i),(1:2)));

  end

end

 

Функция fsk4_demod_thing поясняет процесс интеграции в 4FSK демодуляторе. Считается, что интегратор подобен КИХ-фильтру с прямоугольным окном. Произведено исследование о зависимости и изменении импульсной характеристики фильтра с изменением задержки, а также построены графики этих зависимостей. Благодаря этой функции в сети осуществляется синхронизация.

 

 

function dat = bitreps(in,M)

  dat = zeros(1,length(in)*M);

  for i=1:length(in)

    dat(1+(i-1)*M:i*M) = in(i);

  end

end

 

function syms = mrd4(bits)

  syms = zeros(1,length(bits));

  rd=0;

  lastsym=0;

  for n = (1:length(bits))

    bit = bits(n);

    sp = [1 3]*(bit+1); %уровень бита +1 или +3

    [x,v] = min(abs([rd+sp rd-sp]));

    ssel = [sp -sp]*(v);

    if(ssel == lastsym)ssel = -ssel;

end %никогда 2 одинаковых символа подряд

    syms(n) = ssel; %появляется символ

    rd = rd + ssel; %обновление запуска неравенста

    lastsym = ssel; %запоминание этого символа для след промежутка времени

  end

  end

end

 

В данной функции mrd4 продемонстрирован минимальный запуск 4-х символьного кодера. Это простой кодер, где 1 бит приходится на 1 символ кодировки для 4FSK модемов, построенных на старых FM модуляторах. Устанавливается, что уровни бита +1 или +3. Из условия если следующий символ не равен предыдущему, получается новый символ, затем обновляется запуск неравенства, и новый символ запоминается для следующего промежутка времени.

function syms = mrd8(bits)

  bitlen = length(bits);

  if mod(bitlen,2) == 1

    bits = [bits 0];

  end

 

  syms = zeros(1,length(bits)*.5);

  rd=0;

  lastsym=0;

  for n = (1:2:length(bits))

    bit = (bits(n)*2)+bits(n+1);

    sp = [1 3 7 5]*(bit+1); %уровень бита +1 или +3

    [x,v] = min(abs([rd+sp rd-sp])); %Выберите + n или -n , в зависимости от того сводит к минимуму неравенство

    ssel = [sp -sp]*(v);

    if(ssel == lastsym)ssel = -ssel;

    end %никогда 2 одинаковых символа подряд

    syms((n+1)/2) = ssel; %испускают символ

    rd = rd + ssel; %обновление запуска неравенста

    lastsym = ssel; %запоминание этого символа для след промежутка времени

  end

  end

 

Функция mrd8 подобна предыдущей, однако здесь и используется 8-ми символьный кодер. В данном кодере на один символ приходится 2 бита. Далее алгоритм действий тот же, однако вместо 2-х возможных значений уже имеется 4.

% кодирование "Манчестер 4"

function syms = mane4(bits)

    syms = zeros(1,floor(bits/2)*2);

    for n = (1:2:length(bits))

      bit0 = bits(n);

      bit1 = bits(n+1);

      sel = 2*bit0+bit1+1;

      syms(n:n+1) = [[3 -3];[-3 3];[1 -1];[-1 1]]*( sel(1:2));

    end

end

 

 Данная функция mane4 осуществляет алгоритм кодирования типа «Манчестер-4». При манчестерском кодировании каждый такт делится на две части. Информация кодируется перепадами потенциала в середине каждого такта. Единица кодируется перепадом от низкого уровня сигнала к высокому, а ноль — обратным перепадом. Так как сигнал изменяется, по крайней мере, один раз за такт передачи одного бита данных, то манчестерский код обладает хорошими самосинхронизирующими свойствами. Важным достоинством манчестерского кода является отсутствие необходимости в синхронизующих битах (как в NRZ-коде) и, вследствие этого, данные могут передаваться подряд сколь угодно долго, из-за чего плотность данных в общем потоке кода приближается к 100 %.

 

 

function [bits err ,xphi] = fsk4_demod_fmrid(fsk4_states, rx, enable_fine_timing) % = 0

  %Демодуляция fsk сигнала при помощи аналогового fm модулятора

  rxd = analog_fm_demod(fsk4_states.fm_states,rx);

 

  M = fsk4_states.M;

  verbose = fsk4_states.verbose;

  fine_timing = 61;

 

  %для тонкой оценки синхронизации

  %идеальное время

  %fine_timing = 59;

 

  %Фильтр приподнятого косинуса для избавления от некоторого шума

  rxd = filter(fsk4_states.rx_filter, 1, rxd);

 

  %Попробуйте выяснить, где выборка должна произойти более 30 периодов символа

  diffsel = fold_sum(abs(diff( rxd(3001:3001+(M*30)) )),10);

 

  if verbose

    figure(11);

    plot(diffsel);

    title('Fine timing estimation');

  end

 

  %регулировки точной синхронизации

  [v iv] = min(diffsel);

  if enable_fine_timing

    fine_timing = 59 + iv;

  end

  rxphi = iv;

 

  %простые символы

  sym = rxd(fine_timing:M:length(rxd));

 

  if verbose

    figure(4)

    plot(sym(1:1000));

    title('Sampled symbols')

  end

  %глазковая диаграмма(afsym,2);

  % После отбора проб , то дальние символы имеют тенденцию распространяться о .80

 

  % обратное отображение символов

  % Возьмите гистограмма выборки символов , найти центр наибольшего распределения, а также исправить символ карты , чтобы соответствовать его

  [a b] = hist(abs(sym),50);

  [a ii] = max(a);

  %grmax = abs(b(ii));

  %grmax = (grmax<.65)*.65 + (grmax>=.65)*grmax;

  grmax = .84;

  dmsyms = rot90(fsk4_states.symmap*grmax);

  (dmsyms(2)+dmsyms(1))/2;

 

  if verbose

    figure(2)

    hist(abs(sym),200);

    title('Sampled symbol histogram')

  end

 

  %символы выполняют обратное отображение

  [err, symout] = min(abs(sym-dmsyms));

 

  if verbose

    figure(3)

    hist(symout);

    title('De-mapped symbols')

  end

 

  bits = zeros(1,length(symout)*2);

  %Перевод символов обратно в биты

 

  for i=1:length(symout)

    bits(1+(i-1)*2:i*2) = [[1 1];[1 0];[0 1];[0 0]] *(symout((i),(1:2)));

  end

  end

 

 

fsk4_demod_fmrid – функция выполняет демодуляцию FSK-модулированного сигнала при помощи просто аналогового FM демодулятора. Для начала выбирается время синхронизации. Для тонкой оценки времени выбирается 59. После FM демодулятора используется фильтр приподнятого косинуса для избавления от некоторого шума. Определяется где произойдёт выборка после 30 периодов символа. Строится график оценки точной синхронизации. Затем регулировка точной синхронизации. И в результате на графике можно видеть простые символы. Построение глазковой диаграммы, по гистограмме находится центр наибольшего распределения, символы выполняют обратное отображение, что наблюдается на графике и переводятся обратно в биты.

 

Частотные характеристики фильтра типа «приподнятый косинус»

% из стандарта ETSI TS 102 361-1 V2.2.1 page 111

dmr.tx_filt_resp = @(f) sqrt(1.0*(f<=1920) - cos((pi*f)/1920).*1.0.*(f>1920 & f<=2880));

dmr.rx_filt_resp = dmr.tx_filt_resp;

dmr.max_dev = 1944;

dmr.syms = [-1944 -648 1944 648];

dmr.rs = 4800;

dmr.no_filter = 0;

dmr.demod_fx = @fsk4_demod_fmrid;

global dmr_info = dmr;

 

 

% для неидеального 4fsk фильтра

nfl.tx_filt_resp = @(f) 1;

nfl.rx_filt_resp = nfl.tx_filt_resp;

nfl.max_dev = 7200;

%nfl.syms = [-3600 -1200 1200 3600];

nfl.syms = [-7200,-2400,2400,7200];

nfl.rs = 4800;

nfl.no_filter = 1;

nfl.demod_fx = @fsk4_demod_thing;

global nflt_info = nfl;

 

На данном участке кода задаются параметры фильтров.

% Тестовая вероятность битовой ошибки ----------------------------------------------------------

% Параметры - aEsNodB - EbNo в децибелах

%        - timing_offset - как далеко смещена синхронизация

%        - bitcnt - сколько бит для проверки

%        - demod_fx - функция демодулятора

% Returns - ber - teh measured BER

%         - thrcoh - теория BER когерентного демодулятор

%         - thrncoh - теория BER некогерентного демодулятор

function [ber thrcoh ,hrncoh] = nfbert(aEsNodB,modem_config, bitcnt, timing_offset)% bitcnt=100000 timing_offset = 10

 

  rand('state',1);

  randn('state',1);

 

  %Сколько битов этот тест должен работать?

  bitcnt = 120000;

 

  test_bits = [zeros(1,100) rand(1,bitcnt)>.5]; %Random bits. Pad with zeros to prime the filters

  fsk4_states.M = 1;

  fsk4_states = fsk4_init(fsk4_states,modem_config);

 

  %Установливается значение 0 , для сокращения записи

  fsk4_states.verbose = 1;

  Fs = fsk4_states.Fs;

  Rb = fsk4_states.Rs * 2;  %Умножаем скорость передачи символов на 2, так как мы имеем 2 бита на символ

 

  tx = fsk4_mod(fsk4_states,test_bits);

 

  %добавить шум

  EsNo = 10^(aEsNodB/10);

  EbNo = EsNo;

  variance = Fs/(Rb*EbNo);

  nsam = length(tx);

  noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam));

  rx    = tx*exp(j*pi/2) + noise;

 

  rx    = rx(timing_offset:length(rx));

 

  rx_bits = modem_config.demod_fx(fsk4_states,rx);

  ber = 1;

 

  %нужно учитывать смещение от входных данных для вывода данных

  %Нет обнаружения преамбулы

  ox = 1;

  for offset = (1:100)

    nerr = sum(xor(rx_bits(offset:length(rx_bits)),test_bits(1:length(rx_bits)+1-offset)));

    bern = nerr/(bitcnt-offset);

    if(bern < ber)

      ox = offset;

      best_nerr = nerr;

    end

    ber = min([ber bern]);

  end

  offset = ox;

  disp(ncoarse, timing); %d nerr: %d\n", offset, best_nerr);

 

  % Coherent BER theory

  thrcoh = erfc(sqrt(EbNo));

 

  % Расчет теории BER некогерентного

 

  ms = 4;

  ns = (1:ms-1);

  as = (-1).^(ns+1);

  bs = (as./(ns+1));

 

  cs = ((ms-1)./ns);

 

  ds = ns.*log2(ms);

  es = ns+1;

  fs = exp( -(ds./es)*EbNo );

 

  thrncoh = ((ms/2)/(ms-1)) * sum(bs.*((ms-1)./ns).*exp( -(ds./es)*EbNo ));

 

end

 

Функция nfbert указывает сколько бит отводится на проверку, в данном случае тестовый сигнал содержит 120000. В данной функции происходит расчёт зависимости вероятности ошибки от ОСШ.

 

function rxphi = fine_ex(timing_offset) %временное смещение = 1

  global dmr_info;

 

  global nflt_info;

 

  rand('state',1);

  randn('state',1);

 

  bitcnt = 12051;

  test_bits = [zeros(1,100) rand(1,bitcnt)>.5]; %Случайные биты .

  t_vec = [0 0 1 1];

  %test_bits = repmat(t_vec,1,ceil(24000/length(t_vec)));

 

 

  fsk4_states.M = 1;

  fsk4_states = fsk4_init(fsk4_states,dmr_info);

  Fs = fsk4_states.Fs;

  Rb = fsk4_states.Rs * 2;  %Умножение скорости передачи символов на 2, так как мы имеем 2 бита на символ

 

  tx = fsk4_mod(fsk4_states,test_bits);

 

  %добавление шума

  %EsNo = 10^(aEsNodB/10);

  %EbNo = EsNo

  %variance = Fs/(Rb*EbNo);

  %nsam = length(tx);

  %noise = sqrt(variance/2)*(randn(1,nsam) + j*randn(1,nsam));

  %rx    = tx*exp(j*pi/2) + noise;

  rx    = tx;

  rx    = rx(timing_offset:length(rx));

 

  [rx_bits biterr rxphi] = fsk4_demod_fmrid(fsk4_states,rx);

  ber = 1;

 

  %нужно учитывать смещение входных данных для вывода данных

  %пока не обнаружена преамбула

  ox = 1;

  for offset = (1:100)

    nerr = sum(xor(rx_bits(offset:length(rx_bits)),test_bits(1:length(rx_bits)+1-offset)));

    bern = nerr/(bitcnt-offset);

    if(bern < ber)

      ox = offset;

      best_nerr = nerr;

    end

    ber = min([ber bern]);

  end

  offset = ox;

  disp(ncoarse, timing); %d nerr: %d\n", offset, best_nerr);

 

end

 

Функция fine_ex – в данном пункте кода оценивается точная синхронизация, с учётом смещений происходит оценка параметров синхронизации.

 

function fsk4_rx_phi(socket)

  %pkg load parallel

  offrange = [100:200];

  [a b c phi] = pararrayfun(1.25*nproc(),@nfbert,10*length(offrange),offrange);

 

  close all;

  figure(1);

  clf;

  plot(offrange,phi);

end

 

Функция fsk4_rx_phi – был запущен широкий диапазон смещений и по графику удалось убедиться в том, что точная синхронизация имеет смысл.

 

function fsk4_ber_curves

  global dmr_info; %глобальная переменная

  global nxdn_info;

  global nflt_info;

 

  EbNodB = 1:20;

  bers_tco  = ones(1,length(EbNodB)); %= bers_real = bers_tnco = bers_idealsim

 

  %vectors of the same param to pass into pararrayfun

  dmr_infos = repmat(dmr_info,1,length(EbNodB));

  nflt_infos = repmat(nflt_info,1,length(EbNodB));

  thing = @fsk4_demod_thing;

 

  %Прекрасные инновации Брэди , чтобы использовать все ядра и реально ускорить процесс моделирования

 

  %практика

   % параллельная нагрузка

    bers_idealsim = pararrayfun(floor(1.25*nproc()),@nfbert,EbNodB,nflt_infos);

    [bers_real,bers_tco,bers_tnco] = pararrayfun(floor(1.25*nproc()),@nfbert,EbNodB,dmr_infos);

  %catch

  %  printf("You should install package parallel. It'll make this run way faster\n");

  %  for ii=(1:length(EbNodB));

      %[bers_real(ii),bers,tco(ii),bers_tnco(ii)] = nfbert(EbNodB(ii));

  %  end

  %end_try_catch

 

  close all

  figure(1);

  clf;

  semilogy(EbNodB, bers_tnco,'r;4FSK non-coherent theory;')

  hold on;

 

  semilogy(EbNodB, bers_tco,'b;4FSK coherent theory;')

  semilogy(EbNodB, bers_real ,'g;4FSK DMR simulation;')

  semilogy(EbNodB, bers_idealsim, 'v;FSK4 Ideal Non-coherent simulation;')

  hold off;

  grid('minor');

  axis([min(EbNodB) max(EbNodB) 1E-5 1])

  legend('boxoff');

  xlabel('Eb/No (dB)');

  ylabel('Bit Error Rate (BER)')

 

end

 

При запуске функции fsk4_ber_curves сравниваем теоретическую производительность модема 4FSK с нашим моделированным DMR модемом. Чтобы ускорить процесс моделирования вводится параллельная нагрузка. Получаем в результате графики зависимости реальной вероятности битовой ошибки от ОСШ.

 В статье "Исследование принципа построения L1 уровня цифровой мобильной связи DMR", часть 3 будут рассмотрены результаты моделирования и выполнения кода.

 

Список используемой литературы:

1.Simulating%20the%20DMR%20Modem%20_%20Rowetel.html

2.ETSI TS 102 361-1 V2.2.1 (2013-02) "Electromagnetic compatibility and Radio spectrum Matters (ERM); Digital Mobile Radio (DMR) Systems; Part 1: DMR Air Interface (AI) protocol"