Разработка системы определения местоположения устройства

Размещено на http://www.allbest.ru/

РАЗРАБОТКА СИСТЕМЫ ОПРЕДЕЛЕНИЯ МЕСТОПОЛОЖЕНИЯ УСТРОЙСТВА

Дубовик Николай Николаевич

Туманов Владислав Михайлович



Алгоритмическое обеспечение системы определения местоположения устройства

Первым этапом разработки программного обеспечения системы является анализ математических методов определения местоположения устройства внутри здания. Следующим этапом реализации системы является алгоритмическое проектирование и разработка программного обеспечения.

Алгоритмы, представленные ниже наглядно показывают то, как работают ключевые элементы разрабатываемой системы.

Рисунок 1. Алгоритм работы фазового метода определения времени распространения сигнала

Первый основной алгоритм -- алгоритм фазового метода определения времени распространения сигнала. На первом этапе после инициализации алгоритма проверяется, получены ли данные в виде CSI характеристик от каждой из антенн на всех установленных частотах. Если они не получены, то начинаются прыжки между частотами. Останавливаясь на каждой непроверенной частоте, снимаются характеристики CSI с требуемого количества пакетов, вычисляется фаза прибытия и записывается соответствующее уравнение для вычисления времени. После того, как все частоты проверены, собирается системы уравнений, полученных на этих частотах, и вычисляется искомое время распространения сигнала. Блок-схема алгоритма представлена на рисунке 1.

Другим основным алгоритмом является алгоритм поиска угла прихода сигнала. Блок схема алгоритма представлена на рисунке 2.

Рисунок 2. Алгоритм определения угла



Алгоритм работает следующим образом. На первом этапе работы алгоритм получает характеристики CSI, собранные ранее в момент переключения между частотами. Данные от каждой антенны о полученных ей пакета на всех соответствующих поднесущих записываются в общую матрицу CSI. Затем происходит сглаживание и обработка матрицы в соответствии с математическими формулировками.

Рисунок 3. Алгоритм работы системы определения местоположения устройства внутри здания



После этого вычисляется матрица собственных векторов, к которой затем применяется алгоритм MUSIC. Пики спектра MUSIC являются возможными углами прихода сигнала. На последнем этапе возможные углы объединяются в группы, и затем в соответствии с ранее вычисленным временем прихода сигнала выбирается наиболее вероятная группа углов.

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

Программное обеспечение системы определения местоположения устройства

После того, как разработаны основные алгоритмы работы системы, следует их программная реализация. Поскольку одной из целей работы является разработка и исследования прототипа системы определения местоположения устройства внутри различных зданий, то реализация готового устройства и соответствующего программного продукта не требуется. Для удобства исследований и настройки, а также для упрощения реализации, в качестве основной программной среды выбран пакет MATLAB. Некоторые части системы для простоты реализации написаны на языке C++.

Для получения характеристик CSI на физическом уровне используется популярный в подобных исследованиях инструмент Linux 802.11 CSI Tool, развернутый на ОС Ubuntu Linux, ядро 3.5.7. Так же для обеспечения работы протокола смены частоты, и для использования таймера высокого разрешения был внедрен программный патч для драйвера iwlwifi.

Ниже рассматривается программная реализация основных функций, использующихся во время работы системы.

Изменение несущей частоты по стандарту Wi-Fi

Для того чтобы точно определить время распространения, реализуемая система снимает данные CSI на различных частотах. В обычном режиме устройства меняют частоты только в определенных случаях. Однако стандарт 802.11 регламентирует как наличие, так и состав служебных пакетов, которые позволяют как передатчику, так и приемнику перемещаться с одной частоты на другую. Благодаря программному патчу, драйвер iwlwifi в соответствии с алгоритмом работы самостоятельно посылает пакет с информацией о смене канала и частоты.

Фрейм анонса смены канала использует формат фрейма Action, и показан на рисунке 4

Рисунок 4. Формат фрейма анонса смены канала

Поле Category устанавливается в 0 (соответствует работе со спектром).

Поле Action Value устанавливается в 4 (соответствует фрейму анонса смены канала).

Поле Channel Switch Announcement element показано на рисунке 5.



Рисунок 5. Формат элемента анонса смены канала

Поле Length устанавливается в 3.

Поле Channel Switch Mode указывает любые ограничения на передачу до переключения канала. Точка доступа должна установить поле режима переключения каналов на 0 или 1 при передаче, в зависимости от условий. Поле New Channel Number должно быть установлено на номер канала, на который происходит переключение.

Поле Channel Switch Count указывает количество отправленных пакетов до тех пор, пока устройство, отправляющее элемент объявления переключателя канала, не переключится на новый канал; или должно быть установлено в 0. Значение 1 указывает, что смена должна произойти непосредственно перед следующим пакетом. Значение 0 указывает, что смена может произойти в любое время после передачи фрейма, содержащего элемент.

Решение задачи определения расстояния

Для поиска решения, удовлетворяющего всем временным уравнениям, полученным алгоритмом перемещения по частотам используется китайская теорема об остатках. Листинг псевдокода, реализующего поиск решения по данной теореме представлен в таблице 1.

Листинг псевдокода решения задачи определения расстояния

Дано: зависимость t от частоты f и фазы h

i пар фаза/частота

Вывод: x являющееся наименьшим общим кратным для всех частот

function chinese_remainder(int *n, int *a, int len)

{

float p, prod = 1, sum = 0;

int i = 1;

for (i = 0; i < len; i++) prod *= n[i];

for (i = 0; i < len; i++) {

p = prod / n[i];

sum += a[i] * remainderf (p, n[i]) * p;

}

return sum % prod;

}

Получение обратной NDFT

Для того, чтобы определить верное время прихода сигнала в условиях многолучевого распространения воспользуемся свойствами обратного дискретного преобразования Фурье. Так как несущие частоты распределены неравномерно, то необходимо использовать специальное преобразование -- NDFT, которое затем инвертируется. Для обеспечения конечного числа последовательностей вводятся дополнительные ограничения. Искомое время распространения в таком случае будет являться первым пиком спектра. Листинг псевдокода программы, реализующей получение обратного NDFT представлен в таблице 2.

Листинг псевдокода нахождения обратного NDFT

Дано: измеренные каналы h

Неравнораспределенная матрица дискретного преобразования фурье F

б -- параметр разреженности; е -- параметр сходимости

Вывод: инвертированная NDFT, p

P[0] = 0, t = 0, г = 1/||F||2

while converged = false do

p[t+1] =SPARSIFY(p[t] ? гF?(Fp[t] ? h), гб)

if ||p[t+1] -- p[t]||2 < е then

converged = true

p = p[t+1]

else

t = t + 1

end if

end while

function SPARSIFY(p, t)

for i = 1, 2, ...length(p) do

if |p[i]| < t then

p[i] = 0

else

p[i] = p[i]*|p[i]|?t/|p[i]|

end if

end for

end function

Реализация алгоритма MUSIC

Для получения из характеристик CSI значений угла прихода сигнала используется функционал программного пакета MATLAB, а именно функция pmusic.

Функция задается как [S, w] = pmusic (x, p) и реализует алгоритм MUSIC (Multiple Signal Classification) и возвращает S, оценку псевдоспектра входного сигнала x и вектор w нормированных частот, по которым оценивается псевдоспектр. Псевдоспектр рассчитывается с использованием оценок собственных векторов корреляционной матрицы, связанных с входными данными x, где x задается как вектор строки или столбца, представляющий одно наблюдение сигнала или как прямоугольный массив, для которого каждая строка x представляет собой отдельное наблюдение за сигналом (например, каждая строка является одним выходом массива датчиков), так что x'*x является оценкой корреляционной матрицы.

Второй входной аргумент p представляется или как скалярное целое число, в этом случае размер подпространства сигнала равен p; или как двухэлементный вектор, в этом случае p(2), второй элемент p, представляет собой порог, который умножается на лmin, наименьшее оцениваемое собственное значение корреляционной матрицы сигнала. Собственные значения ниже порога лmin*p(2) присваиваются шумовому подпространству. В этом случае p(1) задает максимальную размерность сигнального подпространства.

Значения S и w имеют одинаковую длину. В общем случае длина FFT и значения входа x определяют длину вычисленного S и диапазон соответствующих нормированных частот.

Листинг кода в MATLAB, реализующий функцию pmusic представлен в таблице 3.

Листинг функции pmusic

function varargout = pmusic( x, p, varargin )

narginchk(2, 10);

[md, msg, msgobj] = music(x, p, varargin{:});

if ~isempty(msg), error(msgobj); end

% compute the pseudospectrum

range = 'half';

if strcmpi(md.range, 'twosided') || (length(md.nfft)) > 1, range = 'whole'; end

[Sxx, w] = pseudospectrum(md.noise_eigenvects, md.eigenvals, md.nfft, md.Fs, ...

range, md.EVFlag);

% Assign outputs

if nargout==0

titlestring = 'Pseudospectrum Estimate via ';

if md.EVFlag, titlestring = [titlestring, 'Eigenvector Method'];

else titlestring = [titlestring, 'MUSIC'];

end

w = {w};

if ~isempty(md.Fs), w = [w {'Fs', md.Fs}]; end

hps = dspdata.pseudospectrum(Sxx, w{:}, 'SpectrumRange', range);

if md.centerdc

centerdc(hps);

end

plot(hps);

title(titlestring);

ylabel('Power (dB)');

return

end

% center dc if needed

if md.centerdc

[Sxx, w] = psdcenterdc(Sxx, w, [], md);

end

% Cast to enforce precision rules

if isa(Sxx, 'single')

w = single(w);

end

switch nargout

case 1, varargout = {Sxx};

case 2, varargout = {Sxx, w};

case 3, varargout = {Sxx, w, md.noise_eigenvects};

case 4, varargout = {Sxx, w, md.noise_eigenvects, md.eigenvals};

end

function [Sxx, w] = pseudospectrum(noise_eigenvects, eigenvals, nfft, Fs, range, EVFlag)

% compute weights

if EVFlag,

% Eigenvector method, use eigenvalues as weights

weights = eigenvals(end-size(noise_eigenvects, 2)+1:end); % Use the noise subspace eigenvalues

else

weights = ones(1, size(noise_eigenvects, 2));

end

% Compute the freq. response of each noise subspace eigenfilter via freqz.

% Use freqz to compute the freq vector only when computing the "whole"

% spectrum.

den = 0;

for n = 1:size(noise_eigenvects, 2),

% Don't call freqz with Fs=[], because it will default Fs to 1! Use

% Fs={}, which gets ignored in varargin.

if isempty(Fs), sampleRate = {}; else sampleRate{1} = Fs; end

% Compute the freq vector directly in Hz to avoid roundoff errors later.

% Only the "whole" freq vector computed by freqz is valid for spectra.

[h, w] = freqz(noise_eigenvects(:, n), 1, nfft, 'whole', sampleRate{:});

den = den + abs(h).^2./weights(n);

end

s = 1./den; % This is the pseudospectrum

% Calculate the spectrum over the frequency range requested by the user.

[Sxx, w] = computespectrumrange(s, w, range, nfft);

function [Sxx, w] = computespectrumrange(Sxx, w, range, nfft)

%Return the correct spectrum range based on the user input argument RANGE.

% Convert input row vectors to columns (if not a matrix).

if any(size(Sxx)==1), Sxx = Sxx(:); end

w = w(:);

if strcmpi(range, 'half'),

if rem(nfft, 2), select = 1:(nfft+1)/2; % ODD

else select = 1:nfft/2+1; % EVEN

end

Sxx = Sxx(select,:); % Take only [0, pi] or [0, pi)

w = w(select);

end

В результате работы функции оценивается псевдоспектр сигнала, пик которого соответствует углу прихода сигнала.

Задача определения позиции

навигация местоположение устройство программный

По итогам работы предыдущих алгоритмов и функций система определяет несколько различных возможных координат местонахождения искомого устройства, лежащих в некоторой окрестности. Для того, чтобы определить итоговые наиболее вероятные координаты, используется оптимизация по методу наименьших квадратов, которая так же ограничивается геометрическими условиями. Для решения этой задачи оптимизации используется функция fmincon программного пакета MATLAB.

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



Где b и beq -- векторы, A и Aeq -- матрицы, c(x) и ceq(x) -- функции, возвращающие векторы, а f(x) -- функция, возвращающая скаляр. F(x), c(x) и ceq(x) могут быть нелинейными функциями, x, lb и ub могут быть введены как векторы или матрицы.

X = fmincon (fun, x0, A, b) начинает с x0 и пытается найти минимальное значение x функции, описанной в fun с учетом линейных неравенств A*x ? b, при этом x0 может быть скаляром, вектором или матрицей.

В результате работы этой функции выводятся значения координат, которые система принимает как искомые координаты местоположения устройства, которые требовалось определить.



Выводы

Статья посвящена алгоритмической и программной реализации системы определения местоположения устройства внутри зданий.

В начале статьи описываются основные алгоритмы работы системы, а именно алгоритм определения времени распространения сигнала с помощью измерений на различных частотах, алгоритм определения угла прихода сигнала. Затем дан общий алгоритм системы определения местоположения устройства.

Во второй части статьи показаны основные аспекты программной реализации системы, используемые программные пакеты и функции. В качестве среды для обработки и исследований сигналов выбран пакет MATLAB, недостающий функционал реализован с помощью программ на языке C++. Так же дано описание стандартных фреймов Wi-Fi, с помощью которых обеспечивается смена каналов и частот.

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

Размещено на Allbest.ru

...