Ідентифікація

 

 

1.2 Побудова моделі за допомогою поліноміальної інтерполяції поліномом 7-го степеня

Для апроксимації сигналів по відомим вхідним та вихідним параметрам використовують метод найменших  квадратів. В даному випадку рівняння апроксимуючої функції буде:

де  – невідомі коефіцієнти.

Запишемо умову метода найменших квадратів:

тобто

 

Продиференціювавши даний вираз по і прирівнявши його до нуля розв’язуємо слідкуючу систему рівнянь:

 

 

Використовуючи математичний пакет Matlab, виконаємо розрахунок. Нижче наведено лістинг програмного модуля:

clear;

load wt1.mat;

for i=1:8;

    for j=1:8;

        A(i,j)=sum(x.^(i+j-2));

    end;

end;

for i=1:8;

    B(i)=sum(x.^(i-1).*y);

end;

a=B/A;

t=1:150;

f=a(8)*t.^7+a(7)*t.^6+a(6)*t.^5+a(5)*t.^4+a(4)*t.^3+a(3)*t.^2+a(2)*t+a(1);

hPlot=plot(t,f);

hold on;

plot(x,y,'red');

xlabel('t')

ylabel('U(t)')

axis([0,150,-0.2,2.2]);

set (hPlot,'LineWidth',2);

 

Розв’язавши систему, отримаємо  такі значення коефіцієнтів:

 

a₀ = 1,29151435944999;      

a₁ =-0,0272224047072026;

a₂ =0,000865020069750201;

a₃ = -4,08236935634828e-06;

a₄ =-2,21522202872908e-07;

a₅ =3,80702250957075e-09; 

a₆ =-2,28826295649309e-11 ;

a₇ =4,78851337798592e-14.

 

Нижче наведемо графічну модель (рис. 1.1):

Рис. 1.1 - Поліноміальна модель вхідного сигналу U(t) поліномом n-го порядку

 

1.3 Побудова моделі за допомогою полінома Чебишева

За методом найменших  квадратів рівняння апроксимуючої  функції для  полінома Чебишева буде мати вигляд:

 

де  – невідомі коефіцієнти.

 

Запишемо умову метода найменших квадратів:

тобто

 

Продиференціювавши даний вираз по і прирівнявши його до нуля розв’язуємо слідуючу систему рівнянь:

 

 

В системі  функції    вибрані таким чином:

 

        

 

В даних виразах коефіцієнти  визначаються за формулами:

 

    

 

Використовуючи математичний пакет Matlab, виконаємо розрахунок. Нижче наведено лістинг програмного модуля:

clear;

load wt1.mat;

Cheb_por=7;

n=110;

fi=ones(n,7);

for i=1:n

fi(i,2)=x(i)-(1/n)*sum(x);

end;

for p=3:Cheb_por;i=1:n;

beta(p)=-sum(x(i).*fi(i,p-1).^2)/sum(fi(i,p-1).^2);

gama(p)=-sum(x(i).*fi(i,p-1).*fi(i,p-2))/sum(fi(i,p-1).^2);

fi(i,p)=(x(i)+beta(p)).*fi(i,p-1)+gama(p).*fi(i,p-2);

end;

A=fi'*fi;

B=fi'*y;

a=A\B;

for i=1:n;j=1:7;

Cheb(i)=a(j)'*fi(i,j)';

end;

hPlot=plot(x,Cheb);

hold on;

plot(x,y,'red');

xlabel('t');

ylabel('U(t)');

set(hPlot,'LineWidth',2);

 

Розв’язавши систему, такі значення коефіцієнтів:

        a₁= 0,974882715534506;

a₂=-0,00347405766549091 ;

a₃=7,36966019328142e-05;

a₄=2,76744747005052e-06;

a₅=-2,26238197448055e-08;

a₆=-4,75654175893397e-10;

a₇=2,57745717420492e-12.

 

На рис. 1.2 наведена модель вхідного сигналу U(t), побудована поліномом Чебишева.

 

Рис. 1.2 - Модель вхідного сигналу U(t), побудована поліномом Чебишева

 

1.4 Побудова моделі за допомогою перетворення Фур’є

Наведемо коротку характеристику перетворень Фур’є.

Вираз виду:

 

називаються тригонометричними  поліномами (ступеня π). Виявляється, що будь-яка безперервна (і не тільки безперервна) функція на проміжку (-π,π) може бути апроксимована з попередньо заданою точністю тригонометричними поліномами з відповідними коефіцієнтами a₀ ,…,a_n, та b₀ ,…, b_n.

Лінійні комбінації функцій  виду і також називають тригонометричними поліномами, мають період 2πТ. Тому можна так переформулювати наше твердження: будь-яка періодична функція може бути апроксимована тригонометричними поліномами.

Припустимо, що значення функції f(x), виміряються в деякому експерименті й за даними y_1,…,y_m нам вдалося чисельно оцінити коефіцієнти Фур'є функції f(x).

Часто виявляється, що функція:

 

 

Тоді для знаходження  коефіцієнтів Фур’є використовується наступний метод. Нехай функція  , що має період , представлена рядом Фур’є:

 

Якщо формально помножити  дану рівність на або і почленно проінтегрувати від до , то зайві гармоніки зникнуть і отримаємо формулу для коефіцієнтів:

 

,  

 

В нашому випадку рівняння апроксимуючої лінії буде:

 

де m =7.

Запишемо формули для  розрахунку коефіцієнтів:

 

   

 

Використовуючи математичний пакет Matlab, виконаємо розрахунок. Нижче наведено лістинг програмного модуля:

 

>> clear;

>> load wt1.mat;

>> n=110;

>> Fur_por=4;

>> m=Fur_por;

>> for k=1:m+1

A=(1/n)*sum(y);

B(k)=(2/n)*sum(y.*cos(k*x));

C(k)=(2/n)*sum(y.*sin(k*x));

end

>> for k=1:m+1

for i=1:n

FS(i,1)=A+sum(B(k)*sin(k*x(i))+C(k)*cos(k*x(i)));

end

end

>> plot(x,y,'red')

>> hold on;

>> hPlot=plot(x,FS);

>> xlabel('t')

>> ylabel('U(t)')

>> axis([0,150,-0.2,2.2]);

>> set(hPlot,'LineWidth',2);

Розв’язавши систему, отримаємо  таку модель:

 

Рис. 1.3 - Модель вхідного сигналу U(t), побудована за допомогою перетворення  Фур’є

 

і такі значення коефіцієнтів:

 

 

 

 

1.5 Статистична обробка даних

Розрахуємо коефіцієнти  дисперсії, кореляції та коваріації.

Визначимо коефіцієнти дисперсії  значень часу та функції, за такими формулами:

    

 

Коефіцієнт коваріації за такою формулою:

 

 

Коефіцієнт кореляції  за такою формулою:

 

 

Використовуючи математичний пакет Matlab, виконаємо розрахунок і нижче наведено лістинг програмного модуля:

>> clear;

>> load wt1.mat;

>> n=110;

>> i=1:n;

>> tser=sum(x)/n;

>> fser=sum(y)/n;

>> cov=sum((x-tser).*(y-fser))/(n-1)

>> Dt=sum((x-tser).^2)/(n-1);

>> Du=sum((y-fser).^2)/(n-1)

>> r=cov/sqrt(Dt*Du)

Після виконаного розрахунку отримаємо:

 

cov = -0.2308;

Dt =1.4933e+003;

Du = 0.0414;

r = -0.0294.

 

1.6 Знаходження періодограми сигналів

Періодограма сигналу – це оцінка спектральної щільності потужності (СЩП) сигналу по певній кількості відліків реалізації випадкового процесу, заснована на обчисленні квадрата модуля перетворення Фур'є послідовності даних з використанням статистичного усереднення:

 

 

де X_T (iω) – перетворення Фур'є функції x(t) на кінцевому часовому інтервалі, T_r –інтервал фінітності, E – оператор статистичного усереднення (математичне очікування).

Слід зауважити, що дане співвідношення виконується тільки при нескінченному числі використовуваних відрахунків, тому при будь-якому кінцевому N періодограмна оцінка спектральної щільності потужності виявляється зміщеною – виходить, що всередині суми співвідношення кореляційна функція сигналу множиться на трикутну вагову функцію.

Для побудови періодограми сигналу розрахуємо середній спектр потужності , що розраховується за такими залежностями:

де n – загальна кількість точок.

 

Використовуючи математичний пакет Matlab, виконаємо розрахунок і побудуємо періодограму.

Нижче наведено лістинг програмного  модуля:

 

>> clear;

>> load matlabwt1_1.mat;

>> n=106;

>> for k=1:n

w(k,1)=(2*pi*(k))/n;

U(k,1)=(1/sqrt(n))*sum(y.*exp(-j*w(k).*x));

end

>> Ui=imag(U);

>> Ur=real(U);

>> subplot(1,2,1);

>> plot(w,Ui);

>> title('a)');

>> xlabel('Нормована частота(Гц)');

>> ylabel('Спектр щільності(Дб/Гц)');

>> grid on

>> subplot(1,2,2);

>> plot(w,Ur);

>> title('б)');

>> xlabel('Нормована частота(Гц)');

>> ylabel('Спектр щільності(Дб/Гц)');

>> grid on

 

Рис. 1.4 - Періодограма вхідного сигналу U(t), де а) – уявна частина,

б) – дійсна частина 

 

1.7 Вибір оптимальної моделі

Оптимальна модель вибирається  згідно критерію, що відповідає похибці  побудови моделі до вихідної характеристики в співвідношенні до кількості аналізованих даних:

Використовуючи математичний пакет Matlab, виконаємо розрахунок і запишемо результати до таблиці 5.2.

Нижче наведено лістинг програмного  модуля:

 

for i=1:a;n=104;

    Spol_Vt=sqrt(sum((y(i)-f(i)).^2/(n-2)))

end;

SCheb_Ut=sqrt((sum(y-Cheb').^2)/(n-2))

SFur_Ut=sqrt((sum(y-FS).^2)/(n-2))

 

Таблиця 1.2

Оптимальні значення моделей

Вид функції	Поліном	Поліном Чебишева	Функція Фур’є
Значення критерію	0.0397	3.5254e-016	0.0109

 

Як видно з таблиці  1.2, значення критерію мінімальне у тому випадку, коли апроксимуючою функцією виступає поліном Чебишева, отже модель побудована даним чином буде оптимальною.

 

 

 

 

 

 

Розділ  2. Обробка вхідного сигналу V(t), знаходження апроксимуючої функції та вибір оптимальної моделі

 

2.1 Табулювання сигналу

Табулювання графіка функції (сигналу) проводимо тим же способом як і для вхідного сигналу U(t), що описаний в другому розділі.

Після того як ми отримаємо  дві матриці даних: матриця даних зміни часу – xt та  матриця даних зміни значення сигналу залежно від часу – y(vt), складемо таблицю табулювання.

 

Таблиця 2.1

Табулювання вхідного сигналу V(t)

Номер точки	Сигнал  V(t)		Номер точки	Сигнал  V(t)
	t	V(t)		t	V(t)
1	4,30220687	1,2786299	15	24,75142581	0,954858844
2	5,51301588	0,66844599	16	25,55863182	0,830331516
3	6,58929056	1,11363119	17	26,76944084	1,449854971
4	7,80009958	0,49099455	18	27,98024986	0,627974609
5	10,2217176	1,61796686	19	29,05652454	0,914387462
6	12,5088013	1,08249935	20	30,53640222	0,802312867
7	13,7196103	0,98599068	21	31,61267691	0,936179745
8	14,6613507	1,17589485	22	32,82348592	0,886368813
9	15,8721597	1,37825176	23	33,8997606	1,194573949
10	17,217503	1,08872572	24	35,24510396	0,780520585
11	18,6973807	1,51211864	25	36,32137864	0,839671066
12	20,7153958	0,5003341	26	38,60846234	1,368912208
13	22,4643421	0,889482	27	39,95380569	0,543918662
14	23,2715481	1,01400932	28	41,03008037	1,586835032

 

Продовження таблиці 2.1

Номер точки	Сигнал  V(t)		Номер точки	Сигнал  V(t)
	t	V(t)		t	V(t)
29	42,3754237	0,82410515	58	78,16155688	1,166555301
30	43,3171641	1,26306398	59	79,10329722	0,805426051
31	45,7387821	1,45608134	60	80,44864057	1,428062689
32	47,0841255	0,93617974	61	81,65944959	1,082499355
33	48,1604001	1,07938617	62	83,94653329	1,163442118
34	49,5057435	1,02957524	63	86,23361699	1,505892269
35	50,8510868	0,90193473	64	88,65523502	0,571937311
36	51,6582928	1,15721575	65	89,59697537	1,350233109
37	52,6000332	1,18212122	66	91,07685305	0,8147656
38	54,0799109	0,97353794	67	92,15312773	1,418723139
39	55,4252542	1,1291971	68	92,82579941	1,113631187
40	56,5015289	1,82966332	69	93,36393675	0,911274279
41	57,8468723	1,07938617	70	94,44021143	1,46853407
42	58,9231469	1,26306398	71	95,92008912	0,830331516
43	60,2684903	0,93929293	72	96,72729513	0,749388753
44	61,2102306	0,92684019	73	98,20717282	1,34089356
45	62,4210397	0,98599068	74	99,2834475	0,684011906
46	63,766383	1,23193215	75	100,6287909	1,275516712
47	65,2462607	0,76495467	76	101,7050655	1,194573949
48	66,0534667	0,53457911	77	102,7813402	1,026462057
49	67,2642757	1,08872572	78	104,2612179	0,976651126
50	68,2060161	0,96108521	79	105,2029582	1,670890978
51	69,6858938	1,15098938	80	106,5483016	0,973537943
52	70,7621684	1,10117845	81	107,893645	1,020235691
53	72,1075118	1,27551671	82	109,104454	0,799199684
54	72,9147178	1,33778038	83	111,3915377	0,917500645
55	74,5291298	0,68712509	84	112,333278	1,574382299
56	75,7399388	0,64042734	85	113,6786214	0,942406111
57	76,8162135	1,21325305	86	114,8894304	0,877029264