Tiểu luận xử lý số tín hiệu THIẾT KẾ BỘ LỌC CHẮN DẢI PP LẤY MẪU TẦN SỐ - Pdf 23

Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
LỜI MỞ ĐẦU
Xử lý tín hiệu số (Digital Signal Processing – DSP) đã trở thành một môn học cơ
sở cho nhiều ngành khoa học, kỹ thuật như: Điện, Điện Tử, Tin học, Viễn thông, Tự động
hoá Xử lý tín hiệu số được ứng dụng rộng rãi trong nhiều lĩnh vực và thiết bị như: CD,
VCD, DVD, camera, scanner, y khoa , trong các hệ thống truyền hình số, thông tin địa
lý, bản đồ số, viễn thông v.v
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 phương pháp đều có những đặc điểm và
ưu khuyết điểm riêng.
Trong khuôn khổ của bài tiểu luận môn học, Tôi xin phép được trình bày bài toán
thiết kế bộ lọc FIR chắn dải – BSF (Band Stop Filte) phương pháp lấy mẫu tần số. Nội
dung tiểu luận được chia thành 5 chương:
Chương 1: Bài toán thiết kế
Chương 2: Cơ sở lý thuyết của bộ lọc FIR
Chương 3: Phương pháp thiết kế
- Thiết kế lọc FIR bằng phương pháp lấy mẫu tần số
Chương 4:Chương trình thiết kế.
Tôi xin trân trọng cảm ơn Thầy giáo TS. Ngô Văn Sỹ đã tận tình hướng dẫn, truyền
đạt những kiến thức quý giá, cung cấp tài liệu tham khảo và chỉ bảo cho tôi phương pháp
làm việc khoa học.
Trong quá trình làm tiểu luận, tuy đã hết sức cố gắng song chắc chắn không tránh
khỏi những sai sót. Rất mong nhận được sự góp ý của Thầy, các anh chị và các bạn học

πω
πω
πω
Trong đó:
P
δ
: Dung sai dải thông
S
δ
: Dung sai dải chắn
PP 21
,
ωω
: Tần số giới hạn dải thông
SS 21
,
ωω
: Tần số giới hạn dải chắn
A
S
: Suy hao dải chắn
R
p
: Độ gợn sóng trong dải thông
Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 2 / 21
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
CHƯƠNG 2: CƠ SỞ LÝ THUYẾT BỘ LỌC FIR
Chương này sẽ trình bày sơ lược cơ sở lý thuyết thiết kế bộ lọc số FIR, cùng với
các phương pháp thiết kế nó, từ đó tìm hiểu thiết kế bộ lọc thông dãi theo cấu trúc FIR
bằng phương pháp cửa sổ và phương pháp lấy mẫu tần số.

DAC
Lọc
(khôiphục)
Th
A
Th ra
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
o Lọc thông thấp LPF.
o Lọc thông cao HPF.
o Lọc thông dãi BPF.
o Lọc chắn dãi BSF.
- Căn cứ vào độ dài của đáp ứng xung h(n) của bộ lọc:
o Đáp ứng xung hửu hạn (lọc FIR).
o Đáp ứng xung vô hạn (lọc IIR).
- Dựa vào tính nhân quả (khả năng thực hiện bộ lọc):
o Bộ lọc lý tưởng.
o Bộ lọc thực tế.
2.2. Bộ lọc FIR
2.2.2. Cấu trúc bộ lọc FIR
Một bộ lọc đáp ứng xung hữu hạn với hàm hệ thống có dạng:


=
−−


=+++=
1
0
1


Mnxbnxbnxbny
N

Đây chính là tích chập tuyến tính của các dãy hữu hạn.
Bậc của bộ lọc là M-1, trong khi chiều dài của bộ lọc là M (bằng với số lượng các
hệ số). Các cấu trúc bộ lọc FIR luôn luôn ổn định, và tương đối đơn giản hơn so với các
cấu trúc bộ lọc IIR. Hơn thế nữa, các bộ lọc FIR có thể được thiết kế để có một đáp ứng
pha tuyến tính và đó là điều cần thiết trong một số ứng dụng.
Chúng ta sẽ xem xét lần lượt các cấu trúc của bộ lọc FIR sau đây:
2.2.2.1 Cấu trúc dạng trực tiếp
Phương trình sai phân được thực hiện bởi một dãy liên tiếp các bộ trễ do không có
đường phản hồi:
)1()1()()(
110
+−++−+=

Mnxbnxbnxbny
N


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 2.2 với M = 5:
Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 4 / 21
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
2.2.2.2 Cấu trúc dạng ghép tầng:
Hàm hệ thống H(z) được biến đổi thành các tích của các khâu bậc 2 với các hệ số
thực. Các khâu này được thực hiện ở dạng trực tiếp và bộ lọc tổng thể có dạng ghép tầng
của các khâu bậc 2.


1
1
1
10
1)( 


=
−−
++=
K
1k
2
2,k
1
1,k0
)zBzB1(b
Trong đó






=
2
M
K
, B
k,1

1
b
1
z
-
1
b
2
z
-
1
b
3
z
-
1
b
4
y(n)
x(n)
Hình 2.2: Cấu trúc lọc FIR dạng trực tiếp
B
1,1
z
-
1
z
-
1
z

xứng trong phương trình (1.9), ta có:
)1()2()1()()(
0110
+−++−++−+= MnxbMnxbnxbnxbny 

++−+−++−+= )]2()1([)]1()([
10
MnxnxbMnxnxb
Sơ đồ khối thực hiện phương trình sai phân trên được mô tả trong hình 2.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
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.
2.2.2. Các đặc tính của bộ lọc FIR pha tuyến tính
Trong phần này chúng ta sẽ thảo luận về hình dạng của đáp ứng xung, đáp ứng tần số
trong hàm hệ thống của các bộ lọc FIR pha tuyến tính.
Cho h(n), trong đó 0 ≤ n ≤ M – 1, là đáp ứng xung có chiều dài N thì hàm truyền hệ
thống là:
∑∑

=
−−−−

=

==
1
0
1)1(
1

1
z
-
1
y(n)
b
0
z
-
1
b
1
z
-
1
b
2
b
3
y(n)
x(n)
z
-
1
z
-
1
z
-
1

,10),1()(

=−≤≤−−=
M
MnnMhnh
α
Do đó h(n) là đối xứng theo α, là chỉ số đối xứng. Có hai kiểu đối xứng:
• M lẻ: Trong trường hợp này,
2
1−
=
M
α
là một số nguyên. Đáp ứng xung ứng
xung được mô tả bằng hình 2.6 dưới đây:
• M chẵn:Trong trường hợp này,
2
1−
=
M
α
là một số nguyên. Đáp ứng xung ứng
xung được mô tả bằng hình 2.7 dưới đây:
Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 7 / 21
Hình 2.7: Đáp ứng xung đối xứng, M chẵn
Hình 2.6 Đáp ứng xung đối xứng M lẻ
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
Ta cũng có bộ lọc FIR pha tuyến tính loại hai nếu ta yêu cầu đáp ứng pha
( )
ω

±=β


Đ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 −

. Một lần nữa chúng ta lại có 2 kiểu, cho M lẻ và M chẵn.
• M lẻ: Trong trường hợp này,
2
1M −

là một số nguyên. Đáp ứng xung được
mô tả bằng hình 2.8 dưới đây:
Lưu ý rằng mẫu h(α) tại
2
1M −

phải bằng 0, nghĩa là,
0
2
1M
h =







π
±=β=
αω−βωω
Trong đó H
r
(e
j
ω
) là hàm đáp ứng độ lớn chứ không phải là hàm đáp ứng biên độ.
Đáp ứng độ lớn là một hàm thực, có thể vừa dương vừa âm, không giống đáp ứng biên độ
luôn luôn dương. Đáp ứng pha kết hợp với đáp ứng biên độ là một hàm không liên tục,
trong khi kết hợp với đáp ứng độ lớn là một hàm tuyến tính liên tục.
• Bộ lọc FIR pha tuyến tính Loại-1 (Type 1): Đáp ứng xung đối xứng, M lẻ
Trong trường hợp này
0=β
,
2
1M −

là một biến nguyên, và
( ) ( )
n1Mhnh −−=
,
1Mn0
−≤≤
, thì ta có thể chứng tỏ rằng:
( )
( )
( )
2/1Mj

( )








= n
2
1M
h2na
với
2
3M
n1

≤≤
• Bộ lọc FIR pha tuyến tính Loại-2 (Type 2): Đáp ứng xung đối xứng, M chẵn
Trong trường hợp này
0=β
,
( ) ( )
n1Mhnh −−=
,
1Mn0 −≤≤
, nhưng
2
1M −





−ω=

Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 9 / 21
Hình 2.9 Đáp ứng xung phản đối xứng, M
chẵn
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
Trong đó:
( )






−= n
2
M
h2nb
với
2
M
, ,2,1n =
So sánh hai công thức trên ta có:
( )

=

=












−π=π

=
mà không cần quan
tâm đến b(n) hoặc h(n). Do đó chúng ta không thể sử dụng loại này (h(n) đối xứng, M
chẵn) đối với bộ lọc thông cao hoặc bộ lọc chắn dải.
• Lọc FIR pha tuyến tính Loại-3 (Type 3): Đáp ứng xung phản đối xứng, M lẻ
Trong trường hợp này ta có
2
π

,
2
1M −

là một biến nguyên,
( ) ( )




π

=
ω






ω=

2
1M
2
j
2/1M
0n
j
ensinnc)e(H
Trong đó
( )






r

mà không cần quan tâm c(n) hoặc
h(n). Hơn thế nữa,
je
2
j
=
π
, điều đó có nghĩa là
( )
ω
r
jH
là thuần ảo. Do đó, loại bộ lọc này
không thích hợp đối với việc thiết kế bộ lọc thông thấp hoặc thông cao. Tuy nhiên, điều
này thích hợp đối với việc xấp xỉ các bộ vi phân và bộ biến đổi Hilbert số lý tưởng.
• Lọc FIR pha tuyến tính Loại-4(Type 4):Đáp ứng xung phản đối xứng, M chẵn
Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 10 / 21
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
Trong trường hợp này
2
π

,
( ) ( )
n1Mhnh −−−=
,
1Mn0
−≤≤













−ω=

2
1M
2
j
2/M
1n
j
e
2
1
nsinnd)e(H
Trong đó:
( )




r
2
1
nsinnd)(H
Lưu ý: Tại
π=ω
,
0)0(H
r
=

je
2
j
=
π
. Do vậy, loại này cũng thích hợp cho việc
thiết kế các bộ vi phân số và bộ biến đổi Hilbert số.
Bảng sau đây mô tả khả năng thích hợp trong việc thiết kế các bộ lọc và các bộ biến
đổi Hilbert số, bộ vi phân số của 4 loại lọc FIR pha tuyến tính đã nêu:
Type LPF HPF BPF SBF Hilbert Differentiator
FIR Type 1
   
FIR Type 2
 
FIR Type 3
  
FIR Type 4
   
2.3. Khái niệm lọc FIR chắn dải

ω

2
ω
1
ω
2
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
H(e
j
ω
) : đáp ứng biên
ω
c1 ,
ω
c2
: tần số cắt
-Đáp ứng xung của bộ lọc lý tưởng
2.3.2.Bộ lọc chắn dải thực tế
Các bộ lọc chắn dải thực tế được đặc trưng bởi các thông số kỹ thuật trong miền tần số liên tục ω
Đồ thị đáp ứng biên độ:
Để đặc trưng cho bộ lọc thực tế, người ta sử dụng các tham số sau:
1. Loại bộ lọc: thông thấp, thông cao, thông dải, dải chặn
2. Tần số giới hạn dải thông ω
p1
, ω
p2
Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 12 / 21
h
d

2
1
h
d
(n) =
n
n
n
n
n
n
c
cc
c
cc
1
11
2
22
sinsinsin
ω
ω
π
ω
ω
ω
π
ω
π
π

s1
- ω
p1
,∆ω
2
= ω
p2
- ω
s25. Độ gợn sóng dải thông: δ
1
. Trong dải thông, đặc tính biên độ tần số
( )
ω
j
eH
phải
thoả mãn điều kiện:
( )
( )
( )
11
11
δδ
ω
+≤Η≤−
j
e

ω
j
eH
dạng hình chữ nhật nên không thể thực
hiện được trên thực tế. Cơ sở của phương pháp lấy mẫu tần số là xấp xỉ đặc tính tần số
)(
ω
j
N
eH
của bộ lọc số cần tổng hợp theo đặc tính tần số
)(
ω
j
eH
của bộ lọc số lý tưởng
cùng loại.
Việc xấp xỉ được thực hiện bằng cách lấy mẫu tần số qua DFT, tức là làm cho các
mẫu của
)(
ω
j
N
eH

)(
ω
j
eH
bằng nhau tại các tần số rời rạc ω

, sai số xấp xỉ giữa
)(
ω
j
N
eH

)(
ω
j
eH
bằng 0, còn tại các tần số ở giữa khoảng kω
1
và (k+1)ω
1
thì sai số
xấp xỉ là hữu hạn. Sai số xấp xỉ sẽ giảm nhỏ nếu giảm tần số lấy mẫu cơ bản ω
1
= (2π /
N), điều đó tương ứng với tăng độ dài N của đặc tính xung h(n)
N
của bộ lọc số được tổng
hợp.
Ta có thể xác định được đặc tính biên độ tần số
)(
ω
j
N
eH
của 4 loại bộ lọc số FIR

kj
N
πω
ω
ω
(1)
Đặc tính pha θ(ω) của bộ lọc FIR pha tuyến tính loại 1 và loại 2:
αωωωϕ
−=







−= .
2
1
)(
N
N
, với
ωα
.
2
1




β
=

ωα
.
2
1







=
N
Để đánh giá sai số khi xấp xỉ đặc tính tần số
)(
ω
j
N
eH
theo đặc tính tần số bộ lọc
lý tưởng
)(
ω
j
eH
, người ta dung hàm sai số:
)()()(

tính cần tổng hợp bằng biểu thức nội suy sau:


=







−=
1
0
2
sin
)(
)1(
N
)sin(0.5N
)(
N
k
N
kj
N
N
k
kA
eH

Bước 3: Kiểm tra đặc tính biên độ tần số
)(
ω
j
N
eH
có đạt các tiêu chí kỹ thuật đã
cho δ1, δ2, ω
c
, Δω hay không? Hoặc
p
jj
eE )E(e )(
max
ωω

Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 15 / 21
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
Nếu đạt tất cả các tiêu chí kỹ thuật đã cho thì giảm số điểm lấy mẫu N và thực
hiện lại các bước trên cho đến khi chọn được N
min
đảm bảo đạt tất cả các chỉ tiêu kỹ thuật
đã cho.
Nếu không đạt thì tăng số điểm lấy mẫu N và thực hiện lại các bước trên cho đến
khi chọn được N
min
để
)(
ω
j

2
)0(
)(
N
k
N
k
N
N
n
N
k
kA
NN
A
nh
π
- Đối với bộ lọc loại 2: để tìm h(n)
N
có thể tính IDFT theo công thức:


=






+−+=






+−=
2
1
1
1
3
)12(sin.)()1(
2
)(
N
k
N
k
N
n
N
k
kA
N
nh
π
- Đối với bộ lọc loại 4: để tìm h(n)
N
có thể tính IDFT theo công thức:


k
N
n
N
n
N
k
kA
N
N
A
N
nh
π
Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 16 / 21
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
CHƯƠNG 4: CHƯƠNG TRÌNH THIẾT KẾ
BÀI TOÁN :
4.1. Thiết kế bằng phương pháp lấy mẫu tần số:
Các bước tính toán thiết kế lọc FIR chắn dải với các thông số sau:
dBA
dbR
dBR
dBA
Ss
pp
pp
SS
40;6,0
5,0;7,0

πππωωω
πππωωω
65,02/6,07,02/
35,02/)4,03,0(2/
222
111
=+=+=
=+=+=
SPC
SPC
Chúng ta chọn M=63 do đó chúng ta có mẫu tần số tại
p1
ω
đó là tại k = 9:
9
63
2
3.0
1
π
πω
==
p
Và mẫu tiếp theo tại
s1
ω
đó là tại k= 13:
13
63
2

==
p
Do đó chúng ta có 9 mẫu trong stopband [
p1
0
ωω
≤≤
] và 6 mẫu trong passband [
ss 21
ωωω
≤≤
] và 9 mẫu trong stopband [
πωω
≤≤
p2
]. Chúng ta có:
( )
KH
r
=[ones(1,9),zeros(1,13),ones(1,19),zeros(1,13),ones(1,9)]
Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 17 / 21
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
Khi chọn M=63,
31
2
163
=

=
α

% [Hr,w,a,L] = Hr_Type1(h)
% Hr = Amplitude Response
% w = 500 frequencies between [0 pi] over which Hr is
computed
% a = Type-1 LP filter coefficients
% L = Order of Hr
% h = Type-1 LP filter impulse response
M = length(h);
L = (M-1)/2;
a = [h(L+1) 2*h(L:-1:1)]; % 1x(L+1) row vector
n = [0:1:L]; % (L+1)x1 column vector
w = [0:1:500]'*pi/500;
Hr = cos(w*n)*a';
 freqz_m(b,a): Tính đáp ứng tần số của bộ lọc FIR có đáp ứng xung h.
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));
pha=angle(H);

subplot(2,2,4);plot(w/pi,db); axis([0,1,-100,10]); grid
title('Magnitude Response '); xlabel('Frequency in pi units');
ylabel('dB');
set(gca,'XTickMode','manual','XTick',[0,0.3,0.4,0.6,0.7,1])
set(gca,'YTickMode','Manual','YTick',[-16;0]);
set(gca,'YTickLabelMode','manual','YTickLabels',['16';' 0'])
 Kết quả chạy chương trình
 Theo ý tưởng cơ bản :
Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 19 / 21
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
Hình 4.3: Mô tả đáp ứng sau khi thiết kế
Nhận xét:
- Lỗi xấp xỉ là hiệu của đáp ứng lý tưởng và đáp ứng thực tế bằng không tại các tần
số được lấy mẫu.
- Lỗi xấp xỉ ở tất cả các tần số khác nhau phụ thuộc vào hình dạng 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.
Quan sát độ thị hình 4.3 thì độ suy hao dải chắn tối thiểu khoảng 42dB thỏa mãn
với giá trị As=40dB đã cho. Do vậy phương pháp xứ lý này chấp nhận được.
TÀI LIỆU THAM KHẢO
1. Hồ Văn Sung, Xử lý số tín hiệu, NXB Giáo Dục,2003
Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 20 / 21
Tiểu luận xử lý số tín hiệu: THIẾT KẾ BỘ LỌC CHẮN DẢI - PP LẤY MẪU TẦN SỐ
2. Lý thuyết và bài tập xử lý tín hiệu số-Tống Văn On, NXB Lao Động - Xã Hội,2002
3.Xử lý tín hiệu và lọc số-Nguyễn Quốc Trung, NXB Khoa Học và Kỹ Thuật,2002
4.Xử lý tín hiệu số-Quách Tuấn Ngọc, NXB Giáo Dục,1995
5.Digital Signal Processing Using Matlab V4-Vinay K.Ingle-John G.Proakis.
Học Viên: Vũ Bảo Toàn – Lớp: Cao học Tự Động Khóa _K24 Trang 21 / 21


Nhờ tải bản gốc

Tài liệu, ebook tham khảo khác

Music ♫

Copyright: Tài liệu đại học © DMCA.com Protection Status