СПОСОБ АНАЛИЗА МНОГОЧАСТОТНЫХ СИГНАЛОВ, СОДЕРЖАЩИХ СКРЫТЫЕ ПЕРИОДИЧНОСТИ
METHOD OF ANALYSING MULTIPLE FREQUENCY SIGNALS CONTAINING HIDDEN PERIODICITY

Изобретение относится к области обработки информации и измерительной техники и может быть использовано для контроля работоспособности электротехнических и электромеханических устройств. Способ может быть применен для определения математической модели детерминированного сигнала, обладающего определенной периодичностью, на основе дискретной информации о нем, и реализован с использованием ЭВМ в автоматическом режиме, в реальном масштабе времени.
BAX
BAX
Известен способ спектрального анализа периодических многочастотных сигналов (патент РФ №2335778, МПК G01R 23/16, опубл. 2008.10.10), основанный на определении мгновенной спектральной плотности для последовательности частот. Затем определяют экстремумы характеристики распределения мгновенной спектральной плотности, по которым определяют частоты и далее амплитуды, для определения фазы формируют опорный синусоидальный сигнал, строят вольт-амперную характеристику для исходного сигнала, многократно сдвигают ее по фазе, определяя площадь вольт-амперной характеристики F, фазу каждой частотной составляющей сигнала находят из условия F=0. Далее по полученным значениям амплитуд, круговых частот и фаз судят о спектральном составе исходного сигнала.
Данный способ имеет следующие недостатки, снижающие быстродействие и точность метода:
- громоздкость и неточность метода Фурье при расчете спектральной плотности;
- необходимость процедуры перебора фаз и частот при оценке параметров сигнала;
-4
- для реализации метода необходимы большое количество измерений сигнала (порядка 10000) и малый шаг дискретизации (порядка 10).
Известен способ спектрального анализа сигналов (патент РФ №2229139, МПК G01R 23/16, опубл. 2004.05.20), основанный на перемножении опорного бинарного зондирующего и анализируемого сигналов, суммировании полученных значений за период анализируемого сигнала и расчете постоянной составляющей произведения опорного и анализируемого сигналов на каждой частоте опорного сигнала при переборе фазы опорного сигнала от 0 до 180°. По максимуму постоянной составляющей произведения сигналов определяют частоту и фазу гармонической составляющей анализируемого сигнала, а также рассчитывают ее амплитуду.
Данный способ имеет следующие недостатки, снижающие быстродействие и точность метода:
- громоздкость и неточность метода Фурье при расчете спектральной плотности;
- процедура перебора фазы опорного сигнала;
-3
- для реализации метода необходимы большое количество измерений сигнала (более 200) и малый шаг дискретизации (меньше либо равен 10).
И
Известен способ спектрального анализа сигналов (патент РФ №2229140, МПК G01R 23/16, опубл. 2004.05.20), основанный на перемножении анализируемого сигнала и опорных синусоидальных и косинусоидальных сигналов. Суммируют полученные значения на интервале Т анализируемого сигнала и рассчитывают мгновенную спектральную плотность на каждой частоте. Затем рассчитывают амплитудное значение и фазовый угол каждой гармонической составляющей.
Данный способ имеет следующие недостатки, снижающие быстродействие и точность метода:
- громоздкость и неточность метода Фурье при расчете спектральной плотности;
-3
- для реализации метода необходимо большое количество измерений сигнала (более 200) и малый шаг дискретизации (меньше либо равен 10);
- необходимость перебора круговых частот опорного сигнала, производящаяся до появления соответствия с частотами анализируемого сигнала.
Наиболее близким к предлагаемому способу является способ выявления срытых периодичностей сигнала на основе непрерывных цепных дробей (Серебренников М.Г., Первозванский А.А. Выявление скрытых периодичностей. - М.: Наука, 1965), сущность которого состоит в следующем. Последовательность результатов измерений сигнала в равноотстоящие промежутки времени делят пополам
1
2
1
Далее записывают эту последовательность в обратном порядке, причем числа, стоящие в первых λ столбцах, складывают, а стоящие в следующих λ столбцах - вычитают. В результате получаются коэффициенты непрерывных цепных дробей, используемые для построения рекуррентных рядов S и S, причем коэффициенты первой серии должны быть взяты справа налево, а второй серии - слева направо. Ряд S
2
делят на p, в результате имеют ряд S. После этого делят 1-у на ряд S, чтобы образовалось частное вида 1+у+qy и остаток вида
2
Если коэффициенты первого остатка не малы, то его необходимо разделить на p'y. Полученный ряд обозначают S'. Далее, деля S на S' так, чтобы частное имело вид 1+у+q'y, получают второй остаток. Если остаток равен нулю, то действия считаются законченными, в противном случае - процесс деления продолжается (число гармоник процесса n равно числу делений). В результате получаются последовательности величин р, р', р",… и q, q', q",….
2
1
Поступая с рядом S так же, как указано для S, следует в качестве первого делимого взять 1+у. В результате получаются последовательности величин (p), (p'), (p"),… и (q), (q'), (q"),…
Преобразуя р, р', р",…, q, q', q",…, (р), (р'), (р"),… и (q), (q'), (q"),… в зависимости от числа выявленных гармоник, находят модель сигнала в виде скрытых периодичностей:
i
i
1
2
где A - амплитуда i-й гармоники, w - круговая частота i-й гармоники. Причем для нахождения круговых частот достаточно воспользоваться только рядом S или S, а для получения амплитуд необходимо пользоваться двумя рядами.
Этот способ имеет следующие недостатки:
- итерационная процедура определения модели сигнала на основе перебора непрерывных цепных дробей требует выполнения значительного числа операций, что снижает быстродействие и точность расчетов;
- наличие областей неопределенности при нахождении параметров гармоник приводит к невозможности получения модели сигнала.
Предлагаемым изобретением ставится задача выявления скрытых периодичностей многочастотного сигнала, позволяющая автоматически определить структуру и неизвестные параметры математической модели сигнала, исключая итерационную процедуру корректировки модели сигнала на основе перебора непрерывных цепных дробей, значительным образом упрощая и ускоряя процесс получения модели сигнала, что дает возможность использовать данный метод в реальном масштабе времени, контролировать изменения структуры и параметров модели и тем самым повышать точность и достоверность результатов моделирования.
Предлагаемый способ выявления скрытых периодичностей обладает рядом преимуществ, которые выражаются в том, что обеспечивается быстродействие, универсальность реализации способа, простота и высокая точность вычислений.
Способ анализа многочастотных сигналов, содержащих скрытые периодичности, с использованием непрерывных цепных дробей путем измерения сигнала в равноотстоящие промежутки времени, отличается тем, что сигнал подают с датчика анализируемых сигналов в идентификатор непрерывной цепной С-дроби, в котором последовательно проводят обработку значений сигнала по формуле
-1
0
до выполнения правила останова, где α(n)=δ(n) - дельта функция Дирака, α(n)=x(n) - измерения сигнала, m=1, 2, 3,…, n=0, 1, 2,… с последующим восстановлением прогнозирующей модели сигнала в форме скрытой периодичности.
Изобретение поясняется на фигурах 1-4.
Структурная схема системы, изображенная на фиг.1 и реализующая предлагаемый способ, содержит датчик анализируемого сигнала (ДАС), к которому последовательно подсоединены блок 1 - идентификатор непрерывной цепной С-дроби, блок 2 - восстановитель модели сигнала, блок 3 - восстановитель модельных значений сигнала.
м
С выхода ДАС анализируемый сигнал x(k) поступает на вход блока 1 - идентификатора непрерывной цепной С-дроби. В блоке 1 рассчитывается идентифицирующая матрица (5), то есть производится последовательная обработка значений сигнала с помощью формулы (6) до выполнения правила останова, строится непрерывная цепная С-дробь и определяется модель сигнала в форме дискретной передаточной функции (ДПФ) формирующего объекта. Далее в блоке 2 -восстановителя модели сигнала - определяют параметры сигнала (круговые частоты, амплитуды) и его прогнозирующую модель, по которой судят о наличии скрытых периодичностей сигнала. Затем прогнозирующая модель поступает на вход блока 3 -восстановителя модельных значений сигнала, в котором определяется модельный сигнал х (k).
Предлагаемый способ осуществляется следующим образом: с ДАС результаты измерений сигнала в равноотстоящие промежутки времени с шагом дискретизации Δt поступают на вход блока 1, где рассчитывается идентифицирующая матрица:
,
m
где элементы α(nΔt) последовательно определяются с помощью формулы
-1
0
причем α(nΔt)=δ(nΔt), α(nΔt)=x(nΔt) являются начальными условиями при построении матрицы, где δ(t) - дельта функция Дирака, m=1, 2, 3,…, n=0, 1, 2,….
Элементы первого столбца идентифицирующей матрицы (5) порождают непрерывную цепную С-дробь
сворачивая которую, определяют модель формирующего объекта в форме ДПФ.
При аппроксимации дробно-рациональной функции в матрице (5) наблюдается появление нулевой строки, номер которой позволяет определить число периодических компонент. А именно, если в идентифицирующей матрице j-я строка является нулевой, то число периодических компонент в сигнале равно .
i
i
i
Если в некоторой i-той строке (i=0, 1, 2,…) матрицы (5) конечное число k первых элементов равно нулю, а последующие элементы отличны от нуля, то необходимо осуществить сдвиг влево на k элементов до появления в нулевом столбце ненулевого элемента и далее продолжить определение других элементов матрицы (5) по формуле (6). Для i-той строки при восстановлении непрерывной С-дроби (7) элемент α (0) умножается на .
i
i
Полученная ДПФ (7) поступает на вход блока 2, в котором определяются параметры гармоник - круговые частоты и амплитуды. Для этого определяются полюса ДПФ z. Если ДПФ содержит только комплексные полюса, то сигнал является периодическим или почти периодическим. В случае наличия комплексных полюсов приступают к нахождению круговых частот w из выражения
i
где z=u+iv - полюса ДПФ.
Амплитуды находятся как решение системы из n уравнений
На вход блока 3 прогнозирующая модель сигнала поступает в виде скрытых периодичностей:
Таким образом, предлагаемый способ анализа многочастотных сигналов, содержащих скрытые периодичности, отличается от известного тем, что используют последовательную процедуру выполнения операций и математическую формулу вида (6), которые позволяют автоматически определять количество и параметры скрытых периодичностей сигнала, исключая итерационную процедуру определения математической модели сигнала на основе перебора непрерывных цепных дробей. Способ позволяет выявлять наличие скрытых периодичностей для большего класса сигналов (периодических и почти периодических) на основе минимального количества наблюдений. Предлагаемый способ приводит к существенному упрощению и ускорению процесса выявления скрытых периодичностей за счет исключения большого объема вычислительных операций, что позволяет в конечном итоге достоверно прогнозировать значения физического процесса, принимать адекватные решения по его контролю, управлению и диагностике.
Пример 1.
Вибрация многомоторного винтового самолета с несинхронизированными двигателями описывается с помощью модели почти периодического сигнала
.
С ДАС на вход блока 1 - идентификатора непрерывной цепной С-дроби - поступают значения сигнала с шагом дискретизации Δt=0.1, k=0,1,2,…. График сигнала приведен на фиг.2. В блоке 1 измерения преобразовываются в непрерывную дробь путем расчета идентифицирующей матрицы на основе формулы (6)
Восьмая строка нулевая. Число периодических компонент в сигнале равно 8/4=2. Модель сигнала в форме ДПФ на выходе блока 1 имеет вид
.
Блок 2 - восстановитель модели определяет ее полюса, круговые частоты и амплитуды:
,
1
2
w=9.4248, w=4.4429,
1
A=2,
2
A=1.5.
На выходе блока 2 снимают прогнозирующую модель, содержащую две скрытых периодичности:
м
х(kΔt)=2sin(9.4348 kΔt)+1.5sin(4.4429 kΔt)
м
м
м
и передают ее на вход блока 3 - восстановителя модельных значений сигнала, на выходе которого снимают модельный сигнал х(k). Результаты вычислений приведены в Табл.1, где х(kΔt) - модельные значения сигнала, рассчитанные с помощью предлагаемого способа, e(kΔt)=x(kΔt)-х(kΔt) - погрешности модельных значений.
Таким образом, по предлагаемому способу точно восстановлена прогнозирующая модель сигнала, что в конечном итоге позволяет получить наилучший прогноз значений сигнала вибрации многомоторного винтового самолета и определить его скрытые периодичности.
Пример 2.
Измерения многочастотного сигнала напряжения на выходе электрического генератора с шагом Δt=0.001 (график сигнала приведен на фиг.3)
x(kΔt)=4sin(300·π·kΔt)+2sin(200·π·kΔt)+sin(380·π·kΔt) поступают с ДАС на вход блока 1, в котором рассчитывается идентифицирующая матрица. Двенадцатая строка в матрице нулевая, число периодических компонент сигнала равно 3. Модель формирующего объекта в форме ДПФ снимают на выходе блока 1
и передают на вход блока 2 - восстановителя прогнозирующей модели сигнала, в котором определяются круговые частоты
1
2
3
w=1193.804, w=942.477, w=628.319
и амплитуды как решение системы уравнений:
На выходе блока 2 прогнозирующая модель с 3-мя скрытыми гармониками имеет вид:
м
х(kΔt)=sin(1193.804kΔt)+4sin(942.477kΔt)+1.999sin(628.319kΔt).
м
м
Результаты вычислений модельного сигнала напряжения на выходе блока 3 приведены в Табл.2, где х(kΔt) - модельные значения сигнала, рассчитанные с помощью предлагаемого способа, e(kΔt)=x(kΔt)-х(kΔt) - погрешности модельных значений. Сравнение экспериментальных и модельных значений сигнала напряжения позволяет сделать заключение о точном (с точностью до вычислительных погрешностей) модельном его восстановлении.
Пример 3.
Сравнительный анализ предлагаемого способа с прототипом.
Имеются измерения анализируемого сигнала по прототипу, например электрического тока в цепи питания
0
4
1
с шагом дискретизации Δt=1, k=0,…,9 (график сигнала приведен на фиг.4). Были сняты десять наблюдений сигнала k=0,…9 (таблица 3, 1-й столбец). В качестве x выбирается величина, равная x=-657.79913. Далее записывается последовательность из тех же чисел, но в обратном порядке. Затем первые пять чисел складывались, а следующие пять - вычитались. В результате получились коэффициенты непрерывных цепных дробей, значения которых приведены в 2-м столбце Табл.3, причем коэффициенты первой серии должны быть взяты в обратном порядке, а второй серии - в порядке следования. Тогда ряд S имеет вид
1
2
3
S=-940.8425+1196.2099у - 517.1016у+1038.4708у+…
1
В результате деления S на р=-940.8425 получается ряд
2
3
S=1-1.2714у+0.5496у-1.1038у+…
Деля 1-у на ряд S, рассчитывается частное
2
1+у+0.2714у,
так что q=0.2714, а первый остаток от деления равен
2
3
4
-1.2045у+2.2260у-0.7411у+…,
откуда р'=-1.2045. Так как коэффициенты первого остатка не малы, то следует продолжать действие дальше.
2
Производится деление всех членов первого остатка на р'у
2
3
S'=1-1.8480у+0.6152у-0.8057у+…
2
2
Деля S на S' так, чтобы частное имело вид 1+у+q'у, рассчитывается второе частное 1+у+0.5766у, а второй остаток
3
4
5
1.1952у-0.0673у+0.5706у+…
Коэффициенты остатка не равны нулю, процесс деления должен продолжаться. Однако, согласно алгоритму, число скрытых периодичностей процесса n равно числу делений, следовательно, количество данных не достаточно для выявления скрытых периодичностей.
0
9
1
Были сняты двадцать наблюдений сигнала N=0,…19 (Табл.3, 3-й столбец). В качестве x выбирается величина, равная x=-462. Рассчитываются коэффициенты непрерывных цепных дробей первой и второй серии (4-й столбец Табл.3). Используя коэффициенты первой серии, составляется ряд S, который делится на р=-1001.1864, и получается ряд
2
3
S=1-1.1597у+0.3431у-1.0350у+…
Деля 1-у на ряд S, рассчитывается частное
2
1+у+0.1597у,
так что q=0.1597, а первый остаток от деления равен
2
3
4
-1.1579у+2.1399у-0.7124у+…,
откуда р'=-1.1579.
2
Так как коэффициенты первого остатка не малы, то следует продолжать действие дальше. Производится деление всех членов первого остатка на р'у:
2
3
S'=1-1.8480у+0.6152у+0.3895у+…
2
2
Деля S на S', получается второе частное 1+у+0.6884у, второй остаток равен нулю. Для ряда S используются коэффициенты второй серии:
2
2
3
S=-77.1864+714.7438у+116.6291у-791.0643у+…
Так как р=-77.1864, то
2
3
S=1-9.2600у-1.5110у+10.2488y+…
Результатом деления 1+у на ряд S является частное
2
1+у+10.2600у,
так что q=10.2600, а первый остаток от деления равен
2
3
4
95.5181у+14.5141у-103.2409у+…,
откуда р'=95.5181. Так как коэффициенты первого остатка не малы, то следует продолжать действие дальше.
2
Производится деление всех членов первого остатка на р'у:
2
3
S'=1+0.1519у-1.0808у-0.0761у+…
2
Результатом деления S на S' является второе частное 1+у -9.4119у, второй остаток равен нулю.
Таким образом, рассчитаны следующие значения величин p, q, p', q':
Для ряда S<sub>1</sub>
Для ряда S<sub>2</sub>
р==-1001.1864
p'=-1.1579
р=-77.1864
p'=95.5181
q=0.1597
q'=0.6884
q=10.2600
q'=-9.4119
Так как сигнал содержит две гармоники, то круговые частоты находятся по формулам
, ,
где . Параметры сигнала равны:
Для ряда S<sub>1</sub>
Для ряда S<sub>2</sub>
k<sub>1</sub>=0.6840
k<sub>2</sub>=-1.5321
k<sub>1</sub>=12.7643
k<sub>2</sub>=-13.6124
w<sub>1</sub>=1.2217
w<sub>2</sub>=2.4435
1
2
1,2
Для ряда S рассчитанные круговые частоты совпадают с истинными значениями. Для ряда S найти круговые частоты не представляется возможным, так как не выполняется условие -2<k<2. Полученные значения параметров попадают в область неопределенности, следовательно, найти значения амплитуд гармоник невозможно. Получить модель сигнала не представляется возможным.
Совокупность действий для получения модели сигнала с помощью предлагаемого способа.
С выхода ДАС поступают измерения анализируемого сигнала (Δt=1). В блоке 1 рассчитывается на основе (6) идентифицирующая матрица:
На вход блока 2 подается ДПФ:
.
1
2
1
2
Значения параметров ДПФ обрабатываются в блоке 2 и определяются круговые частоты w=2.4435, w=1.2217 и амплитуды A=593, А=462. На выходе блока 2 снимают прогнозирующую модель сигнала
м
х(kΔt)=593sin(2.4435kΔt)+462sin(1.2217kΔt),
м
м
содержащую две скрытых периодичности, и подают ее на вход блока 3. Результаты расчетов модельных значений сигнала в блоке 3 приведены в Табл.4, где х(kΔt) -модельные значения сигнала, рассчитанные с помощью предлагаемого способа, e(kΔt)=x(kΔt)-х(kΔt) - погрешности модельных значений. Сравнение экспериментальных и модельных значений сигнала электрического тока позволяет сделать заключение о точном (с точностью до вычислительных погрешностей) модельном его восстановлении.