Способ идентификации мультисинусоидальных цифровых сигналов
METHOD FOR IDENTIFYING MULTISINUSOIDAL DIGITAL SIGNALS

Изобретениие относится к области обработки информации и измерительной техники и может быть использовано для контроля работоспособности электротехнических и электромеханических устройств. Способ может быть применен для выделения полезного мультисинусоидального сигнала и для определения его математической модели на основе цифровых отсчетов, и реализован с использованием ЭВМ в автоматическом режиме, в реальном масштабе времени.
Известен СПОСОБ ВЫДЕЛЕНИЯ ПОЛЕЗНОГО СИГНАЛА РЕАЛИЗУЕМОГО ПРОЦЕССА (патент РФ №2394216, G01M 7/00, G01R 23/16, опубл. 10.07.2010 http://www.freepatent.ru/patents/2394216), основанный на узкополостной фильтрации сигнала. Узкополосную фильтрацию этого сигнала путем автоматической частотной настройки соответствующего фильтра проводят с последующим детектированием отфильтрованного сигнала и его сравнением с пороговой величиной, после чего формируют сигнал управления компьютерной обработкой процесса, а затем осуществляют широкополосную фильтрацию реализуемого процесса.
Данный способ имеет следующие недостатки:
- о спектральном составе сигнала можно судить только с помощью визуального анализа графиков, что не позволяет произвести автоматизацию расчетов и снижает быстродействие и точность метода;
- после фильтрации шумовой составляющей определяются только пороговые значения периодического сигнала, а число гармоник и их параметры (амплитуда, частота) остаются неизвестными.
Известен способ определение наличия гармонических составляющих и их частот в дискретных сигналах (патент РФ №2498324, МПК G01R 23/00, опубл. 10.11.2013 http://www.freepatent.ru/patents/2498324), использующий быстрое преобразование Фурье входных сигналов. Определяют комплексно-сопряженные значения прямого преобразования Фурье анализируемого дискретного сигнала, попарно умножают полученные комплексные сигналы прямого преобразования Фурье анализируемого дискретного сигнала, тем самым формируют m сигналов и подвергают их обратному преобразованию Фурье, определяют частотно-временную автокорреляционную функцию и строят ее график, по которому судят о наличии гармонических составляющих и их частотах. Данный способ авторы предлагают использовать для анализа сигналов различного происхождения в задачах неразрушающего контроля и диагностики оборудования, однако такие сигналы имеют в своем составе шумовые составляющие. Данный способ не предусматривает предварительную фильтрацию шумов, что может впоследствии привести к ошибочным оценкам структуры и параметров цифровых сигналов. Кроме того, данный способ имеет следующие недостатки, снижающие быстродействие и точность метода:
- большое количество операций для осуществления быстрого преобразования Фурье, формирование большого количества дополнительных сигналов (более 1000), вследствие чего наблюдается громоздкость расчетов и неточность вычислений;
-4
- для реализации метода необходимо большое количество измерений сигнала (более 16000) и малый шаг дискретизации (менее 10);
- визуальный анализ графиков частотно-временной корреляционной функции для определения количества гармоник и частот сигнала, не позволяющий произвести автоматизацию расчетов;
- нет возможности оценить амплитуды гармонических составляющих.
Известен способ гармонического анализа периодического многочастотного сигнала (патент РФ №2435168, МПК G01R 23/16, опубл. 27.11.2011, http://www.freepatent.ru/patents/2435168), использующий итерационный процесс определения необходимой формы весовой функции. В результате многократных оценок частот периодического сигнала форма весовой функции подбирается такой, чтобы на частоте любого из анализируемых составляющих сигнала слагаемые спектра других составляющих сигнала по амплитуде были равны нулю, причем кратность нуля, определяемая порядком производных модуля спектра, может быть заданной степени. Амплитуду и фазу каждого составляющего сигнала оценивают после определения точных значений частот.
Данный способ имеет следующие недостатки, снижающие быстродействие и точность метода:
- итерационная процедура подбора весовой функции и проведение дополнительных циклов гармонического анализа требуют значительного числа операций, поэтому обладают недостаточной быстродействием и наличием погрешности расчетов;
- громоздкость и неточность нахождения Фурье-образа спектра;
- визуальный анализ лепестков спектра значительно затрудняет процедуру оценки частот сигнала и не позволяет автоматизировать расчеты;
- для реализации метода необходимо большое количество измерений сигнала (более 100 периодов колебаний сигнала).
Наиболее близким к предлагаемому способу является способ идентификации мультисинусоидальных цифровых сигналов (патент РФ 2703933, МПК G01R 23/16, опубл. 22.10.2019), включающий использование непрерывных цепных дробей путем измерения сигнала в равноотстоящие промежутки времени, подачи их на идентификатор непрерывной С-дроби с последующим восстановлением числа гармоник, круговых частот, амплитуд, прогнозирующей модели и модельных значений сигнала, отличающийся тем, что при обработке значений сигнала в идентификаторе непрерывной С-дроби на основе правила останова судят о наличии смещения сигнала, а при его выявлении выполняют последующее восстановление его значения и далее с учетом полученных результатов судят о начальных фазовых сдвигах в гармониках сигнала.
Этот способ может применяться только для не зашумленных мультисинусоидальных сигналов, что значительно снижает качество идентификации, поскольку в реальных системах контроля работоспособности электротехнических и электромеханических устройств измерительная информация обычно искажена помехами. Для мультисинусоидальных сигналов, содержащих шумовые составляющие, способ приводит к получению ошибочной модели полезного сигнала.
Предлагаемым изобретением ставится задача идентификации мультисинусоидального цифрового сигнала, позволяющая осуществить фильтрационную обработку сигнала с помощью структурного преобразователя, в идентификаторе непрерывной цепной С-дроби вести подсчет числа гармонических составляющих сигнала, оценивать параметры сигнала (амплитуда, частота) поэтапно по мере убывания значений их амплитуд в восстановителе модели гармоники с последующим восстановлением остатка сигнала, определять наличие смещения и восстанавливать его значение, расширяя тем самым область применения метода идентификации для различных типов мультисинусоидальных сигналов, автоматически определяя структуру и параметры полезного сигнала, значительным образом повышая степень автоматизации, точность, достоверность и быстродействие процесса идентификации, что позволяет использовать данный метод как способ текущей идентификации в реальном масштабе времени, определять изменение структуры модели и, тем самым, повышать точность, надежность и качество результатов моделирования.
Техническим результатом предлагаемого способа идентификации мультисинусоидальных цифровых сигналов является то, что обеспечивается быстродействие, универсальность реализации способа, простота и высокая точность вычислений в условиях наличия помех в сигналах.
Изобретение поясняется на фигурах 1-8.
м
Структурная схема системы, реализующая предлагаемый способ, изображена на фиг. 1. Она содержит датчик анализируемого сигнала (ДАС), с выхода которого сигнал x(kΔt) поступает на вход последовательно соединенного блока 1 - структурного преобразователя, в котором осуществляется предварительная фильтрационная обработка значений сигнала на основе трех последовательных структурных преобразований. Результаты третьего структурного преобразования сигнала последовательно передаются на вход блока 2 - идентификатора непрерывной цепной С-дроби, который включает в себя 2 подблока. В подблоке 2.1 счетчика гармоник рассчитывается матрица-идентификатор, в которой 4-ая строка будет являться нулевой, и далее запускается подсчет числа гармоник i=i+1 (первоначально i=0). В подблоке 2.2 аппроксиматора модели определяется модель гармоники в форме дискретной передаточной функции (ДПФ) объекта-идентификатора. Далее в блоке 3 - восстановителя модели гармоники определяют параметры гармоники - круговую частоту и амплитуду. Полученную модель гармоники подают на вход блока 4 - восстановителя сигнала остатка, в котором из исходного сигнала вычитают модельные значения сигнала выявленной гармоники, и далее результаты сигнала остатка подаются на вход блока 1. Поскольку число гармоник исходного сигнала неизвестно, все действия повторяют до выполнения правила останова в блоке 2 - появления 4-ой ненулевой строки в матрице-идентификаторе. Выявление гармоник в блоке 2 осуществляется последовательно по мере убывания значений их амплитуд. Результат выполнения правила останова поступает на вход блока 5 - анализатора смещения, в котором определяют наличие смещения и выполняют восстановление его значения. Далее в блоке 6 - восстановителя модели сигнала определяют прогнозирующую модель мультисинусоидального сигнала. Затем прогнозирующая модель поступает на вход блока 7 - восстановителя модельных значений сигнала, в котором определяется модельный сигнал x(kΔt).
Предлагаемый способ осуществляется следующим образом: с ДАС результаты измерений сигнала x(kΔt) (k=0, 1, 2, …) в равноотстоящие промежутки времени с шагом дискретизации Δt поступают на вход блока 1, в котором производится расчет первого структурного преобразования
второго структурного преобразования
и третьего структурного преобразования
1
2
3
где N - число измерений сигнала x(kΔt); N, N - число расчетных значений второго и третьего структурных преобразований соответственно.
Результаты расчетов значений третьего структурного преобразования (3) передают на вход блока 2, где заносят их в 0-ую строку матрицы-идентификатора:
Последующие строки матрицы-идентификатора определяются с помощью формулы
где , , m=1, 2, 3, …, n=0, 1, 2, …
Если анализируемый сигнал является мультисинусоидальным и представим в виде модели
i
i
где μ - постоянное смещение, n - число гармоник, ω - круговая частота i-ой гармоники, C - амплитуда i-ой гармоники, a(kΔt) - шумовая составляющая, k=0, 1, 2, …, то в матрице (4) 4-ая строка будет являться нулевой. Запускается подсчет числа гармоник i=i+1 (первоначально i=0).
Первый столбец матрицы-идентификатора позволяет построить непрерывную цепную С-дробь в виде:
1
сворачивая которую, определяют модель гармоники в форме ДПФ объекта-идентификатора. Полученную модель передают на вход блока 3, где приступают к нахождению круговой частоты ω данной гармоники, имеющей максимальную амплитуду:
1
где z=u±iv - полюс ДПФ (7). Тогда максимальную амплитуду гармоники находят по формуле:
3
где k=1, …, N. Выявленная гармоника с максимальной амплитудой будет иметь вид:
1
Поскольку число гармоник исходного сигнала неизвестно, необходимо провести следующую процедуру. Полученную модель гармоники (10) передают на вход блока 4, где значения (10) вычитаются из исходного сигнала, получая на выходе блока сигнал остатка x(kΔt)
После этого с сигналом (11) следует повторить всю процедуру для выявления остальных гармоник в нем. Результаты сигнала остатка (11) поступают на вход блока 1. Следующей гармоникой будет являться та, у которой значение амплитуды является наибольшим из оставшихся. Процедура повторяется до тех пор, пока в матрице-идентификаторе (4) элементы 4-й строки будут обращаться в нуль. Таким образом, критерием останова процедуры выявления гармонических компонент в исходном сигнале служит появление ненулевых элементов в 4-й строке матрицы-идентификатора.
После выполнения правила останова модель сигнала будет иметь вид
Она подается на вход блока 5, где определяют значение смещения
Если μ≠0, то фиксируют наличие смещения и принимают его равным (13). Результаты поступают на вход блока 6, где окончательно формируют прогнозирующую модель мультисинусоидального сигнала в виде (14)
если μ≠0, или в виде (12), если μ=0.
Идентифицированная модель мультисинусоидального сигнала поступает на вход блока 7, который восстанавливает значения модельного сигнала
Таким образом, предлагаемый способ идентификации мультисинусоидальных цифровых сигналов, отличается от известного тем, что восстановление гармоник в идентификаторе непрерывной дроби осуществляется последовательно по мере убывания значений их амплитуд до выполнения правила останова с предварительной фильтрационной обработкой значений сигнала, а после выполнения правила останова судят о наличии смещения и выполняют восстановление его значения.
Способ позволяет структурно преобразовать сигнал с целью фильтрации помех, автоматически оценить неизвестную заранее структуру математической модели мультисинусоидального сигнала, а также все его параметры. Вследствие использования метода процесс идентификации сигнала упрощается, увеличивается быстродействие, повышается точность, качество и достоверность результатов исследований. Универсальность и работоспособность метода продемонстрирована на примере различного рода мультисинусоидальных сигналов. В конечном итоге, предлагаемый способ позволяет достоверно прогнозировать значения анализируемых физических процессов различных типов, принимать адекватные решения по их контролю, управлению и диагностике.
Пример 1
График выходного сигнала спирального дозатора при подаче на вход йодида калия приведен на фиг. 2. К закономерным колебаниям выходного сигнала дозатора могут добавляться особые условия (зависание материала в бункере, очаговые уплотнения материала, влажность и т.п.), приводящие к стохастическим изменениям выходного сигнала дозатора. Пусть модель выходного сигнала имеет вид
1
С выхода ДАС поступают измерения анализируемого сигнала (шаг дискретизации Δt=0.5, число измерений N=240), в блоке 1 над сигналом берутся три последовательных структурных преобразования. График третьего структурного преобразования приведен на фиг. 3.
Результаты расчетов значений третьего структурного преобразования (3) поступают на вход блока 2 и заполняют 0-ую строку матрицы-идентификатора:
Ее 4-ая строка обращается в нулевую. В подблоке 2.2 находят ДПФ объекта-идентификатора
1
1
1
Полученную модель передают на вход блока 3, где приступают к нахождению круговой частоты ω гармоники с максимальной амплитудой по формуле (8) ω=2.00330. Сама амплитуда согласно (9) получается равна C=0.49973.
1
1
На вход блока 4 подают полученную модель гармоники s(kΔt)=0.49973sin(2.00330⋅kΔt), получая на выходе блока 4 сигнал остатка x(kΔt)=x(kΔt)-0.49973sin(2.00330⋅kΔt). Измерения сигнала остатка поступают на вход блока 1, где с сигналом осуществляются три последовательных структурных преобразования. Результаты расчетов значений третьего структурного преобразования (3) поступают на вход блока 2 в 0-ую строку матрицы-идентификатора:
В 4-й строке матрицы-идентификатора появляются ненулевые элементы, что служит критерием останова процедуры выявления гармонических компонент в сигнале. Тогда модель сигнала будет иметь вид
Она подается на вход блока 5, где фиксируют наличие смещения и определяют его значение μ=8.02655 согласно (13). Результаты поступают на вход блока 6, где окончательно формируют прогнозирующую модель сигнала в виде
-2
Максимальная абсолютная погрешность оценки параметров имеет порядок 10, структура сигнала (15) восстановлена верно. Идентифицированная модель сигнала поступает на вход блока 7, который восстанавливает значения модельного сигнала
Результаты расчетов модельных значений сигнала в блоке 7 приведены в Таблице 1, где  - модельные значения сигнала, рассчитанные с помощью предлагаемого способа,  - погрешности модельных значений. Сравнивая экспериментальные и модельные значения сигнала на выходе спирального дозатора и принимая во внимание наличие стохастических помех и погрешностей, можно сделать заключение о высокой точности идентификации выходного сигнала (с точностью восстановления до вычислительных погрешностей).
Пример 2.
Выходной вибросигнал роторного оборудования горнодобывающего производства имеет вид (график сигнала приведен на фиг. 4):
1
С выхода ДАС поступают измерения анализируемого сигнала (шаг дискретизации Δt=0.5, число измерений N=100), в блоке 1 над сигналом берутся три последовательных структурных преобразования. Далее рассчитывается матрица-идентификатор, в которой 0-ая строка заполняется значениями третьего структурного преобразования (график третьего структурного преобразования приведен на фиг. 5), поступающих с выхода блока 1:
Ее 4-ая строка содержит нулевые значения. В подблоке 2.2 находят ДПФ объекта-идентификатора
1
1
Далее модель ДПФ объекта-идентификатора передают на вход блока 3, где приступают к нахождению параметров гармоники с максимальной амплитудой: круговая частота ω=4.002, амплитуда C=7.778.
На вход блока 4 подают полученную модель гармоники  получая на выходе блока 4 сигнал остатка  Сигнала остатка поступает на вход блока 1, где с сигналом осуществляются три последовательных структурных преобразования. Результаты расчетов значений третьего структурного преобразования (фиг. 6) поступают на вход блока 2 и заполняется матрица-идентификатор:
В подблоке 2.2 определяют ДПФ объекта-идентификатора
2
2
Далее модель ДПФ объекта-идентификатора передают на вход блока 3, где находят параметры второй гармоники: круговая частота ω=1.411, амплитуда С=1.888.
2
2
1
На вход блока 4 передают модель второй гармоники s(kΔt)=1.888sin(1.411⋅kΔt), получая на выходе сигнал остатка x(kΔt)=x(kΔt)-1.888sin(1.411⋅kΔt), который предается на вход блока 1. Результаты расчетов третьего структурного преобразования (фиг. 7) поступают на вход блока 2 и заполняется матрица-идентификатор:
В подблоке 2.1 счетчик гармоник фиксирует наличие третьей гармоники и в подблоке 2.2 восстанавливают ее модель в форме ДПФ объекта-идентификатора
3
3
ДПФ объекта-идентификатора передают на вход блока 3, где находят параметры третьей гармоники: круговая частота ω=2.017, амплитуда С=1.026.
В блок 4 поступает модель третьей гармоники  получая на выходе сигнал остатка  Сигнал остатка поступает на вход блока 1, где производится расчет третьего структурного преобразования (фиг. 8). В блоке 2 заполняется матрица-идентификатор:
В 4-й строке матрицы-идентификатора появляются ненулевые элементы, что служит критерием останова процедуры выявления гармонических компонент в сигнале. Тогда модель сигнала будет иметь три гармоники:
Она подается на вход блока 5, где фиксируют наличие смещения и определяют его значение μ=10.047 согласно (13). Результаты поступают на вход блока 6, где окончательно формируют прогнозирующую модель мультисинусоидального сигнала
-2
Максимальная абсолютная погрешность оценки параметров имеет порядок 10, структура модели мультисинусоидального сигнала (16) восстановлена верно. Идентифицированная модель сигнала поступает на вход блока 7, который восстанавливает значения модельного сигнала  Результаты расчетов модельных значений сигнала в блоке 7 приведены в Таблице 2, где  - модельные значения сигнала, рассчитанные с помощью предлагаемого способа,  - погрешности модельных значений. Сравнивая экспериментальные и модельные значения вибросигнала на выходе роторного оборудования горнодобывающего производства, и принимая во внимание наличие стохастических помех и погрешностей, можно сделать заключение о высокой точности идентификации мультисиносидального выходного сигнала (с точностью восстановления до вычислительных погрешностей).