Введение

Сигналы могут рассматриваться как импульсная характеристика авторегрессивной линейной модели и, таким образом, могут быть смоделированы с использованием таких инструментов, как ar.

Данные для сигналов могут быть инкапсулированы в iddata путем установки выходных данных объекта в значения сигнала и оставления ввода пустым. Например, если x(t) представляет моделируемый сигнал, затем соответствующий iddata объект может быть создан как: data = iddata(x,[],T);, где T - время выборки x.

Стандартные средства идентификации, такие как n4sid, ssest, ar и arx может использоваться для оценки характеристик данных «только для вывода». Эти модели оцениваются по их спектральной оценочной способности, а также их способности предсказывать будущие значения сигнала из измерения их прошлых значений.

Анализ данных

Этот пример мы начинаем с загрузки данных для сигнала тока от трансформатора:

load current.mat

Теперь мы упаковываем текущие данные (i4r) в iddata объект. Время выборки составляет 0,001 секунды (1 ms).

i4r = iddata(i4r,[],0.001)  % Second argument empty for no input

i4r =

Time domain data set with 601 samples.
Sample time: 0.001 seconds             
                                       
Outputs      Unit (if specified)       
   y1                                  
                                       

Давайте проанализируем эти данные. Сначала взгляните на данные:

plot(i4r)

Закрытый вид данных показан ниже:

plot(i4r(201:250))

Затем вычисляем необработанную периодограмму сигнала:

ge = etfe(i4r)
spectrum(ge)

ge =
IDFRD model.
Contains spectrum of signal in the "SpectrumData" property.
 
Output channels: 'y1'
Status:                                        
Estimated using ETFE on time domain data "i4r".

Эта периодограмма обнаруживает несколько гармоник, но не очень гладкая. Сглаженную периодограмму получают:

ges = etfe(i4r,size(i4r,1)/4);
spectrum(ge,ges);
legend({'ge (no smoothing)','ges (with smoothing)'})

Сконфигурируйте график для использования линейной шкалы частот и единиц Гц:

h = spectrumplot(ges);
opt = getoptions(h);
opt.FreqScale = 'linear';
opt.FreqUnits = 'Hz';
setoptions(h,opt);
axis([0 500,-5 40])
grid on, legend('ges')

Мы ясно видим доминирующую частотную составляющую 50 Гц и ее гармоники.

Проведем спектральный анализ данных с использованием spa, которое использует окно Ханна для вычисления спектральных амплитуд (в отличие от etfe который просто вычисляет необработанную периодограмму). Стандартная оценка (с размером окна% по умолчанию, который не подстраивается под резонансные спектры) дает:

gs =  spa(i4r);
hold on
spectrumplot(gs);
legend({'ges (using etfe)','gs (using spa)'})
hold off

Мы видим, что потребуется очень большое окно запаздывания, чтобы увидеть все тонкие резонансы сигнала. Стандартный спектральный анализ работает плохо. Нам нужна более сложная модель, например, обеспечиваемая методами параметрического авторегрессионного моделирования.

Параметрическое моделирование текущего сигнала

Теперь вычислим спектры параметрическими AR-методами. Модели 2-го 4-го и 8-го порядка получены:

t2 = ar(i4r,2);
t4 = ar(i4r,4);
t8 = ar(i4r,8);

Рассмотрим их спектры:

spectrumplot(t2,t4,t8,ges,opt);
axis([0 500,-8 40])
legend({'t2 (2nd order AR)','t4 (4th order AR)','t8 (8th order AR)','ges (using spa)'});

Мы видим, что параметрические спектры не способны улавливать гармоники. Причина в том, что AR-модели уделяют слишком много внимания более высоким частотам, которые трудно смоделировать. (См. Ljung (1999) Пример 8.5).

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

Каким будет полезный заказ? Мы можем использовать arxstruc чтобы определить это.

V = arxstruc(i4r(1:301),i4r(302:601),(1:30)'); % Checking all order up to 30

Для интерактивного выбора наилучшего порядка выполните следующую команду: nn = selstruc(V,'log');

Как показано на рисунке выше, наблюдается резкое падение n=20. Поэтому давайте выберем этот порядок для следующих обсуждений.

t20 = ar(i4r,20);
spectrumplot(ges,t20,opt);
axis([0 500 -25 80])
legend({'ges (using spa)','t20 (20th order AR)'});

Все гармоники сейчас подобраны, но почему уровень упал? Причина в том, что t20 содержит очень тонкие, но высокие пики. С грубой сеткой частотных точек в t20 мы просто не видим истинных уровней пиков. Мы можем проиллюстрировать это следующим образом:

g20c = idfrd(t20,(551:650)/600*150*2*pi); % A frequency region around 150 Hz
spectrumplot(ges,t20,g20c,opt)
axis([0 500 -25 80])
legend({'ges (using spa)','t20 (20th order AR)','g20c (resp. around 150 Hz)'});

Как показывает этот график, модель t20 является довольно точным; при построении графика на тонкой частотной сетке он фиксирует гармоники сигнала достаточно точно.

Моделирование только гармоник нижнего порядка

Если мы в первую очередь заинтересованы в нижних гармониках и хотим использовать модели более низкого порядка, нам придется применить предварительную фильтрацию данных. Мы выбираем фильтр Баттерворта 5-го порядка с частотой отсечки 155 Гц. (Это должно охватывать режимы 50, 100 и 150 Гц):

i4rf = idfilt(i4r,5,155/500); % 500 Hz is the Nyquist frequency
t8f = ar(i4rf,8);

Теперь сравним спектр, полученный из отфильтрованных данных (модель 8-го порядка), с спектром для нефильтрованных данных (8-го порядка) и с периодограммой:

spectrumplot(t8f,t8,ges,opt)
axis([0 350 -60 80])
legend({'t8f (8th order AR, filtered data)',...
   't8 (8th order AR, unfiltered data)','ges (using spa)'});

Мы видим, что с отфильтрованными данными мы достаточно хорошо подбираем первые три пика в спектре.

Мы можем вычислить численные значения резонансов следующим образом: Корни дискретизированной синусоиды частоты, скажем, om, расположены на окружности блока при exp(i*om*T), T время выборки. Таким образом, мы действуем следующим образом:

a = t8f.a % The AR-polynomial
omT = angle(roots(a))'
freqs = omT/0.001/2/pi';
% show only the positive frequencies for clarity:
freqs1 = freqs(freqs>0) % In Hz

a =

  Columns 1 through 7

    1.0000   -5.0312   12.7031  -20.6934   23.7632  -19.6987   11.5651

  Columns 8 through 9

   -4.4222    0.8619


omT =

  Columns 1 through 7

    1.3591   -1.3591    0.9620   -0.9620    0.3146   -0.3146    0.6314

  Column 8

   -0.6314


freqs1 =

  216.3063  153.1036   50.0665  100.4967

Таким образом, мы находим первые три гармоники (50, 100 и 150 Гц) достаточно хорошо.

Мы также могли бы проверить, насколько хорошо модель t8f способен прогнозировать сигнал, скажем, 100 мс (100 шагов) вперед, и оценивать подгонку на выборках, 201 к 500:

compare(i4rf,t8f,100,compareOptions('Samples',201:500));

Как было замечено, модель первых 3 гармоник является довольно хорошей для прогнозирования будущих выходных значений, даже на 100 шагов вперед.

Моделирование только гармоник высшего порядка

Если бы нас интересовали только четвертая и пятая гармоники (около 200 и 250 Гц), мы бы продолжили полосовую фильтрацию данных в этот более высокий частотный диапазон:

i4rff = idfilt(i4r,5,[185 275]/500);
t8fhigh = ar(i4rff,8);
spectrumplot(ges,t8fhigh,opt)
axis([0 500 -60 40])
legend({'ges (using spa)','t8fhigh (8th order AR, filtered to high freq. range)'});

Таким образом, мы получили хорошую модель в t8fhigh для описания 4-й и 5-й гармоник. Таким образом, мы видим, что при правильной предварительной фильтрации могут быть построены параметрические модели низкого порядка, которые дают хорошее описание сигнала в нужных частотных диапазонах.

Заключения

Какая модель лучшая? Вообще, модель более высокого порядка давала бы более высокую верность. Чтобы проанализировать это, рассмотрим, что даст модель 20-го порядка с точки зрения ее способности оценивать гармоники:

a = t20.a  % The AR-polynomial
omT = angle(roots(a))'
freqs = omT/0.001/2/pi';
% show only the positive frequencies for clarity:
freqs1 = freqs(freqs>0) %In Hz

a =

  Columns 1 through 7

    1.0000    0.0034    0.0132    0.0012    0.0252    0.0059    0.0095

  Columns 8 through 14

    0.0038    0.0166    0.0026    0.0197   -0.0013    0.0143    0.0145

  Columns 15 through 21

    0.0021    0.0241   -0.0119    0.0150    0.0246   -0.0221   -0.9663


omT =

  Columns 1 through 7

         0    0.3146   -0.3146    0.6290   -0.6290    0.9425   -0.9425

  Columns 8 through 14

    1.2559   -1.2559    1.5726   -1.5726    1.8879   -1.8879    2.2027

  Columns 15 through 20

   -2.2027    2.5136   -2.5136    3.1416    2.8240   -2.8240


freqs1 =

  Columns 1 through 7

   50.0639  100.1139  149.9964  199.8891  250.2858  300.4738  350.5739

  Columns 8 through 10

  400.0586  500.0000  449.4611

Мы видим, что эта модель очень хорошо воспринимает гармоники. Эта модель будет предсказывать 100 шагов вперед следующим образом:

compare(i4r,t20,100,compareOptions('Samples',201:500));

Теперь у нас 93% соответствует t20, в отличие от 80% для t8f.

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

Дополнительная информация

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