Giáo trình Matlab cơ bản - Chương 3: Nội suy và xấp xỉ hàm

g một  tín hiệu  tương  tự,  tuỳ  thuộc vào 
kích thước DFT, chu kì  lấy mẫu, khoảng  lấy mẫu và đệm zero. So sánh phố 
khi giảm chu kì lấy mẫu T từ 0.1s đến 0.05s 
 232
3. Nội suy bằng các dùng biến đổi Fourrier rời rạc: Ta dùng DFS/DFT để nội 
suy dãy x[n] nhận được từ kết quả lấy mẫu tín hiệu ở khoảng cách cách đều. 
Thủ tục gồm hai bước: lấy N điểm FFT X(k)  của x[n] và dùng công thức: 
j2 kt /NT
|k| N/ 2
N / 2 1
j2 kt /NT
k 1
1xˆ(t) X(k)e
N
1 X(0) 2 Real X(k)e X(N/ 2)cos( /T)
N
π
<
− π
=
=
⎧ ⎫⎡ ⎤= + + π⎨ ⎬⎣ ⎦⎩ ⎭
∑
∑
%
  (5) 
Ta xây dựng hàm nội suy interpdfs(): 
function [xi,Xi] = interpdfs(T, x, Ws, ti) 
%T : chu li lay mau 
%x : thu tu roi rac hoa 
%Ws: tan so dung chuan (1.0 = pi[rad]) 
%ti: khoang thoi gian noi suy 
if nargin < 4 
    ti = 5;  
end 
if nargin  1 
    Ws = 1;  
end 
N = length(x); 
if length(ti) == 1 
    ti = 0:T/ti:(N‐1)*T; %khoang con duoc chia cho ti 
end 
ks = ceil(Ws*N/2); 
Xi = fft(x); 
Xi(ks + 2:N ‐ ks) = zeros(1,N ‐ 2*ks ‐ 1); %pho da loc 
xi = zeros(1,length(ti)); 
for k = 2:N/2 
    xi = xi+Xi(k)*exp(j*2*pi*(k ‐ 1)*ti/N/T); 
end 
xi = real(2*xi+Xi(1)+Xi(N/2+1)*cos(pi*ti/T))/N; %pt.(.5) 
Để nội suy ta dùng chương trình ctfourier.m: 
clear, clf 
 233
w1 = pi;  
w2 = .5*pi; %hai tan so 
N = 32;  
n = [0:N ‐ 1];  
T = 0.1;  
t = n*T; 
x = sin(w1*t)+0.5*sin(w2*t)+(rand(1,N) ‐ 0.5); %0.2*sin(20*t); 
ti = [0:T/5:(N ‐ 1)*T]; 
subplot(411), plot(t,x,ʹk.ʹ) %so lieu ban dau 
title(ʹSo lieu ban dau va ket qua noi suyʹ) 
[xi,Xi] = interpdfs(T,x,1,ti); 
hold on, plot(ti,xi,ʹrʹ) %tai tao tin hieu 
k = [0:N ‐ 1]; 
subplot(412), stem(k,abs(Xi),ʹk.ʹ) %pho ban dau 
title(ʹPho ban dauʹ) 
[xi,Xi] = interpdfs(T,x,1/2,ti); 
subplot(413), stem(k,abs(Xi),ʹr.ʹ) %pho da  loc 
title(ʹPho da locʹ) 
subplot(414), plot(t,x,ʹk.ʹ, ti,xi,ʹrʹ) %tin hieu da loc 
title(ʹTin hieu da locʹ) 
§9. XẤP XỈ HÀM BẰNG PHƯƠNG PHÁP BÌNH PHƯƠNG BÉ NHẤT 
1. Khái niệm chung: Trong các mục trước ta đã nội suy giá trị của hàm. Bài 
toán đó là cho một hàm dưới dạng bảng số và phải tìm giá trị của hàm tại một 
giá trị của đối số không nằm trong bảng. 
  Trong thực tế, bên cạnh bài toán nội suy ta còn gặp một dạng bài toán 
khác. Đó là tìm công thức thực nghiệm của một hàm.  
Nội dung bài toán là từ một loạt các điểm cho trước (có thể là các giá trị 
của một phép đo nào đó) ta phải tìm một hàm xấp xỉ các giá trị đã cho. Ta sẽ 
dùng phương pháp bình phương tối thiểu để giải bài toán.  
Giả  sử  có mẫu quan  sát  (xi, yi)  của hàm y =  f(x). Ta  chọn hàm  f(x)  có 
dạng:  
  f(x) = a0f0(x) + a1f1(x) + a2f2(x)...            (1) 
Trong đó các hàm f0(x), f1(x), f2(x) v.v. là (m+1) hàm độc lập tuyến tính mà ta 
có thể chọn tuỳ ý và các hệ số ai là tham số chưa biết mà ta phải xác định dựa 
 234
vào hệ hàm đã chọn và các điểm quan sát. Sai số giữa trị đo được và trị tính 
theo (1) là : 
  ei = yi ‐ f(xi)                   (2) 
Sai  số này  có  thể âm hay dương  tuỳ  từng giá  trị  của yi. Khi dùng phương 
pháp bình phương bé nhất ta xét bình phương của sai số tại một điểm: 
  [ ]2 2i i ie y f(x )= −                    (3) 
Với n điểm tổng bình phương của sai số sẽ là : 
  [ ]{ }n n 22i i 0 0 i 1 1 i m m i
i 1 i 1
S e y a f (x ) a f (x ) a f (x )
= =
= = − + + ⋅ ⋅ ⋅ +∑ ∑  
Rõ ràng S là hàm của các giá trị cần tìm ai và chúng ta sẽ chọn các ai sao 
cho S đạt giá trị min, nghĩa là các đạo hàm 
ia
S
∂
∂ phải bằng không.  
Ta sẽ xét các trường hợp cụ thể. 
2. Hàm xấp xỉ có dạng đa thức: Trong trường hợp tổng quát ta chọn hệ hàm 
xấp xỉ là một đa thức, nghĩa là: 
  f(x) = a0 + a1x + a2x2 +∙∙∙+ amxm 
Vậy hàm S là : 
  ( )22 mi 0 1 2 mS y a a x a x a x= − + + + ⋅⋅⋅ +    
Theo điều kiện đạo hàm  0
a
S
i
=∂
∂ ta nhận được hệ phương trình: 
⎪⎪
⎪⎪
⎪⎪
⎪
⎩
⎪⎪
⎪⎪
⎪⎪
⎪
⎨
⎧
=+⋅⋅⋅++
⋅⋅⋅
=+⋅⋅⋅++
=+⋅⋅⋅++
=+⋅⋅⋅++
=+⋅⋅⋅++
∑∑∑ ∑
∑∑∑ ∑
∑∑∑ ∑
∑∑∑ ∑
∑∑ ∑
=== =
−
−
=== =
+
−
+
=== =
+
−
+
=== =
−
+
== =
−
−
n
1i
i
m
i
n
1i
m
i0
n
1i
n
1i
1m2
i1m
m2
im
n
1i
i
3
i
n
1i
3
i0
n
1i
n
1i
2m
i1m
3m
im
n
1i
i
2
i
n
1i
2
i0
n
1i
n
1i
1m
i1m
2m
im
n
1i
ii
n
1i
i0
n
1i
n
1i
m
i1m
1m
im
n
1i
i0
n
1i
n
1i
1m
i1m
m
im
yxxaxaxa
yxxaxaxa
yxxaxaxa
yxxaxaxa
ynaxaxa
Đây là một hệ phương trình tuyến tính. Giải nó ta nhận được các gía trị ai. 
Ta xây dựng hàm polynomfit() thực hiện thuật toán trên: 
 235
function x = polyfits(xData, yData, m) 
%Dung de tinh he so cua da thuc xap xi 
% Cu phap: x = polyfits(xData, yData, m) 
m = m+1; 
A = zeros(m);  
b = zeros(m, 1);  
s = zeros(2*m‐1, 1); 
for i = 1:length(xData) 
    temp = yData(i); 
    for j = 1:m 
        b(j) = b(j) + temp; 
        temp = temp*xData(i); 
    end 
    temp = 1; 
    for j = 1:2*m‐1 
        s(j) = s(j) + temp; 
        temp = temp*xData(i); 
    end 
end 
for i = 1:m 
    for j = 1:m 
        A(i, j) = s(i+j‐1); 
    end 
end 
x = A\b; 
% Sap xep lai he so tu so mu cao nhat 
x = flipdim(x, 1); 
Để  xấp  xỉ  một  dãy  số  liệu  bằng  hàm  đa  thức  ta  dùng  chương  trình 
ctpolynomfit.m:  
clear all, clc 
xData = [0 1 2 3 4]; 
yData = [1  8 24  63 124]; 
x = polyfits(xData, yData, 3); 
y = 0:0.1:4; 
 236
z = polyval(xʹ, y); 
hold on 
plot(y, z,ʹ‐bʹ, xData, yData, ʹroʹ); 
3.Hàm dạng Aecx: Khi các số liệu thể hiện một sự biến đổi đơn điệu ta dùng 
hàm xấp xỉ là y = Aecx. Lấy logarit hai vế ta có : 
  lny = lnA + cxlne 
Theo điều kiện đạo hàm  0
a
S
i
=∂
∂ ta có hệ phương trình : 
⎪⎪⎩
⎪⎪⎨
⎧
=+
=+
∑ ∑∑
∑ ∑
= ==
= =
n
1i
n
1i
ii
n
1i
i
2
i
n
1i
n
1i
ii
ylnxxAlnxc
ylnAlnnxc
Giải hệ phương trình này ta có các hệ số A và c.  
Ta xây dựng hàm expfit() để xấp xỉ 
function [c,A] = expfit(x, y) 
a = sum(x); 
b = size(x,2); 
c = sum(log(y)); 
d = sum(x.^2); 
e = sum(x.*log(y)); 
d1 = a*a ‐ d*b; 
d2 = c*a ‐ e*b; 
d3 = a*e ‐ c*d; 
c = d2/d1; 
A = exp(d3/d1); 
Ta dùng chương trình ctexpfit.m để xấp xỉ dãy số liệu đã cho 
clear all, clc 
x = [1.2  2.8  4.3  5.4   6.8  7.9]; 
y = [7.5  16.1  38.9  67  146.6 266.2]; 
[c, A] = expfit(x, y); 
t = 0:0.1:8; 
z = A*exp(c*t); 
plot(t, z, ʹ‐bʹ, x, y, ʹroʹ); 
 237
4. Hàm dạng Axq: Khi các số liệu thể hiện một sự biến đổi đơn điệu ta  cũng 
có thể dùng hàm xấp xỉ là y = Axq. Lấy logarit hai vế ta có: 
  lny = lnA + qlnx 
Theo điều kiện đạo hàm triệt tiêu ta có hệ phương trình : 
⎪⎪⎩
⎪⎪⎨
⎧
=+
=+
∑ ∑∑
∑ ∑
= ==
= =
n
1i
n
1i
ii
n
1i
ii
2
n
1i
n
1i
ii
ylnxlnxlnAlnxlnq
ylnAlnnxlnq
Giải hệ phương trình này ta có các hệ số 
A và q. 
Ta xây dựng hàm powerfit() để xấp xỉ: 
function [q, A] = powerfit(x, y) 
a = sum(log(x)); 
b = size(x, 2); 
c = sum(log(y)); 
d = sum(log(x).^2); 
e = sum(log(x).*log(y)); 
d1 = a*a ‐ d*b; 
d2 = c*a ‐ e*b; 
d3 = a*e ‐ c*d; 
q = d2/d1; 
A = exp(d3/d1); 
Ta dùng chương trình ctpowerfit.m để xấp xỉ dãy số liệu đã cho: 
clc 
x = [  1  2  3  4      5]; 
y = [1.5  15.1     52.5    130.5  253]; 
[q,A] = powerfit(x, y) 
t = 0.1:0.1:5; 
z = exp(log(A)+q*log(t)); 
plot(t, z, ʹ‐bʹ, x, y, ʹroʹ); 
5. Hàm lượng giác: Khi quan hệ y = f(x) có dạng tuần hoàn ta dùng hàm xấp 
xỉ là tổ hợp tuyến tính của các hàm sin và cosin dạng: 
 238
  ∑ ∑
= =
ω+ω+= n
1i
n
1i
ii0 )xisin(b)xicos(aa)x(f  
Để đơn giản trước hết ta xét hàm chỉ có một số hạng sin‐cos, nghĩa là : 
  xsinbxcosaa)x(f 110 ω+ω+=  
Hàm S sẽ có dạng : 
n 2
i 0 1 1
i 1
S y (a a cos x b sin x)
=
= − + ω + ω⎡ ⎤⎣ ⎦∑  
Theo  điều kiện  đạo hàm  triệt  tiêu  ta  có hệ phương  trình  đối với  các hệ  số 
dạng: 
i i 0 i
2
i i i i 1 i i
2
i i i i 1 i i
n cos x sin x a y
cos x cos x cos x sin x a y cos x
sin x cos x sin x sin x b y sin x
⎡ ⎤ ⎡ ⎤ω ω ⎡ ⎤⎢ ⎥ ⎢ ⎥⎢ ⎥ω ω ω ω = ω⎢ ⎥ ⎢ ⎥⎢ ⎥⎢ ⎥ ⎢ ⎥⎢ ⎥ω ω ω ω ω⎣ ⎦⎣ ⎦ ⎣ ⎦
∑ ∑ ∑∑ ∑ ∑ ∑∑ ∑ ∑ ∑
Do: 
i isin x cos x0 0
n n
ω ω= =∑ ∑  
2 2
i i
i i
sin x cos x1 1
n 2 n 2
cos x sin x
0
n
ω ω= =
ω ω =
∑ ∑
∑  
nên hệ phương trình có dạng đơn giản : 
0 i
1 i i
1 i i
n 0 0 a y
0 n 2 0 a y cos x
0 0 n 2 b y sin x
⎡ ⎤⎡ ⎤ ⎡ ⎤ ⎢ ⎥⎢ ⎥ ⎢ ⎥ = ω⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ω⎣ ⎦ ⎣ ⎦ ⎣ ⎦
∑∑∑
Giải hệ ta có : 
  i0 1 i i 1 i i
y 2 2a a y cos x b y sin x
n n n
= = ω = ω∑ ∑ ∑  
Trong trường hợp tổng quát, một cách tương tự ta có: 
∑∑∑ ω=ω== xisinyn
2bxicosy
n
2a
n
ya ii0  
Ta xây dựng hàm sinfit() để xấp xỉ: 
function [a, b, c, omega] = sinfit(x, y, T) 
%T la chu ki 
omega = 2*pi/T; 
n = size(x,2); 
 239
a = sum(y)/n; 
b = (2/n)*sum(y.*cos(omega*x)); 
c = (2/n)*sum(y.*sin(omega*x)); 
Ta dùng chương trình ctsinfit.m để tính: 
c ear all, clc 
x = [0 0.15  0.3  0.45  0.6  0.75  0.9  1.05  1.2  1.3]; 
y = [2.2   1.595 1.031 0.722 0.786 1.2  1.81  2.369 2.678 2.614]; 
T = 1.5; 
[a, b, c, omega] = sinfit(x, y, T) 
t = 0.:0.01:1.5; 
z = a + b*cos(omega*t) + c*sin(omega*t); 
plot(t, z,ʹ‐bʹ, x, y, ʹroʹ); 
6. Hàm hữu tỉ: Khi quan hệ y = f(x) có dạng đường cong bão hoà hay dạng 
arctan, tan v.v ta dùng hàm xấp xỉ là hàm hữu tỉ dạng đơn giản: 
xb
axy +=  
Lấy nghịch đảo của nó ta có : 
a
1
x
1
a
b
y
1 +=  
Đặt 1/y = Y, 1/x = X, b/a = B và 1/a = A phương trình trên sẽ có dạng: 
    Y = A + BX 
và là một đa thức bậc một. Do vậy ta có hệ phương trình đối với các hệ số A 
và B là: 
⎪⎪⎩
⎪⎪⎨
⎧
=+
=+
∑ ∑∑
∑ ∑
= ==
= =
n
1i
n
1i ii
2
i
n
1i i
n
1i
n
1i ii
yx
1
x
1B
x
1A
y
1
x
1BnA
và từ đó tính được a và b.  
Ta xây dựng hàm racfit() để xấp xỉ: 
function [a, b] = racfit(x, y) 
a1 = size(x, 2); 
b1 = sum(1./x); 
c1 = sum(1./y); 
 240
d1 = sum(1./x.^2); 
e1 = sum((1./x).*(1./y)); 
del = a1*d1 ‐ b1*b1; 
del1 = c1*d1 ‐ e1*b1; 
del2 = a1*e1 ‐ b1*c1; 
A = del1/del; 
B = del2/del; 
a = 1/A; 
b = B/A; 
Để xấp xỉ ta dùng chương trình ctracfit.m: 
clear all, clc 
x = [1  2  3  4  5]; 
y = [0.3333333    0.5  0.6    0.66666  0.7142857]; 
[a, b] = racfit(x, y) 
t = 0.:0.01:5; 
z = a*t./(b+t) 
plot(t, z,ʹ‐bʹ, x, y, ʹroʹ);