ĐẠI HỌC ĐÀ NẴNG
TRƯỜNG ĐẠI HỌC BÁCH KHOA
oOo
TIỂU LUẬN MÔN HỌC
XỬ LÝ TÍN HIỆU SỐ
ĐỀ TÀI:
THIẾT KẾ BỘ LỌC VI PHÂN THEO CẤU TRÚC FIR
BẰNG PHƯƠNG PHÁP LẤY MẪU TẦN SỐ
GVHD : TS. NGÔ VĂN SỸ
SVTH : Phạm Đức Linh
: Nguyễn Văn Lượm
: Hoàng Hữu Vấn
Lớp : 07CLC2
Đà Nẵng, tháng 11/2010
2
LỜI MỞ ĐẦU
Xử lý tín hiệu số (Digital Signal Processing – DSP) là một lĩnh vực mới, đã phát
triển mạnh mẽ về mặt lý thuyết và công nghệ. Ngày nay, nó đã trở thành môn học
không thể thiếu cho nhiều ngành học về khoa học và kỹ thuật. Các thuật toán của nó
đã mang lại những thuận tiện cho công nghệ và ứng dụng mới, nhất là các lĩnh vực về
kỹ thuật điện tử và điều khiển tự động.
Phép xử lý cơ bản nhất của DSP là lọc, và các hệ thống được đề cập đến nhiều nhất
trong xử lý tín hiệu số là các bộ lọc số (Digital Filter).
Nếu xét về đáp ứng xung có thể chia các bộ lọc số thành 2 loại chính là bộ lọc có
đáp ứng xung hữu hạn FIR (Finite Impulse Response) còn gọi là lọc không đệ quy, và
bộ lọc có đáp ứng xung vô hạn IIR (Infinte Impulse Response) còn gọi là lọc đệ quy.
Xét về đáp ứng tần số biên độ có thể chia các bộ lọc, FIR hay IIR, thành 4 loại cơ
bản: thông thấp, thông cao, thông dải và chắn dải. Các bộ lọc này có thể được thiết kế
bằng những phương pháp sau đây: Phương pháp cửa sổ (Window Design Techniques),
Phương pháp lấy mẫu tần số (Frequency Sampling Design Techniques) và Phương
pháp xấp xỉ tối ưu cân bằng gợn sóng (Optimal Equiripple Design Techniques). Mỗi
tần số ở mỗi chu kỳ về số lượng và đáp ứng pha của bộ lọc mong muốn. Nói chung
đáp ứng pha cần tuyến tính trong dải thông.
- Trong trường hợp các bộ lọc FIR , Có thể đạt được chính xác yêu cầu về pha tuyến
tính.
- Trong trường hợp các bộ lọc IIR, một dải thông có pha tuyến tính là rất khó đạt.
Do đó, chúng ta chỉ xét các chỉ tiêu về biên độ.
Có 2 nhóm chỉ tiêu:
• Các chỉ tiêu tuyệt đối (Absolute Specifications) : Đưa ra một tập các yêu cầu
trên hàm đáp ứng biên độ |H(e
jw
)|. Những chỉ tiêu này được sử dụng chung
cho các bộ lọc FIR.
• Các chỉ tiêu tương đối (Relative Specifications - DB) : Đưa ra các yêu cầu
tính theo decibels (dB), được cho bởi :
Phương pháp này được là một phương pháp phổ biến nhất trong thực tế và được sử
dụng cho cả bộ lọc FIR và IIR.
Trong những
phần tiếp sau đây, để mô tả các chỉ tiêu của bộ lọc chúng ta sẽ xem xét việc thiết kế
một bộ lọc thông thấp như một ví dụ cơ sở để thiết cho bộ lọc FIR.
1.2. Các chỉ tiêu
Trên hình (1.1) là mô tả các chỉ tiêu của bộ lọc FIR thông thấp (Low Pass Filter):
1.2.1.Các chỉ tiêu tuyệt đối
•Band [0, w
p
] được gọi là dải thông, và δ
1
là dung sai (gợn sóng) được chấp
nhận trong đáp ứng dải thông lý tưởng.
0
)(
: Độ gợn sóng trong dải thông tính theo dB.
A
s
: Suy hao trong dải chắn tính theo dB.
Quan hệ giữa các chỉ tiêu này như sau:
Tại sao tập trung vào bộ lọc thông thấp?
Các chỉ tiêu trên được đưa ra đối với bộ lọc FIR thông thấp, và tất nhiên đối với các
bộ lọc khác như thông cao HPF (High Pass Filter), thông dải BPF (Band Pass Filter)
đều có thể được định nghĩa tương tự. Tuy nhiên, các tham số thiết kế quan trọng nhất
là các dung sai dải tần và các tần số cạnh-dải. Bởi vậy, trong phần 1 về cơ sở lý thuyết
này chúng ta chỉ tập trung vào bộ lọc FIR thông thấp. Việc thiết kế cụ thể cho bộ lọc
FIR thông dải bằng kỹ thuật cửa sổ sẽ được phát triển trên cơ sở lọc thông thấp và sẽ
được mô tả chi tiết trong phần 2.
Các thuận lợi trong thiết kế và thi hành lọc số FIR
Việc thiết kế và thực hiện lọc FIR có những thuận lợi sau đây:
(1.2)
(1.3)
0
1
1
2
10
log20 >
+
−=
δ
δ
s
A
• Bộ lọc pha tuyến tính không có méo trễ nhóm và chỉ bị trễ một khoảng không
đổi.
• Đối với bộ lọc có chiều dài M (hoặc bậc M-1) số phép toán có bậc M/2 như đã
khảo sát trong thi hành pha tuyến tính.
1.3. Cấu trúc của bộ lọc FIR
Cho h(n), n=0,1,…,M-1 là đáp ứng xung có chiều dài M. Thì hàm truyền hệ thống là
một bộ lọc đáp ứng xung hữu hạn với hàm hệ thống có dạng:
∑
−
=
−−
−
−
=+++=
1M
0n
n
n
M1
1M
1
10
zbzbzbb)z(H
(1.4)
Như vậy đáp ứng xung h(n) là:
−≤≤
=
(1.7)
Do mẫu thức bằng đơn vị nên ta chỉ có một cấu trúc dạng trực tiếp duy nhất. Cấu
trúc dạng trực tiếp được cho trong hình (1.2) với M = 5:
1.3.2. Cấu trúc dạng ghép tầng
b
0
z
-
1
b
1
z
-
1
b
2
z
-
1
b
3
z
-
1
b
4
y(n)
x(n)
Hình (1.2) Cấu trúc lọc FIR dạng trực tiếp
b
b
1bzbzbb)z(H
(1.8)
∏
=
−−
++=
K
1k
2
2,k
1
1,k0
)zBzB1(b
trong đó
=
2
M
K
, B
k,1
và B
0110
+−++−++−+=
++−+−++−+= )]2Mn(x)1n(x[b)]1Mn(x)n(x[b
10
Sơ đồ khối thực hiện phương trình sai phân trên được mô tả trong hình (1.4) dưới
đây đối với cả M lẻ và M chẵn:
Đối với M lẻ: M = 7, còn đối với M chẵn: M = 6 ( hình 1.4 )
Rõ ràng, với cùng một bậc của bộ lọc (cùng M) cấu trúc pha tuyến tính sẽ tiết kiệm
được 50% các bộ nhân so với cấu trúc dạng trực tiếp.
B
1,1
z
-
1
z
-
1
z
-
1
y(n)
x(n)
B
2,1
B
3,1
b
0
B
0n
n1M)1M(
1M
0n
n
z)n(hzz)n(h)z(H
(1.12)
có (M-1) điểm cực ở gốc (trivial poles) và M-1 điểm không nằm ở vị trí bất kỳ trên
mặt phẳng z. Đáp ứng tần số là:
π≤ω<π−=
∑
−
=
ω−ω
,e)n(h)e(H
1M
0n
njj
(1.13)
1.4.1. Đáp ứng xung h(n)
Chúng ta có thể đưa ra ràng buộc pha tuyến tính:
π≤ω<π−αω−=∠
ω
,)e(H
j
(1.14)
trong đó: α là một hằng số trễ pha. Ta đã biết rằng h(n) phải đối xứng, nghĩa là:
2
1M
,1Mn0),n1M(h)n(h
π≤ω<π−
(1.16)
Đáp ứng pha là đường thẳng nhưng không đi qua gốc. Trong trường hợp này α
không phải là hằng số trễ pha, nhưng:
α−=
ω
∠
ω
d
)e(Hd
j
(1.17)
là hằng số, chính là trễ nhóm (α là một hằng số trễ nhóm). Trong trường hợp này, các
tần số được làm trễ với một tốc độ không đổi.
Đối với kiểu pha tuyến tính này, có thể thấy rằng:
1Mn0),n1M(h)n(h −≤≤−−−=
và
2
,
2
1M π
±=β
−
=α
(1.18)
Điều này có nghĩa rằng đáp ứng xung h(n) là phản đối xứng (antisymmetric). Chỉ số
đối xứng vẫn là
2
1M −
=α
=α
không phải là một số nguyên. Đáp
ứng xung được mô tả trong hình 1.8.
1.4.2. Đáp ứng tần số H(e
j
ω
)
Như vậy, khi tổ hợp hai loại đối xứng và phản đối xứng với M chẵn và M lẻ, ta có
bốn kiểu lọc FIR pha tuyến tính. Đáp ứng tần số của mỗi kiểu có biểu thức và hình
dạng riêng. Để nghiên cứu các đáp ứng pha của các kiểu này, ta viết biểu thức của
H(e
j
ω
) như sau:
2
1M
,
2
;e)e(H)e(H
)(jj
r
j
−
=α
π
±=β=
αω−βωω
(1.19)
trong đó H
r
−
=
ω
ω=
∑
(1.20)
trong đó:
( )
−
=
2
1M
h0a
với mẫu ở chính giữa
(1.21)
( )
, nhưng
2
1M −
=α
không phải là một biến nguyên, thì ta có thể chứng tỏ rằng:
( )
( )
2/1Mj
2/M
1n
j
e
2
1
ncosnb)e(H
−ω−
=
ω
−ω=ω
2/M
1n
r
2
1
ncosnb)(H
(1.24)
Lưu ý: Tại
π=ω
, ta có
( )
0
2
1
ncosnb)(H
2/M
1n
( ) ( )
n1Mhnh −−−=
,
1Mn0 −≤≤
, và
0
2
1M
h =
−
thì ta có thể chứng tỏ:
( )
( )
ω
−
−
= n
2
1M
h2nc
với
2
M
, ,2,1n =
(1.26)
So sánh (1.24) và (1.18), ta có:
( ) ( )
( )
∑
−
=
ω=ω
2/1M
0n
r
nsinncH
(1.27)
Lưu ý: Tại
0
=ω
và
n1Mhnh −−−=
,
1Mn0 −≤≤
, nhưng
2
1M −
=α
không phải là một biến nguyên, thì ta có thể chứng tỏ rằng:
11
( )
ω
−
−
π
=
ω
trong đó:
( )
−= n
2
M
h2nd
với
2
M
, ,2,1n =
(1.29)
So sánh (1.27) và (1.18), ta có:
( )
∑
=
FIR Type 2
FIR Type 3
FIR Type 4
12
PHẦN 1
1.PHUƠNG PHÁP CỬA SỔ - PHẠM ĐỨC LINH
a)LÝ THUYẾT
Ý tưởng cơ bản: chọn một bộ lọc chọn tần lý tưởng (mà đáp ứng xung
luôn luôn phi nhân quả, dài vô hạn) và cắt (lấy cửa sổ) đáp ứng xung của
nó để thu được bộ lọc FIR có pha tuyến tính và nhân quả.
Để thu được bộ lọc FIR có pha tuyến tính và nhân quả h(n) có độ dài M, ta cần
có
Thao tác này được gọi là "lấy cửa sổ".
Nhận xét:
+ Do cửa sổ w(n) có chiều dài M hữu hạn, đáp ứng của nó có một peaky
main lobe mà độ rộng tỷ lệ với 1/M, và side lobes của nó có chiều cao
thấp hơns.
2
1
,
0
10)(
)(
−
=
jw
)
+ Main lobe sinh ra một dải chuyển tiếp trong H(e
jw
) mà độ rộng là chấp
nhận được so với độ rộng chuyển tiếp. Độ rộng này tỷ lệ với 1/M. (the
wider the main lobe, the wider will be the transition width).
+Side lobes sinh ra các gợn sóng có hình dạng như nhau trong cả dải
thông và dải chắn.
Ý tưởng thiết kế cửa sổ căn bản
+Với các chỉ tiêu bộ lọc đã cho, chọn chiều dài bộ lọc và một hàm cửa
sổ w(n) với độ rông main lobe hẹp nhất và hệ số suy giảm side lobe bé
nhất có thể được.
+Từ nhận xét 4 nêu trên ta chú ý rằng dung sai dải thông delta1 và dung
sai dải chắn delta2 không thể ấn định một cách độc lập. Ta lấy chung
delta
2
= delta
1
.
Hanning Window
This is a raised cosine window function
Hamming Window
−≤≤
−≤≤
−
−
=
otherwise
Mn
M
n
nw
,0
10,
1
2
cos46.054.0
)(
π
14
A. THIẾT KẾ BỘ LỌC THÔNG DẢI PHƯƠNG PHÁP CỬA
SỔ
%thiet ke bo loc bang phuong phap kaiser window
ws1 = 0.2*pi; wp1 = 0.35*pi;
wp2 = 0.65*pi; ws2 = 0.8*pi;
set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,0.65,0.8,1])
set(gca,'YTickMode','manual','YTick',[-60,0])
set(gca,'YTickLabelMode','manual','YTickLabels',['60';' 0'])
15
2. KẾT QUẢ
2 /
*
(0), 0
( ) ( )
( ) 1,2, 1
j k M
H k
H k H e
H M k k M
π
=
= =
− = −
L
16
THIẾT KẾ BỘ LỌC VI PHÂN THEO CẤU TRÚC FIR BẰNG PHƯƠNG
PHÁP LẤY MẪU TẦN SỐ
(FREQUENCY SAMPLING DESIGN TECHNIQUES)
( HOÀNG HỮU VẤN )
Cách tiếp cập này đưa trên sự kiện là hàm truyền hệ thống H(z) có thể thu được từ
các mẫu H(k) của đáp ứng tần số H(e
jw
H z h n z
M z e
π
−
− −
−
−
= =
−
= =
−
∑ ∑
1
2 /
0
1 ( )
( )
1
jwM
M
jw
jw j k M
k
e H k
H e
M e e
π
−
−
−
H k
k
H k H
H k M
M
π
π
−
=
= =
÷
= −
L
1 2 1
, 0, ,
2 2
( )
1 2 ( ) 1
, 1, , 1
2 2
M k M
k
M
1 2 ( ) 1
, 1, , 1
2 2 2
M k M
k
M
H k
M M k M
k M
M
π π
π π
− −
± − =
÷ ÷ ÷
∠ =
− − −
− ± − = + −
÷ ÷ ÷
lỗi xấp xỉ ở tất cả các tần số khác nhau phụ thuộc vào tần số của đáp ứng tần số lý
tưởng, nghĩa là đáp ứng tần số lý tưởng càng sắc nét thì lỗi xấp xỉ càng lớn
Lối càng lớn khi ở gần cạnh dải và càng bé khi ở bên trong dải
Hai cách tiếp cạnh thiết kế:
Phương pháp thiết kế đơn giản: Sử dụng ý tưởng cơ bản và không đưa ra một
ràng buộc nào về lỗi xấp xỉ, nghĩa là chấp nhận lỗi sinh ra do thiết kế
Phương pháp thiết kế tối ưu: cố gắng tối thiểu hoá lỗi trong dải chắn bằng cách
thay đổi các giá trị của mẫu trong dải chuyển tiếp.
Phương pháp thiết kế đơn giản
Trong phương pháp này ta đặt H(k)=H
d
(e
j2pik/M
), k=0,…, M-1 và sử dụng (7.35)
đến (7.39) để thu được đáp ứng xung h(n).
Ví dụ 7.14
Phương pháp thiết kế đơn giản hiếm khi được sử dụng trong thực tế.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-60
-50
-40
-30
-20
-10
0
10
Độ suy hao dải chắn tối thiểu vào khoảng 16dB, rõ ràng là không thể chấp nhận.
Nếu tăng M, sẽ có các mẫu trong dải chuyển tiếp mà chúng ta không thể biết rõ đáp
ứng tần số.
Phương pháp thiết kế tối ưu
2
0.65 , 1
p
w Ap dB
π
= =
1
0.2 , 60
s
w As dB
π
= =
19
2. Chương trình thiết kế bộ lọc thông dải bằng phương pháp lấy mẫu
tần số
2.1. Bài toán thiết kế
Cho các thông số sau đây:
Cạnh thấp dải chắn:
Cạnh thấp dải thông:
Cạnh cao dải thông:
Cạnh cao dải chắn:
Sử dụng phương pháp lấy mẫu tần số để thiết kế bộ lọc FIR thông dải
Giải :
Chúng ta chọn M=40 để chúng ta có 2 mẫu trong dải chuyển tiếp, đưa vào tần số lấy
mẫu trong dải chuyển tiếp là T
1
và T
2
Các mẫu của đáp ứng độ lớn:
{
n = [1:1:L]; n = n-0.5;
w = [0:1:500]'*pi/500;
Hr = cos(w*n)*b';
freqz_m:
function [db,mag,pha,grd,w]=freqz_m(b,a);
% Modified version of freqz subroutine
%
% [db,mag,pha,grd,w]=freqz_m(b,a);
% db=Relative magnitude in dB computed over 0 to pi radians
% mag=absolute magnitude computed over 0 to pi radians
% grd= Group delay over 0 to pi radians
% w=501 frequency samples between 0 to pi radians
% b=numerator polynomial of H(z) (for FIR: a=h)
% a=demonitor polynomial of H(z) (for FIR: a=[1])
%
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))';w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
20
pha=angle(H);
grd=grpdelay(b,a,w);
Chương trình chính:
% THIET KE BO LOC FIR THONG DAI SU DUNG PP LAY MAU TAN SO
% Cac thong so cua bo loc:
% ws1=0.2pi, wp1=0.35pi, wp2=0.65pi, ws2=0.8pi, Rp=1dB, As=60dB
% T2 = 0.59417456, T1=0.109021
M = 40; alpha = (M-1)/2; l = 0:M-1; wl = (2*pi/M)*l;
T1 = 0.109021; T2 = 0.59417456;
Hrs =
22
C. THIẾT KẾ BỘ LỌC SỐ BẰNG PHƯƠNG PHÁP TỐI ƯU:
(NGUYỄN VĂN LƯỢM)
Việc thiết kế bộ lọc thông thấp FIR bằng các sử dụng kỹ thuật thiết kế dùng cửa sổ sẽ
đơn giản và mộ cách tổng quát ,dẫn đến kết quả là một bộ lọc có phaarmr chất tương đối
tốt.Tuy nhiên,các bộ lọc này không tối ưu:
1. Thứ nhất ,các độ lệch dải thông và dỉa chận ,delta_p và delta_s,xấp xỉ bằng
nhau.mặc dù thường ta cần delta_s nỏ hưn nhiều so với delta_p các thông số này không thể
được điều khiển một cách độc lập trong phương pháp thiết kế sử dụng cửa sổ.Do vậy,với
phương pháp thiết kế sử dụng cửa sổ, vta cần phải thiết kế lố dải thông của bộ lọc để thỏa các
yêu cầu nghiêm ngặt hơn trong dải chận.
2. Thứ hai,với đa số các cửa sổ ,độ gợn sẽ không đều hoặc trong dỉa thông hoặc trong
dải chận và một cách tông quát giảm khi di chuyển từ giải chuyển tiếp.việc cho phắp độ gợn
được phân bố đều trên toàn dỉa sẻ tạo ra độ gợn đỉnh nhỏ hơn.
Bộ lọc có pha tuyến tính đồng độ gợn , mặt khác, sẻ tối ưu theo nghĩa biên độ của độ gợn
được tối thiểu hóa trên tất cả các dải đối với bậc cho trước N của bộ lọc.
Trong đề cập sau đây,ta khảo sát việc thiết kế bộ lọc có pha tuyến tính loại 1.kết quả thiết
kế này có thể dể dàng được sửa đổi để thiết kế các loại bộ lọc có pha tuyến tính loại khác.
Dáp ứng bộ lọc FIR có pha tuyến tính có thể viết như sau:
Trong đó các giá trị từ beta và các biểu thức đối với Hr(w) được cho trong bảng
7.2(P.278)
Sử dụng các biểu thức lượng đơn giản,mỗi biểu thức đối với Hr(w) trên có thể được viết
dưới dạng tích của hàm cố định theo w(Q(w)) và một hàm là tổng của cosine(P(w)).
A. Mục đích của phân tích là để có được công thức chung đối với Hr(w) trong cả bốn
trường hợp .giúp cho định dạng bài toán dễ dàng hơn.Để thiêt lập bài toán thiết kế bộ lọc
như một bài toán xấp xỉ Chebyshev,ta phải xác định đáp ứng độ lớn Hdr(w) và một hàm
trọng số W(w), được xác định trên cả dải thông và dải chắn.
B. Bài toán xắp xỉ chepbyseb:
Hàm trọng là cần thiết để có thể điều khiển độc lập trên delta1 và delta2. lỗi trọng số
sprdr
wwSwwHwHwWwE ∪=∈−=
SwwP
wQ
wH
wQwW
wPwQwHwWwE
dr
dr
∈
−=
−=
,)(
)(
)(
)()(
)]()()()[()(
)]()(
ˆ
)[(
ˆ
)( wPwHwWwE
dr
−=
Khi ấy weight=1 mọi nơi, Và ftype không phải là bộ biến đổi Hilbert hoặc bộ Vi phân
[h]=remez(N,f,m)
h là các hệ số bộ lọc có chiều dài M=N+1
N ký hiệu bậc của bộ lọc
f một mảng ký hiệu các cạnh dải tính theo đơn vị của π.
m – đáp ứng biên độ mong muốn ở mỗi f
Swww
wEwEwE
L
S
ii
∈<<<∀
±=±=−=
+
−
110
1
|)(|max)()(
δ
π
δδ
δδ
δ
δ
δ
δ
δ
2
,
f
f
M
A
R
s
p
p
−
=∆
∆
−
=
+=⇒
+
−=
+
−
=⇒
+
−
−=
−
24
a. Thiết kế bộ lọc thông dải dùng phương pháp tối ưu:
ws1 = 0.2*pi; wp1 = 0.35*pi; wp2 = 0.65*pi; ws2 = 0.8*pi;
Rp = 1.0; As = 60;
delta1=(10^(Rp/20)-1)/(10^(Rp/20)+1);
delta2=(1+delta1)*(10^(-As/20));
delataH=max(delta1,delta2);deltaL=min(delta1,delta2);
axis([0,M-1,-0.4,0.5]);ylabel('h(n)');text(M,-0.4,'n')
set(gca,'XTickMode','manual','XTick',[0:3:M-1])
set(gca,'YTickMode','manual','YTick',[-0.4:0.1:0.5])
%dap ung do lon
subplot(2,2,2);plot(w/pi,db);title('dap ung do lon o db');
axis([0,1,-80,10]); xlabel('frequency in pi units'); ylabel('Decibels')
set(gca,'XTickMode','manual','XTick',f)
set(gca,'YTickMode','manual','YTick',[-60,0]);
set(gca,'YTickLabelMode','manual','YTickLabels',['60';' 0']);grid
%dap ung bien do
subplot(2,2,3);plot(omega/pi,Hr);title('Dap ung bien do');
axis([0 1 -0.1 1.1]); xlabel('frequency in pi units'); ylabel('Hr(w)')
set(gca,'XTickMode','manual','XTick',f)
set(gca,'YTickMode','manual','YTick',[0,1]);grid
%trong so
delta_w = 2*pi/1000; sp_edge1 = ws1/delta_w+1; sp_edge2 = ws2/delta_w+1;
subplot(2,2,4);
sb1w = omega(1:1:ws1i)/pi; sb1e = Hr(1:1:ws1i);
pbw = omega(wp1i:wp2i)/pi; pbe = Hr(wp1i:wp2i)-1;
sb2w = omega(ws2i:501)/pi; sb2e = Hr(ws2i:501);
plot(sb1w,sb1e,pbw,pbe*(delta2/delta1),sb2w,sb2e);
title('Loi trong so');
axis([0,1,-deltaL,deltaL]);
xlabel('X 10^(-3)'); ylabel('Hr(w)')
set(gca,'XTickMode','manual','XTick',f)
set(gca,'YTickMode','manual','YTick',[-deltaL,0,deltaL]);
set(gca,'XGrid','on','YGrid','on')
25
b. Kết quả chương trình