Введение

Сигналы могут быть рассмотрены как импульсная характеристика авторегрессивной линейной модели, и, таким образом, могут быть смоделированы с помощью таких инструментов, как 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 содержит очень тонкий но высокий peaks. С сырой сеткой частотных точек в t20 мы просто не видим истинных уровней peaks. Мы можем проиллюстрировать это следующим образом:

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 является естественным выбором, как с точки зрения захвата гармоник, так и с точки зрения его возможностей предсказания. Для моделей в определенных областях значений мы, однако, можем сделать очень хорошо с моделями более низкого порядка, но затем мы должны предфильтровать данные соответственно.

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

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