T1 

Techniki Obliczeniowe i Symulacyjne
Interpolacja, estymacja średniokwadratowa, optymalizacja 

dr inż. Tomasz Zieliński 

2010.03.01

 

Ćwiczenie 1 (1 pkt) 

Interpolacja. W poniższej tabeli przedstawiono wyniki pomiaru temperatury (N=4): 

LP 

1 2 3 4 

Godzina x

1

=5:00

x

2

=6:00

x

3

=8:00

x

4

=11:00 

Temperatura [

o

C] y

1

=1 y

2

=2 y

3

=7 y

4

=15 

Załóżmy, że wielomian interpolujący funkcję z tablicy ma postać: 

( )

1

0

1

1

N

N

y x

a

a x

a

x

−

−

=

+

+K

                                                 (1) 

Dla N punktów pomiarowych mamy N równań: 

( )

,

1,

,

i

i

y x

y

i

N

=

= K

 

które możemy zapisać razem w postaci macierzowej: 

1

1

1

1

1

1

2

2

2

1

1

0

1

1

1

N

N

N

N

N

N

N

a

y

x

x

y

x

x

a

a

y

x

x

−

−

−

−

⎡

⎤ ⎡

⎤ ⎡

⎤

⎢

⎥ ⎢

⎥ ⎢

⎥

⎢

⎥ ⎢

⎥ ⎢

⎥

=

⎢

⎥ ⎢

⎥ ⎢

⎥

⎢

⎥ ⎢

⎥ ⎢

⎥

⎣

⎦ ⎣

⎦

⎣

⎦

T

L

M

L

M

M

M

L M

L

14442444

3

                                              (2) 

Rozwiąż równanie (2) jedną z metod z Lab2 (czyli oblicz a), a następnie skorzystaj z (1) aby obliczyć 
poszukiwaną temperaturę yz dla interesującej nas godziny xz=9:30. Zwróć uwagę, że układ równań (2) 
jest źle uwarunkowany dla dużego N. Macierz T jest macierzą Vandermonde’a: 

T=vander(x); a=inv(T)*y; yz=polyval(a,xz); 

yz=vander([xz zeros(1,N-1)])*a; yz=((xz*ones(1,N)).^(N:-1:0))*a; 

Ćwiczenie 2 (2 pkt) 

Interpolacja - wielomiany interpolujące Lagrange’a i Newtona. Wyznacz analitycznie jaka była 
najbardziej prawdopodobna temperatura o godzinie 9:30 jeśli zastosujemy metodę interpolacji 

Lagrange’a oraz Newtona (sprawdź dla obu metod). Napisz ogólny program dla obu metod interpolacji, 
zakładając, że danych jest N punktów pomiarowych (x

k

, y

k

) oraz interesuje nas przybliżona wartość dla 

x

 zadanego: 

1

zadane

N

x

x

x

≤

≤

. Porównaj otrzymane wyniki z wynikiem uzyskanym w ćw. 1. 

Ćwiczenie 3 (1 pkt) 

Estymacja. Załóżmy, że mamy dyskretny układ liniowy, którego dane wyjściowe 

y

n

 są wynikiem splotu 

danych wejściowych 

x

n

 

z próbkami odpowiedzi impulsowej układu 

h

n

 (w naszym przypadku tylko dwóch, 

h

0

 i 

h

1

), do którego dodaje się szum układu (kanału transmisyjnego): 

1

0

1

1

1

0

1

1

2

1

1

1

0

0

1

0

2

0

2

2

1

2

2

2

2

1

1

2

1

3

1

3

3

2

3

3

3

2

4

4

4

3

,

,

,

x

x

y

s

x

x

y

s

x

x

y

x

x

h

s

h

y

h

s

x

x

y

s

y

x

x

h

s

h

y

h

s

x

x

y

s

x

x

y

s

x

x

⎡

⎤

⎡ ⎤

⎡ ⎤

⎡

⎤

⎡ ⎤

⎡ ⎤

⎢

⎥

⎢ ⎥

⎢ ⎥

⎡ ⎤ ⎡

⎤ ⎡ ⎤ ⎡ ⎤

⎡ ⎤

⎡ ⎤

⎢

⎥

⎢ ⎥

⎢ ⎥

⎢

⎥

⎢ ⎥

⎢ ⎥

=

+

=

+

=

+

⋅

⎢ ⎥ ⎢

⎥ ⎢ ⎥ ⎢ ⎥

⎢ ⎥

⎢ ⎥

⎢

⎥

⎢ ⎥

⎢ ⎥

⎢

⎥

⎢ ⎥

⎢ ⎥

⎣ ⎦ ⎣

⎦ ⎣ ⎦ ⎣ ⎦

⎣ ⎦

⎣ ⎦

⎢

⎥

⎢ ⎥

⎢ ⎥

⎢

⎥

⎢ ⎥

⎢ ⎥

⎣ ⎦

⎣ ⎦

⎣

⎦

⎣ ⎦

⎣ ⎦

⎣

⎦

y = X h

+ s

 

(3) 

Znamy macierz 

X

 (dane wejściowe) i wektor 

y

 (dane wyjściowe). Chcemy znaleźć 

h

, dla którego 

wartość  błąd=

(

) (

)

min

T

T

=

=

s s

y - Xh

y - Xh

. Obliczamy pochodną błędu względem 

h

 i przyrównujemy 

ją do zera. Z tego równania wyznaczamy optymalną wartość 

h

. 

2 

3 

Minimum błędu kwadratowego otrzymujemy dla: 

(

)

T

T

est

-1

h

= X X

X y

                                                            (4) 

Macierz 

(

)

T

-1

X X

X

jest macierzą pseudo-odwrotną macierzy X. Ponieważ wektor 

h

 składa się z 2 liczb, 

najmniejsze wymiary macierzy 

X

 to 2×2 (układ dwóch równań z dwoma niewiadomymi). Im więcej 

przeprowadzamy pomiarów, tym macierz 

X

 ma więcej wierszy, jej macierz pseudo-odwrotna ma więcej 

kolumn i rozwiązanie 

h

est

 jest dokładniejsze. 

Przyjmij: 

M=2; h =(M:-1:1)’; 

a potem N równe kolejno 2, 3, 4, …25. Dla każdej wartości 

N

 podstaw: 

s=randn(N,1); x=rand(N+(M-1),1);

 

zbuduj macierz 

X

 i rozwiąż równanie (4). Narysuj jak się zmieniały elementy wektora 

h

est

 kiedy liczba 

punktów pomiarowych rosła. Powtórz eksperyment dla M=4. 
 

Ćwiczenie 4 (0.5 pkt) 

Zauważ, że równanie interpolacji (2) może być także przedstawione jako zadanie estymacji 
(aproksymacji) średniokwadratowej, jeśli założymy że zmierzone wartości 

y

 są zaszumione przez szum 

s

: 

1

1

1

1

1

1

1

2

2

2

2

1

1

0

1

1

1

N

N

N

N

N

N

N

N

y

a

s

x

x

y

s

x

x

a

y

a

s

x

x

−

−

−

−

⎡

⎤

⎡

⎤

⎡

⎤ ⎡ ⎤

⎢

⎥

⎢

⎥

⎢

⎥ ⎢ ⎥

⎢

⎥

⎢

⎥

⎢

⎥ ⎢ ⎥

=

+

⎢

⎥

⎢

⎥

⎢

⎥ ⎢ ⎥

⎢

⎥

⎢

⎥

⎢

⎥ ⎢ ⎥

⎣

⎦

⎣

⎦ ⎣ ⎦

⎣

⎦

X

L

M

L

M

M

M

M

L M

L

14442444

3

 

Wówczas optymalne wartości wektora 

a

, minimalizujące błąd średniokwadratowy są równe: 

(

)

T

T

est

-1

a

= X X

X y

 

Oblicz je i porównaj z poprzednio wyznaczonymi w ćw. 1 i 2. 

Ćwiczenie 5 (0.5 pkt) 

Optymalizacja. Załóżmy, że dysponujemy N=1000 próbkami sygnału 

y(n)

, który jest zaszumionym 

sygnałem 

x(t) 

o znanej postaci (

s(t)

 – szum): 

0

( )

( )

( )

cos(2

)

( )

t

y t

x t

s t

Ae

f t

s t

α

π

φ

−

=

+

=

+

+

 

Sygnał ten został spróbkowany z częstotliwością 

f

pr

 

= 1000 Hz (dt=1/fpr; t = dt*(0:N-1)). Załóżmy, 

że: 

A =4.5, f

0

 = 10 Hz, 

α = 0.5*f

0

, 

φ = 0; 

s(t) to szum gaussowski, dla którego stosunek sygnału do szumu jest równy 38 dB (s = awgn(x,38);). 

Wygeneruj sygnał 

y(t) 

oraz za pomocą funkcji Matlaba  fminsearch()  znajdź optymalne wartości 

parametrów {A, α, f

0

, φ} sygnału 

x(t)

 (minimalizujących błąd odległości punktów eksperymentalnych od 

krzywej teoretycznej: 

0

( )

cos(2

)

t

x t

Ae

f t

α

π

φ

−

=

+

 

Skorzystaj z programów cw_7.m i cw_8.m z przedmiotu MiTP-II. 

1

0

1

1

1

0

1

1

2

1

1

1

0

0

1

0

2

0

2

2

1

2

2

2

2

1

1

2

1

3

1

3

3

2

3

3

3

2

4

4

4

3

,

,

,

x

x

y

s

x

x

y

s

x

x

y

x

x

h

s

h

y

h

s

x

x

y

s

y

x

x

h

s

h

y

h

s

x

x

y

s

x

x

y

s

x

x

⎡

⎤

⎡ ⎤

⎡ ⎤

⎡

⎤

⎡ ⎤

⎡ ⎤

⎢

⎥

⎢ ⎥

⎢ ⎥

⎡ ⎤ ⎡

⎤ ⎡ ⎤ ⎡ ⎤

⎡ ⎤

⎡ ⎤

⎢

⎥

⎢ ⎥

⎢ ⎥

⎢

⎥

⎢ ⎥

⎢ ⎥

=

+

=

+

=

+

⋅

⎢ ⎥ ⎢

⎥ ⎢ ⎥ ⎢ ⎥

⎢ ⎥

⎢ ⎥

⎢

⎥

⎢ ⎥

⎢ ⎥

⎢

⎥

⎢ ⎥

⎢ ⎥

⎣ ⎦ ⎣

⎦ ⎣ ⎦ ⎣ ⎦

⎣ ⎦

⎣ ⎦

⎢

⎥

⎢ ⎥

⎢ ⎥

⎢

⎥

⎢ ⎥

⎢ ⎥

⎣ ⎦

⎣ ⎦

⎣

⎦

⎣ ⎦

⎣ ⎦

⎣

⎦

y = X h

+ s

 

 

 

 

 

X 

h 

y 

= 

+ 

s 

(

)

T

T

est

-1

h

= X X

X y

 

 
 

 

 

 

 

 

 
 

 

 

= 

(X

T

X)

-1

X

T

 

h 

= 

X

T

X 

X

T

-1

h 

 

X 

= 

X

T 

X

T

h 

-1

y 

y 

y