Гидроакустический модем в Матлаб

Гидроакустический модем в МатлабГидроакустический модем на сегодняшний день довольно избитая тема. В промышленном плане существует множество готовых решений для различных целей обмена цифровой информацией под водой. Но все решения довольно дорогостоящие, используют классические способы для предачи и приема цифровой информации под водой, например, такие, как BPSK, ASK и др. В данном посте речь пойдет о матлаб модели простейшего ASK гидроакустического модема с нестандартным способом оцифровки данных, позволяющим обойтись без обязательной петли ФАПЧ. В матлаб модели способ оцифровки данных на приеме надежно работает при SNR > 3 и доплеровским сдвигом соответствующим скорости движения передатчика относительно приемника в пределах 60 узлов. При этом скорость обмена информацией может доходить до 2000 бод в секунду.

Матлаб модель устройства предполагает обмен короткими сообщениями объемом 4...8 бит в каждом. Для передачи кодированной информации от одного устройства на другое формируется протокол, в котором содержится один стартовый бит, четыре бита данных и один стоповый бит. Предполагается что усилительное устройство приемника усиливает принимаемый сигнал до состояния ограничения, что дает возможность пропустить его через компаратор и работать с сигналом как с цифровым. Вот как выглядит модель зашумленного сигнала 0ХD5.

А вот сигнал на выходе компаратора.

На данной картинке по сути вообще ничего не понятно. Не понятно где находится старт, где стоп, где собственно биты информации. И тем не менее все прекрасно определяется. Резултат сравнения входной информации и выходной "bitend Ok"

В чем же хитрость? Хитрость в алгоритме или методике подсчета импульсов тактовой или несущей частоты во временных пределах соответствующих импульсам старт, данных и стоп. Предлагаю просто запустить приведенный ниже код в Matlab и увидеть самостоятельно как все крутится. Код снабжен комментариями. Разобраться в алгоритме будет просто.


%% BodyWatcherTest

%% Исходные данные
clear;

%bit_stream = (rand(1, 8) > 0.5);
bit_stream = [1 1 0 1 0 1 0 1];
N = length(bit_stream);

% Несущая частота Гц
f = 36000;
% Кол. выб. на период несущей 
nsp = 16; 
% Частота семплирования Гц 
fs = f*nsp;
% Частота Найквиста
fn = fs/2;    
% Частота доплера Гц
fdop = -1200;
% время выборки
Ts = 1/fs;
% Амплитуда
A = 1;
% количество периодов несущ.для одного бита
Tb = 20;
% количество периодов несущ.для задержки
Td = 0;
% количество бит между старт-стоп
Nb = 8;
% Символьная частота с доплером Гц
fsimd = (f+fdop)/Tb;
fsim = f/Tb;
% частота задержки с долером Гц
fdeld = (f+fdop)/Td;
fdel = f/Td;
% Результир. симв. частота с доплером Гц
fresd = (f+fdop)/(Tb + Td);
fprintf('fresd %d\n', fresd);
% Количество выборок необх. для одного символа с доплером 
tsymbold = f*nsp/fsimd;
tsymbol = f*nsp/fsim;
% Количество выборок необх. для задержки с доплером
tdelayd = f*nsp/fdeld;
tdelay = f*nsp/fdel;
% Количество выборок необх. для символа+задержки с доплером
ttotald = tsymbold + tdelayd;
ttotal = tsymbol + tdelay;
fprintf('ttotald %d\n', ttotald);
% количество точек моделирования
Nl = 3*Tb*nsp + 2*Tb*8*nsp + 2*Tb*nsp;
fprintf('Nl %d\n', Nl);

BodySignalBeg = zeros(1,Nl); 
BodySignalDop = zeros(1,Nl); 

%% Модуляция

t = 0;

% пробел
for i = 1:3
  for j = 1:tsymbold
    BodySignalBeg(t+1) = 0;
    BodySignalDop(t+1) = 0;
    t = t + 1;
  end 
end

% старт
for i = 1:1
  for j = 1:tsymbold
    BodySignalBeg(t+1) = A*sin(2*pi*f*Ts*t);
    BodySignalDop(t+1) = A*sin(2*pi*(f+fdop)*Ts*t);
    t = t + 1;
  end 
  % задержка
  for j = 1:tdelayd
    BodySignalDop(t+1) = 0;
    BodySignalBeg(t+1) = 0;
    t = t + 1;
  end 
end

% данные
for i = 1:Nb
  for j = 1:tsymbold
    if bit_stream(i) == 1
      BodySignalDop(t+1) = A*sin(2*pi*(f+fdop)*Ts*t);
      BodySignalBeg(t+1) = A*sin(2*pi*f*Ts*t);
    % нулевой сигнал  
    else
      BodySignalBeg(t+1) = 0;
      BodySignalBegDop(t+1) = 0;
    end 
    t = t + 1;
  end 
  % задержка
  for j = 1:tdelayd
    BodySignalBeg(t+1) = 0;
    BodySignalDop(t+1) = 0;
    t = t + 1;
  end 
end

% стоп
for i = 1:1
  for j = 1:tsymbold
    BodySignalDop(t+1) = A*sin(2*pi*(f+fdop)*Ts*t);
    BodySignalBeg(t+1) = A*sin(2*pi*f*Ts*t);
    t = t + 1;
  end 
  % задержка
  for j = 1:tdelayd
    BodySignalBeg(t+1) = 0;
    BodySignalDop(t+1) = 0;
    t = t + 1;
  end 
end

%% Канал

% отношение сигнал/шум в децибелах
SNRawgn = 12;
% смесь мод.сигнала с шумом
BodySignalAwgn = awgn(BodySignalDop, SNRawgn);

% фактическая мощность сигнала 
Sbpsk = mean(BodySignalDop.^2);
fprintf('Sbpsk %d\n', Sbpsk);
% фактическая мощность зашумленного сигнала 
SNbpsk = mean(BodySignalAwgn.^2);
fprintf('SNbpsk %d\n', SNbpsk);
% фактическая разнсть мощностей сигнала и заш. сигн 
Nss = mean((BodySignalDop - BodySignalAwgn).^2);
fprintf('Nss %d\n', Nss);
% фактический SNR
SNR = SNbpsk/Nss;
fprintf('SNR %d\n', SNR);

BodySignalAwgn = 3*BodySignalAwgn;

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

BodySignal = BodySignalAwgn;
%BodySignal = BodySignalAwgnFiltr;

for i=1:Nl
  if BodySignal(i) > 3 
    BodySignal(i) = 3;
  elseif BodySignal(i) < -3 
    BodySignal(i) = -3;
  end
end


%% Компаратор

BodySignalComp = reshape( BodySignal > 1.5, [], 1 );


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

ignoreData = 1;
step = 0;
start = 0;
count = 0;
count_test = 0;
fl_start = 0;
t = 1;
tn = 1;
bit = zeros(1,N);
bitend = zeros(1,N);
bit_block = zeros(1,N+1);
step_start = 0;
bd_test = zeros(1,N);
bit_test = zeros(1,N);
count_end = Tb*nsp/10;
cnt = 0;
n = 0;
% количество выборок шума при регистр старта
noise = 7;
cnt_smb = 0;
cnt_smb_en = 0;

for step = 1:Nl

  % фиксация старта
  if (ignoreData)
    % начальные ворота
    if cnt_smb_en 
      cnt_smb = cnt_smb + 1;
      if cnt_smb > tsymbol/4 
        count = 0;
        cnt_smb = 0;
      end  
    end
    % Ждем старта
    if BodySignalComp(step) == 1 
      count = count + 1;
      % разрешение 
      if cnt_smb_en == 0 
        cnt_smb_en = 1;
      end
      if count > noise 
        % Стартовый символ прошел
        ignoreData = 0;
        % флаг старта в 1
        fl_start = 1;
        % фиксация точки старта
        step_start = step;
        count = 0;
      end
    end
  end

  % начало побайтной синхронизации 
  if ignoreData == 0 && t < 9;
    cnt = cnt + 1;
    if mod(cnt, ttotal) == 0 
      n = n + 1;
    end
  end

  % битовые ворота
  if ignoreData == 0 && step > step_start + ttotal
    if step > (step_start + n*ttotal) && step < (step_start + n*ttotal + tsymbol)
      bd_test(step) = 1;
    else 
      bd_test(step) = 0;
    end
  end
  
  % декодер
  if step > (step_start + t*ttotal) && step < (step_start + (t+1)*ttotal) && bit_block(t) == 0 && t < 9 && ignoreData == 0
    if BodySignalComp(step) == 1  
      count = count + 1;
      count_test = count_test + 1;
      if count > count_end 
        bit_block(t) = 1;
        bit(t) = 1;
        bit_test(step) = 1;
        count = 0;
        t = t + 1;
      end
    elseif step > (step_start + t*ttotal + ttotal/2)
      bit(t) = 0;
      bit_block(t) = 1;
      bit_test(step) = 0;
      count = 0;
      t = t + 1;
    end
  % стоп  
  elseif t == 9   
    if BodySignalComp(step) == 1  
      count = count + 1;
      % загрузка результата в bitend
      if count > count_end 
        bitend = bit;
      end
    end
  end

end

%fprintf('start %d\n', start);
fprintf('step_start %d\n', step_start);
%fprintf('bitend %d\n', bitend);
%fprintf('bit_block %d\n', bit_block);
fprintf('count_test %d\n', count_test);
fprintf('count_end %d\n', count_end);
%fprintf('bit_stream %d\n', bit_stream);
if bitend == bit_stream 
  fprintf('bitend Ok \n');
else
  fprintf('bitend Error \n');
end

%% Графика

figure;
plot(BodySignalAwgn, '-b', 'LineWidth', 2);
hold on;
grid on;
plot(BodySignalAwgnFiltr, 'r-', 'LineWidth', 2);
hold on
hold off;
axis([0 Nl -10 10]),xlabel('Выборки'),ylabel('Амплитуда')
title('Сигнал+Шум');legend('Сигнал+Шум', 'Сигнал+Шум+Фильтр');

figure;
plot(BodySignal, '-b', 'LineWidth', 2);
hold on;
grid on;
% plot(BodySignalAwgnFiltr, 'r-', 'LineWidth', 2);
% hold on
% hold off;
axis([0 Nl -5 5]),xlabel('Выборки'),ylabel('Амплитуда')
title('Сигнал+Шум');legend('Сигнал+Шум');

% figure;
% plot(BodySignalBeg, 'b-', 'LineWidth', 2);
% hold on;
% grid on;
% plot(BodySignalDop, 'r-', 'LineWidth', 2);
% hold on
% hold off;
% axis([0 Nl -7 7]),xlabel('Выборки'),ylabel('Амплитуда')
% title('BodySignal доплер');

% Фильтр полосовой на 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
% axis([0 8e4 -80 40]);
% xlabel('частота/Hz');ylabel('амплитуда/dB');title('АЧХ полосового фильтра');

figure;
plot(BodySignalComp, '-g', 'LineWidth', 2);
hold on;
grid on;
plot(bd_test, 'r-', 'LineWidth', 2);
hold on
plot(bit_test, 'b-', 'LineWidth', 2);
hold on
hold off;
axis([0 Nl -1 2]),xlabel('Выборки'),ylabel('Амплитуда')
title('Выход компаратора'); legend('компаратор', 'битовые ворота', 'биты');
	

Переведя все в цифру, тоесть оцифровав принятый сигнал, можно построить обработку на более высоком уровне, применив цифровые фильтры. Результат - работа устройства при соотношении сигнал/шум меньше 3. Выглядит это примерно так.

P.S. На базе данного алгоритма выполнен лабораторный макет устройства без акустической части. Испытания дали положительный результат. Таким образом можно строить недорогие гидроакустические модемы для обмена цифровой информацией под водой при больших скоростях перемещения передатчика относительно приемника. Это могут быть устройства для управления подводными дронами или устройства типа BodyWatcher и им подобные. Дерзайте, кому не лень... И не забывайте оставлять свои комменты. Кому пригодится данный алгоритм не забывайте об интересах автора.

Top.Mail.Ru