Цифровая обработка сигналов

    • Alexey Antipin
      17.07.2014
      10:38

      Измерение вектора ошибки для системы стандарта IEEE 802.15.4 (ZigBee)

      Данный пример демонстрирует, как использовать системный объект COMM.EVM для измерения величины векторной ошибки (EVM) при симуляции передатчика стандарта IEEE® 802.15.4 [ 1 ]. IEEE 802.15.4 - это основной документ спецификации ZigBee.

      Contents

      Величина векторной ошибки

      Велечина векторной ошибки (error vector magnitude, EVM) является мерой различия между идеальным модулированным сигналом и реально переданным сигналом. EVM используется для колличественной оценки точности модуляции передатчиком. Согласно [ 1 ], среднеквадратическая EVM передатчика стандарта 802.15.4 не должна превышать 35%.

      Параметры системы

      Система 802.15.4 для полосы 868 МГц применяет технику расширения спектра методом прямой последовательности (direct sequence spread spectrum, DSSS) и относительную двоичную фазовую манипуляцию (DBPSK) для модуляции последовательности чипов.

      dataRate = 20e3;   % Скорость битового потока в Гц
      M = 2;             % Порядок модуляции (BPSK)
      chipValues = [1;1;1;1;0;1;0;1;1;0;0;1;0;0;0];
                         % Значения чипов для бита равного 0
                         % Значения чипов для 1 имеют противоположные значения
      

      Раздел 6.7.3 из [ 1 ] определяет, что измерения проводятся для 1000 отсчетов синфазной (I) и квадратурной (Q) составляющих сигнала в основной полосе. Для учета задержки, вносимой фильтром, при симуляции передаваемых символов мы добавили 1 дополнительный бит. Передаваемый сигнал интерполируется с коэффициентом интерполяции равным 4. Мы предполагаем отношение сигнал/шум (SNR) равное 60 дБ для учета искажений вносимых передатчиком и тестируемым оборудованием.

      numSymbols = 1000;          % Число символов, необходимое для одного измерения EVM
      numFrames = 100;            % Число фреймов
      nSamps = 4;                 % Число отсчетов, представляющих символ
      filtSpan = 8;               % Длина фильтра в символах
      gain = length(chipValues);  % Усиление за счет расширения (число чипов в символе)
      chipRate = gain*dataRate;   % Скорость потока чипов
      sampleRate = nSamps*chipRate;    % Окончательная частота дискретизации
      numBits = ceil((numSymbols)/gain)+1;
                                  % Число бит, необходимое для одного измерения EVM
      SNR = 60;                   % Отношение сигнал/шум в дБ
      

      Инициализация системы

      Мы можем получить символы с BPSK модуляцией применив простое преобразование: 0 --> +1 и 1 --> -1. Если мы выполним такое же преобразование для последовательности чипов, то получим модулированный сигнал еще до преобразования битов в чипы. В этом случае появляется возможность использовать матричную математику и за счет этого получить эффективный код MATLAB. Спецификация ZigBee также определяет характеристики формирующего импульсы фильтра. Это должен быть фильтр типа приподнятого косинуса с коэффициентом сглаживания (roll-of factor) равным 1.

      % Преобразование последовательности чипов
      chipValues = 1 - 2*chipValues;
      
      % Создание фильтра типа приподнятого косинуса с коэф. сглаживания равным 1
      hTxFilter = comm.RaisedCosineTransmitFilter('RolloffFactor', 1, ...
        'OutputSamplesPerSymbol', nSamps, ...
        'FilterSpanInSymbols', filtSpan);
      hRxFilter = comm.RaisedCosineReceiveFilter('RolloffFactor', 1, ...
        'InputSamplesPerSymbol', nSamps, ...
        'FilterSpanInSymbols', filtSpan, ...
        'DecimationFactor', nSamps);
      

      Измерение величины векторной ошибки

      Расширение Communications System Toolbox™ предоставляет системный объект COMM.EVM для расчета среднеквадратического, максимального и процентного значения EVM. Раздел 6.7.3 документа [ 1 ] определяет метод расчета EVM, при котором средняя ошибка по I и Q отсчетам нормируется к мощности символа. Для сигнала с BPSK мощности символом сигнального созвездия совпадают, поэтому мы можем использовать опцию нормировки 'Peak constellation power' (пиковая мощность сигнального созвездия). Доступные также и другие опции нормировки, которые могут использоваться в других стандартах систем связи. Это средняя мощность сигнального созвездия и средняя мощность эталонного сигнала.

      hEVM = comm.EVM('Normalization', 'Peak constellation power')
      
      hEVM = 
      
        System: comm.EVM 
      
        Properties:
                     Normalization: 'Peak constellation power'
            PeakConstellationPower: 1                         
              MaximumEVMOutputPort: false                     
          XPercentileEVMOutputPort: false                     
       X 
      

      Симуляция

      Первым делом генерируется последовательность случайных бит, затем выполняется дифференциальное кодирование с помощью системного объекта DifferentialEncoder и модуляция BPSK. Расширение спектра модулированных символов осуществляется путем матричного умножения на вектор содержащий последовательность чипов. Затем символы расширенной последовательности проходят через фильтр типа приподнятого косинуса. Объект EVM предполагает, что принятые символы (rd) и идеальные символы (с) синхронизированы и имеют одинаковую частоту дискретизации, поэтому мы понижаем частоту дискретизации принятого сигнала ® и синхронизируем его с идеальным сигналом (с).

      [ 1 ] требует, чтобы для расчета одного значения среднеквадратической EVM использовалось 1000 символов. Чтобы получить точные результаты мы выполняем симуляцию 100 фреймов, каждый из которых содержит 1000 символов, и выбираем максимальное значение из 100 полученных измерений EVM. Можно заметить, что разработанный передатчик удовлетворяет требованиям, описанным в предыдущей секции Величина векторной ошибки

      % Задержки, вносимые фильтрами передатчика и приемника, одинаковы и равны
      % половине длины фильтра. Общая задержка системы равняется их сумме, что
      % равносильно длине одного фильтра.
      refSigDelay = hTxFilter.FilterSpanInSymbols;
      
      % Число символов в одном фрейме
      simNumSymbols = numBits*gain;
      
      % Начальное значение пиковой среднеквадратической EVM
      peakRMSEVM = -inf;
      
      % Создание объекта comm.DifferentialEncoder для дифференциального
      % кодирования
      hdiffenc = comm.DifferentialEncoder;
      
      % Создание объекта comm.AWGNChannel и установка опции NoiseMethod
      % в значение 'Signal to noise ratio (SNR)'
      hChan = comm.AWGNChannel('NoiseMethod', 'Signal to noise ratio (SNR)',...
        'SNR', SNR);
      
      % Цикл обработки пакетов
      for p=1:numFrames
          % Создание случайно последовательности бит
          b = randi([0 M-1], numBits, 1);
          % Дифференциальное кодирование
          d = step(hdiffenc, b);
          % Модуляция
          x = 1-2*d;
          % Преобразование символов в чипы (расширение)
          c = reshape(chipValues*x', simNumSymbols, 1);
          % Формирование импульсов
          cUp = step(hTxFilter, c);
          % Расчет и настройка опции 'SignalPower' для объекта канала
          hChan.SignalPower = sum(cUp.^2)/length(cUp);
          % Добавление шума
          r = step(hChan, cUp);
          % Понижение частоты дискретизации сигнала. Учет задержек, вносимых
          % фильтрами.
          rd = step(hRxFilter, r);
          % Измерение с помощью системного объекта EVM
          rmsEVM = step(hEVM, ...
            complex(rd(refSigDelay+(1:numSymbols))), ...
            complex(c(1:numSymbols)));
          % Обновление пикового значения среднеквадратической EVM
          if (peakRMSEVM < rmsEVM)
              peakRMSEVM = rmsEVM;
          end
      end
      
      % Вывод результатов
      fprintf(' Наихудшее значение RMS EVM (%%): %1.2f\n', peakRMSEVM)
      
       Наихудшее значение RMS EVM (%): 0.19
      

      Комментарии

      Здесь рассмотрено, как использовать объект COMM.EVM для тестирования передатчика ZigBee, в соответствии с указаниями стандарта по расчету EVM. В данном примере использовалась грубая модель, учитывающая только аддитивный белый шум Гаусса, которая показала, что измеренное значение EVM меньше указанных в стандарте 35%.

      Библиография

      1. IEEE Standard 802.15.4, Wireless Medium Access Control (MAC) and Physical Layer (PHY) Specifications for Low-Rate Wireless Personal Area Networks, 2003.
    • Alexey Antipin
      10.02.2014
      18:04
      Начиная с релиза MATLAB R2013a в состав Phased Array System Toolbox входят три графических приложения, позволяющих анализировать сенсорные решетки и радиолокационные сигналы, а также оценивать основные характеристики радаров. С их помощью пользователи могут сразу же приступить к работе с инструментами данного расширения, не тратя времени на изучение командного интерфейса.

      Загрузка плеера
    • Alexey Antipin
      04.02.2014
      14:35

      Данный пример демонстрирует применение адаптивных фильтров для ослабления акустического шума в системах активного шумоподавления.

      Оглавление:

      Активное шумоподавление.

      Системы активного шумоподавления (active noise control) применяются для ослабления распространяющегося по воздуху нежелательного шума с помощью электроакустических приборов: измерительных устройств (микрофонов) и возбудителей сигнала (динамиков). Шумовой сигнал обычно исходит от некоторого устройства, например вращающегося механизма, и имеется возможность измерить шум рядом с его источником. Целью системы активного шумоподавления является создание «анти-шумового» сигнала с помощью адаптивного фильтра, который ослабит шум в определенной тихой области. Эта проблема отличается от обычного адаптивного шумоподавления тем, что: - ответный сигнал не может быть тут же измерен, а доступна только его ослабленная версия; - при адаптации система активного шумоподавления должна учитывать вторичную ошибку распространения сигнала от динамиков до микрофона.

      Более детально задачи активного шумоподавления рассмотрены в книге S.M. Kuo и D.R. Morgan, "Active Noise Control Systems: Algorithms and DSP Implementations", Wiley- Interscience, New York, 1996.

      Путь вторичного распространения.

      Путь вторичного распространения – это путь, который проходит «анти-шумовой» сигнал с выхода динамиков до измеряющего ошибку микрофона, находящегося в тихой зоне. Следующие команды описывают импульсную характеристику пути динамик-микрофон с ограниченной полосой 160-2000 Гц и длиной фильтра равной 0.1 с. Для этой задачи активного шумоподавления мы будем использовать частоту дискретизации равную 8000 Гц.

      Fs     = 8e3; % 8 КГц
      N      = 800; % 800 отсчетов на 8 КГц = 0.1 секунды
      Flow   = 160; % нижняя частота среза: 160 Гц
      Fhigh  = 2000; % верхняя частота среза: 2000 Гц
      delayS = 7;
      Ast    = 20; % подавление 20 дБ
      Nfilt  = 8; % порядок фильтра
      
      % Создание полосового фильтра для имитации канала с ограниченной полосой
      % пропускания
      Fd = fdesign.bandpass('N,Fst1,Fst2,Ast',Nfilt,Flow,Fhigh,Ast,Fs);
      Hd = design(Fd,'cheby2','FilterStructure','df2tsos',...
          'SystemObject',true);
      
      % Фильтрация шума для получения импульсной характеристики канала
      H = step(Hd,[zeros(delayS,1); log(0.99*rand(N-delayS,1)+0.01).* ...
          sign(randn(N-delayS,1)).*exp(-0.01*(1:N-delayS)')]);
      H = H/norm(H);
      
      t = (1:N)/Fs;
      plot(t,H,'b');
      xlabel('Время, с');
      ylabel('Значения коэффициентов');
      title('Импульсная характеристика вторичного пути распространения сигнала');
      

      Определение вторичного пути распространения.

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

      ntrS = 30000;
      s = randn(ntrS,1); % синтез случайного сигнал
      Hfir = dsp.FIRFilter('Numerator',H.');
      dS = step(Hfir,s) + ... % случайный сигнал прошедший через вторичный канал
          0.01*randn(ntrS,1); % шум микрофона
      

      Создание фильтра для оценки вторичного пути распространения.

      В большинстве случаев для адекватного управления алгоритмом длительность отклика фильтра, оценивающего вторичный путь распространения, должна быть короче самого вторичного пути. Мы будем использовать фильтр 250 порядка, что соответствует импульсной характеристике длиной 31 мс. Для этой цели подходит любой алгоритм адаптивной КИХ- фильтрации, но обычно используют нормализованный алгоритм нахождения минимальной среднеквадратической ошибки (normalized LMS-алгоритм) ввиду его простоты и устойчивости.

      M = 250;
      muS = 0.1;
      hNLMS = dsp.LMSFilter('Method','Normalized LMS','StepSize', muS,...
          'Length', M);
      [yS,eS,Hhat] = step(hNLMS,s,dS);
      
      n = 1:ntrS;
      plot(n,dS,n,yS,n,eS);
      xlabel('Число итераций');
      ylabel('Уровень сигнала');
      title('Идентификация вторичного пути распространения с NLMS-алгоритма');
      legend('Ожидаемый сигнал','Сигнал на выходе','Сигнал ошибки');
      

      Точность полученной оценки.

      Как точно оценивается импульсная характеристика вторичного пути? Этот график показывает коэффициенты настоящего пути и пути, рассчитанного алгоритмом. Только конец полученной импульсной характеристики имеет неточности. Эта остаточная ошибка не навредит производительности системы активного шумоподавления во время ее работы над выбранной задачей.

      plot(t,H,t(1:M),Hhat,t,[H(1:M)-Hhat(1:M); H(M+1:N)]);
      xlabel('Время, с');
      ylabel('Значения коэффициентов');
      title('Определение импульсной характеристики вторичного пути распространения');
      legend('Действительная','Оцененная','Ошибка');
      

      Основной путь распространения сигнала.

      Путь распространения шума, который должен быть подавлен, может быть также описан с помощью линейного фильтра. Следующие команды генерируют импульсную характеристику пути источник шума-микрофон с ограниченной полосой 200-800 Гц и имеет длительность отклика равную 0.1 с.

      delayW = 15;
      Flow   = 200; % нижняя частота среза: 200 Hz
      Fhigh  = 800; % верхняя частота среза: 800 Hz
      Ast    = 20;  % подавление 20 дБ
      Nfilt  = 10;  % порядок фильтра
      
      % Создание полосового фильтра для имитации импульсного отклика с
      % ограниченной полосой
      Fd2 = fdesign.bandpass('N,Fst1,Fst2,Ast',Nfilt,Flow,Fhigh,Ast,Fs);
      Hd2 = design(Fd2,'cheby2','FilterStructure','df2tsos',...
          'SystemObject',true);
      
      % Фильтрация шума для получения импульсной характеристики
      G = step(Hd2,[zeros(delayW,1); log(0.99*rand(N-delayW,1)+0.01).*...
          sign(randn(N-delayW,1)).*exp(-0.01*(1:N-delayW)')]);
      G = G/norm(G);
      
      plot(t,G,'b');
      xlabel('Время, с');
      ylabel('Значения коэффициентов');
      title('Импульсная характеристика первичного пути распространения');
      

      Подавляемый шум.

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

      Инициализация системы.

      Самым распространенным алгоритмом для систем активного шумоподавления является LMS- алгоритм с дополнительной фильтрацией выходного сигнала фильтра перед формированием сигнала ошибки (Filtered-x LMS algorithm). Этот алгоритм использует оценку вторичного пути распространения для расчета выходного сигнала, который разрушительно влияет на нежелательный шум в области датчика измерения ошибки. Опорным сигналом является зашумленная версия нежелательного звука, измеренная вблизи его источника. Мы будем использовать управляемый фильтр с длительностью отклика около 44 мс и шагом подстройки равным 0.0001.

      % КИХ фильтр используемый для моделирования первичного пути распространения
      Hfir = dsp.FIRFilter('Numerator',G.');
      
      % Адаптивный фильтр реализующий алгоритм Filtered-X LMS
      L = 350;
      muW = 0.0001;
      Hfx = dsp.FilteredXLMSFilter('Length',L,'StepSize',muW,...
          'SecondaryPathCoefficients',Hhat);
      
      % Синтез шума с помощью синусоид
      A = [.01 .01 .02 .2 .3 .4 .3 .2 .1 .07 .02 .01]; La = length(A);
      F0 = 60; k = 1:La; F = F0*k;
      phase = rand(1,La); % случайная начальная фаза
      Hsin = dsp.SineWave('Amplitude',A,'Frequency',F,'PhaseOffset',phase,...
          'SamplesPerFrame',512,'SampleRate',Fs);
      
      % Проигрыватель аудио для воспроизведения результатов работы алгоритма
      Hpa = dsp.AudioPlayer('SampleRate',Fs,'QueueDuration',2);
      
      % Анализотор спектра
      Hsa = dsp.SpectrumAnalyzer('SampleRate',Fs,'OverlapPercent',80,...
          'SpectralAverages',20,'PlotAsTwoSidedSpectrum',false,...
          'ShowLegend',true);
      

      Симуляция разработанной системы активного шумоподавления.

      Здесь мы сымитируем работу системы активного шумоподавления. Чтобы подчеркнуть разницу первые 200 итераций шумоподавление будет отключено. Звук на микрофоне до подавления представляет характерный «вой» промышленных моторов.

      Результирующий алгоритм сходится примерно через 5 с (имитационных) после включения адаптивного фильтра. Сравнивая спектры сигнала остаточной ошибки и исходного зашумленного сигнала, можно наблюдать, что большая часть периодичных компонент была успешно подавлена. Однако эффективность стационарного шумоподавления может быть неравномерна по всем частотам. Такое часто бывает в реальных системах, применяемых для задач активной борьбы с шумом. При прослушивании сигнала ошибки раздражающий «вой» значительно снижается.

      for m = 1:400
          s = step(Hsin); % генерация синусоид со случайной фазой
          x = sum(s,2);   % генерация шума сложением всех синусоид
          d = step(Hfir,x) + ... % распространение шума через первичный канал
              0.1*randn(size(x)); % добавление шума, сопроводающего процесс измерения
          if m <= 200
              % отключение шумоподавления на первые 200 итераций
              e = d;
          else
              % включение алгоритма шумоподавления
              xhat = x + 0.1*randn(size(x));
              [y,e] = step(Hfx,xhat,d);
          end
          step(Hpa,e);     % воспроизведение сигнала на выходе
          step(Hsa,[d,e]); % спектр исходного (канал 1) и ослабленного (канал 2) сигналов
      end
      release(Hpa); % отключение динамиков
      release(Hsa); % отключение спектроанализатора
      
      Warning: The queue has underrun by 3456 samples. Try increasing queue duration,
      buffer size, or throughput rate. 
      
    • Alexey Antipin
      28.01.2014
      18:57
      Как правило, процесс сбора данных с внешнего устройства состоит из следующих четырех этапов:
      1.   Инициализация: создание объекта устройства.
      2.   Настройка: добавление каналов и управление процессом сбора данных через параметры устройства.
      3.   Выполнение: запуск устройства с последующим сбором или пересылкой данных.
      4.   Завершение: удаление объекта устройства.



      Чтобы убедиться, что основная частота камертона составляет 440 Гц, сигнал оцифровывается и передается  в MATLAB  для анализа. Это установка для примера описанного ниже.

      В этом примере мы проверим, что основная (наименьшая) частота камертона составляет 440 Гц. Для этого будем использовать микрофон и звуковую карту, чтобы записать звуковые данные. А после применим к ним операцию БПФ и найдем частоты камертона.Начнем с записи двух секунд звукового сигнала на одном из каналов звуковой карты. Так как основная частота камертона составляет 440 Гц, то можно выбрать наименьшую частоту дискретизации равную 8000 Гц.
      Затем установим камертон рядом с микрофоном, придадим ему вибрацию и запустим сбор данных. Данная сессия сбора данных описана ниже.

      Инициализация.
      Первый шагом создадим объект с аналоговым входом (AI) для звуковой карты.

      Код
       AI = analoginput('winsound');

      Настройка.
      Затем добавим к AI один канал и установим частоту дискретизации 8000 Гц и время сбора данных 2 с.

      Код
      addchannel(AI, 1);
      Fs = 8000;                     % Частота дискретизации 8000 Гц
      set (AI, 'SampleRate', Fs)
      duration = 2;                  % Сбор будет длиться 2 секунды
      set(AI, 'SamplesPerTrigger', duration*Fs;

      Сбор данных.
      Теперь мы готовы начать сбор данных. По умолчанию, триггер устройства настроен на запуск по команде start. Перед этим нужно ударить по камертону, чтобы он начал вибрировать и издавать звук (свист будет также работать J).

      Код
      start(AI);

      Чтобы списать все данные.
      Код
      data = getdata(AI);

      Завершение работы с устройством.
      Как только все данные получены, можно завершить сессию сбора данных, удалив объект AI из пространства переменных.

      Код
      delete(AI)

      Обработка результатов.
      Теперь определим частотные компоненты сигнала и отобразим результаты. Сначала рассчитаем абсолютное значение от БПФ данных.

      Код
       xfft = abs(fft(data));

      A затем преобразуем его в дБ и выделим реальные частотные компоненты.
      Код
       mag = 20*log10(xfft);
       mag = mag(1:end/2);


      Результаты показывают основную частоту около 440 Гц и первую гармонику около 880 Гц. Простейший путь найти основную частоту:
      Код
      [ymax,maxindex]=max(mag);

      Ответ: 441 Гц.

      Применение других видов оборудования.
      Этот пример можно повторить и с другими устройствами, изменив всего две строчки кода. Например, возьмем многофункциональную карту от National Instruments и создадим для нее объект аналогового ввода:

      Код
      AI=analoginput('nidaq','Dev1');
      addchannel(AI,0)

      Аналогично, если бы мы использовали плату сбора данных компании Measurement Computing, наш код выглядел бы так:
      Код
      AI=analoginput('mcc',8);
      addchannel(AI,1)
    • Alexey Antipin
      04.10.2013
      12:57
      Фильтрация сигналов – это одна из важнейших операций в системах ЦОС. Она применяется повсюду: в системах связи и радиолокации, в медицинском оборудовании и аудиосистемах… Вообще, везде, где стоит вопрос «запихнуть» сигнал из реального мира в компьютер. В одном только мобильном телефоне (или как сейчас модно говорить Smartphone) на скорую руку я насчитал их штук 30. Это и АЦП для радио и звуковых сигналов, и DDC в приемнике (тут их 4), и фильтры типа приподнятого косинуса (служат для формирования импульсов в радиоканал и из него), и фильтр Гаусса в модуляции GMSK, и схемы шумоподавления и автобалансировки уровня громкости при передаче речи, и банки фильтров в эквалайзере проигрывателя музыки и т.д. и т.п. Тут я не говорю еще о множестве фильтров, через которые сигнал со светочувствительной матрицы вашей фотокамеры проходит  прежде, чем отобразиться на дисплее, так как это относится к несколько другому типу (двухмерных) фильтров. О них не пойдет речи в этой статье, но, к слову, в MATLAB есть инструменты и для этого.
      Итак, вопрос о создании фильтров у инженеров возникает весьма часто, и сейчас мы посмотрим, как легко и непринужденно эти задачи решаются в среде MATLAB.
      .
      Загрузка плеера
  • RSS подписка
    Архив материалов
    Май 2017
    Пн Вт Ср Чт Пт Сб Вс
    1 2 3 4 5 6 7
    8 9 10 11 12 13 14
    15 16 17 18 19 20 21
    22 23 24 25 26 27 28
    29 30 31        
    Поиск