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

    • 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 подписка
    Архив материалов
    Май 2019
    Пн Вт Ср Чт Пт Сб Вс
        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    
    Поиск