Xử lý tín hiệu nâng cao - Chương 4: Biến đổi Fourier của tín hiệu rời

Xử lý tín hiệu nâng cao
-Advanced signal processing-
Chương 4
Biến đổi Fourier của tín hiệu rời rạc
Biến đổi Fourier của tín hiệu rời rạc
X Y
T
Miền không 
gian ban đầu
Không gian 
đặc trưng
T-1
Định nghĩa

Biến đổi Fourier của tín hiệu rời rạc được định 
nghĩa như sau:
∑
+∞
−
=
njj enxeX ωω )()(

Toán tử FT:
−∞=n
( )[ ] ( )ωjeXnxFT =
Biến đổi Fourier ngược

 Từ miền tần số tín hiệu cũng có thể biến đổi ngược lại miền 
thời gian bằng phép biến đổi Fourier ngược:
∫
+
−
=
pi
pi
ωω ω
pi
deeXnx njj )(
2
1)(

 Ta sử dụng ký hiệu IFT để biểu diễn biến đổi Fourier 
ngược:
( )[ ] ( )nxeXIFT j =ω
Các phương pháp thể hiện của X(ejω)

Thể hiện dưới dạng phần thực và phần ảo:
( ) Re ( ) Im ( )j j jX e X e j X eω ω ω   = +   
Các phương pháp thể hiện của X(ejω)

Thể hiện dưới dạng module và argument:
[ ])(arg)()( ωωω jeXjjj eeXeX =

Khi đó:
 X(ejω) : Phổ của tín hiệu x(n).
 |X(ejω)| : phổ biên độ của x(n)
 arg[X(ejω)]= ϕ(ω): phổ pha của x(n)
Các phương pháp thể hiện của X(ejω)

Quan hệ giữa phổ pha và phổ biên độ với thành 
phần thực và ảo của X(ejω):
 Phổ biên độ:
[ ] [ ])(Im)(Re)( 22j ωωω jj eXeXeX +=
 Phổ pha:
[ ] [ ][ ])(Re
)(Im)(arg
ω
ω
ω
j
j
j
eX
eX
arctgeX =
Tính chất quan trọng của X(ejω)

Tuần hoàn: Biến đổi Fourier của tín hiệu 
X(ejω) tuần hoàn với chu kỳ 2pi. 

Tính đối xứng:
j jω ω−   =Re ( ) Re ( )
Im ( ) Im ( )
( ) ( )
( ) ( )
j j
j j
j j
X e X e
X e X e
X e X e
X e X e
ω ω
ω ω
ω ω
−
−
−
   
   = −   
=
   ∠ = −∠   
Ví dụ 1

Cho dãy 
 Viết chương trình bằng MATLAB thể hiện trên đồ 
thị phổ tại 500 điểm rời rạc trong khoảng [- pi, pi].
( ) 0,5 ( )nx n u n=
Ví dụ 1

Thực hiện biến đổi Fourier của tín hiệu:
 Áp dụng công thức, sẽ có:
)(5.0)( nunx n=
∞+∞
∑
∑∑
∞
−
−
−
∞−
−
−
==
==
0
0
5.01
1)5.0(
5.0)()(
ω
ω
ωωω
j
nj
njnnjj
e
e
eenxeX
Ví dụ 1

a. Tìm biến đổi Z của tín hiệu x(n) bằng hàm ztrans
syms n positive;
x=0.5.^n;
X=ztrans(x)
X =
z/(z - 1/2)
 Khi đó phổ của tín hiệu x(n) là: 
( )
0.5
j
j
j
eX e
e
ω
ω
ω
=
−
Ví dụ 1

b. Viết chương trình bằng MATLAB thể hiện trên đồ 
thị phổ tại 500 điểm rời rạc trong khoảng [-pi, pi].
w=linspace(-pi,pi,500);
X = exp(j*w) ./ (exp(j*w)- 0.5*ones(1,500));
magX = abs(X); angX = angle(X);
realX = real(X); imagX = imag(X);
subplot(2,2,1); plot(w,magX); grid;
title('Pho bien do'); xlabel('Tan so chuan hoa (pi)');
ylabel('Bien do');
Ví dụ 1
subplot(2,2,3); plot(w,angX); grid;
title('Pho pha'); xlabel('Tan so chuan hoa (pi)');
ylabel('Pho theo radians');
subplot(2,2,2); plot(w,realX); grid;
title('Phan thuc'); xlabel('Tan so chuan hoa (pi)');
ylabel('Bien do');
subplot(2,2,4); plot(w,imagX); grid;
title('Phan ao'); xlabel('Tan so chuan hoa (pi)');
ylabel('Bien do');
Ví dụ 1
-4 -2 0 2 4
0.5
1
1.5
2
Pho bien do
B
i
e
n
d
o
-4 -2 0 2 4
0.5
1
1.5
2
Phan thuc
B
i
e
n
d
o
Tan so chuan hoa (pi)
-4 -2 0 2 4
-1
-0.5
0
0.5
1
Pho pha
Tan so chuan hoa (pi)
P
h
o
t
h
e
o
r
a
d
i
a
n
s
Tan so chuan hoa (pi)
-4 -2 0 2 4
-1
-0.5
0
0.5
1
Phan ao
Tan so chuan hoa (pi)
B
i
e
n
d
o
freqz

Trong Matlab còn có hàm freqz trả về đáp 
ứng tần số của một hệ thống tại một số hữu 
hạn các điểm rời rạc trên vòng tròn đơn vị khi 
biết hàm truyền đạt của nó

Viết chương trình Matlab sử dụng hàm freqz 
để vẽ đồ thị phổ của tín hiệu 
( ) 0,5 ( )nx n u n=
Bài tập

Cho dãy 
 Viết chương trình tính và thể hiện phổ của dãy x(n) 
tại 512 điểm rời rạc trong khoảng [0,pi]
8( ) ( )x n rect n=
Bài tập

Cho phổ X(ejω) có dạng sau: 
 Viết chương trình thể hiện trên đồ thị các hàm phổ
biên độ, phổ pha, phần thực và phần ảo của
/ 2( ) .sin(3 )j jX e eω ω ω−=
X(ejω), tính tại 2001 điểm rời rạc trong khoảng
[-2pi,2pi].
Tín hiệu và hệ thống trong miền tần số rời rạc

Biến đổi Fourier rời rạc thuận DFT
 Cho dãy x(n) có chiều dài hữu hạn, khi đó biến đổi Fourier 
rời rạc thuận được định nghĩa:
( ) ( )
1
0
0 1
0
N
kn
N
n
x n W k N
X k
k
−
=
 ≤ ≤ −
= 
 ≠
∑
với:

2
k
j knj nkn N
NW e e
pi
ω −−
= =
Tín hiệu và hệ thống trong miền tần số rời rạc

Biến đổi Fourier rời rạc thuận DFT
 Ta cũng có thể biểu diễn DFT dưới dạng ma trận:
 Ta khai triển:
( ) ( )1
0
.
N
kn
N
n
X k x n W
−
=
=∑
( ) ( ) ( ) ( ) ( )0 0 0 0= + + + + −
• k = 0:
• k = 1:
• k = 2:
• k = N-1:
0 0 1 2 1N N N NX x W x W x W x N W
( ) ( ) ( ) ( ) ( ) ( )10 1 21 0 1 2 1 NN N N NX x W x W x W x N W −= + + + + −

( ) ( ) ( ) ( ) ( ) ( )2 10 2 42 0 1 2 1 NN N N NX x W x W x W x N W −= + + + + −

( ) ( ) ( ) ( ) ( ) ( ) ( )( )2 1 1 10 11 0 1 2 1N N NNN N N NX N x W x W x W x N W− − −−− = + + + + −

Tín hiệu và hệ thống trong miền tần số rời rạc

 Ký hiệu:
( )
( )
( )
( )
( )
0
1
2
1
X
X
X k X
X N
 
 
 
 =
 
 
 
− 

( )
( )
( )
( )
( )
0
1
2
1
x
x
x n x
x N
 
 
 
 =
 
 
 
− 

( )
( )
( ) ( ) ( )( )
0 0 0 0
10 1 2
2 10 2 4
1 2 1 1 10
N N N N
N
N N N N
N
N N N N N
N N N N
W W W W
W W W W
W W W W W
W W W W
−
−
− − − −
 
 
 
 
=  
 
 
  
    

 Ta có thể viết lại cho gọn dạng biểu diễn theo ma trận như 
sau:
N N N N
( ) ( ). NX k x n W=
Biến đổi Fourier rời rạc thuận DFT
function [Xk] = dft(xn,N)
% Tim bien doi Fourier roi rac thuan
% ---------------------------------------------
% [Xk] = dft(xn,N)
% Xk = day cac he so DFT tren doan 0<=k<=N-1
% xn = day huu han N diem
% N = chieu dai DFT
%
n = [0:1:N-1];
k = [0:1:N-1];
WN = exp(-j*2*pi/N);
nk = n' * k;
WNnk = WN .^ nk; % ma tran DFT
Xk = xn * WNnk;
Biến đổi Fourier rời rạc thuận DFT

 Dựa vào hàm này, có thể tính DFT 20 điểm của tín hiệu 5( ) ( )x n rect n=
% Tinh DFT 20 diem cua day x(n)
L = 5; N = 20;
n = [0:N-1];
xn = [ones(1,L), zeros(1,N-L)];
k = n;
Xk = dft(xn,N);
magXk = abs(Xk);
subplot(4,2,1); stem(n,xn);
axis([min(n),max(n)+1,-0.5,1.5]);
title('Sequence x(n)');
xlabel('n'); ylabel('x(n)');
subplot(4,2,3); stem(k,magXk);
axis([min(k),max(k)+1,-0.5,5.5]);
title('DFT of SQ. wave: L=5, N=20');
xlabel('k'); ylabel('X(k)');
Ví dụ 2

Thực hiện biến đổi Fourier của tín hiệu: 
x(n)={1,2↑,3,4,5} với n=[-1:3]
 Áp dụng công thức, có:
2 3
( ) ( )
2 3 4 5
j j n
j j j j
X e x n e
e e e e
ω ω
ω ω ω ω
+∞
−
−∞
− − −
=
= + + + +
∑
Ví dụ 2 (tiếp)

Trên Matlab:
n=-1:3; x=1:5; k=0:500; w=(pi/500)*k;
X=x*(exp(-j*pi*(n'*k)/500));
subplot(2,2,1); plot(k,abs(X));
title('Bien do'); grid;
subplot(2,2,2); plot(k,real(X));
title('Phan thuc'); grid;
subplot(2,2,3); plot(k,imag(X));
title('Phan ao'); grid;
subplot(2,2,4); plot(k,angle(X));
title('Pha'); grid;
Ví dụ 2 (tiếp)
0 200 400 600
0
5
10
15
Bien do
0 200 400 600
-5
0
5
10
15
Phan thuc
0 200 400 600
-10
-5
0
5
Phan ao
0 200 400 600
-4
-2
0
2
4
Pha
Bài tập

Viết chương trình tính và thể hiện trên đồ thị 
biến đổi Fourier rời rạc của các dãy sau:
 Dãy có chiều dài 40 và 
 Dãy có chiều dài 80 và 
5( ) ( )x n rect n=
( ) ( )x n rect n=
 Dãy có chiều dài 100 và 
5
8( ) ( )x n rect n=
Biến đổi Fourier rời rạc ngược IDFT
( ) ( )
1
0
1 0 1
0
N
kn
N
k
X k W n N
x n N
n
−
−
=
 ≤ ≤ −
= 
 ≠
∑
với:
2
k
j knj nkn N
NW e e
pi
ω −−−
= =
Biến đổi Fourier rời rạc ngược IDFT
function [xn] = idft(Xk,N)
% Tim bien doi Fourier roi rac nguoc
% ---------------------------------------------
% [xn] = idft(Xk,N)
% xn = day co chieu dai huu han tren doan 0<=n<=N-1
% Xk = day cac he so DFT tren doan 0<=k<=N-1
% N = chieu dai DFT
%
n = [0:1:N-1];
k = [0:1:N-1];
WN = exp(-j*2*pi/N);
nk = n' * k;
WNnk = WN .^ (-nk); % ma tran IDFT
xn = (Xk * WNnk)/N;
Các tính chất của biến đổi Fourier

Tuyến tính: Giả sử ta có hai tín hiệu x1(n) và x2(n) và 
biến đổi Fourier tương ứng là:
FT[x1(n)]=X1(ejω)
FT[x2(n)]=X2(ejω)

Khi đó 1 tín hiệu là tổ hợp tuyến tính của tín hiệu x1 và x2: 
x(n)=a*x1(n)+b*x2(n) và biến đổi Fourier của x(n) là X(ejω) 
thì:
X(ejω)=a* X1(ejω)+b* X2(ejω)
Các tính chất của biến đổi Fourier

Tính chất trễ:
 Ta có:
( )[ ] ( )ωjeXnxFT =
 Khi đó:
( )[ ] )(00 ωω jnj eXennxFT −=−
Các tính chất của biến đổi Fourier

Trễ tần số:
 Ta có:
( )[ ] ( )ωjeXnxFT =
 Khi đó:
( ) ( )( )0o jj nFT x n e X e ω ωω −  = 
Các tính chất của biến đổi Fourier

Liên hợp phức:
( )[ ] ( )ωjeXnxF −= **
Nhân chập và tích đại số

Nhân chập:
FT[x1(n)⊗ x2(n)]=X1(ejω)* X2(ejω)

Tích đại số:
FT[x1(n). x2(n)]=X1(ejω) ⊗ X2(ejω)
Biến đổi Fourier nhanh FFT

 Trong biểu thức DFT ta thấy có N phương trình, trong mỗi
phương trình có N phép nhân. Do đó, để tính DFT cần N2
phép nhân. Thuật toán biến đổi Fourier nhanh FFT cho phép
ta khắc phục được nhược điểm này, nghĩa là cho phép giảm
số phép nhân xuống khi tính DFT.

 Chương trình dưới đây biểu diễn trên đồ thị biểu đồ thể hiện
mối quan hệ giữa chiều dài dãy N, N biến thiên từ 1 đến 2048,
với thời gian thực hiện biến đổi Fourier của hàm MATLAB.
Biến đổi Fourier nhanh FFT
% Do thi thoi gian tinh FFT
Nmin = 1;
Nmax = 2048;
fft_time = zeros(1,Nmax-Nmin+1);
for n = Nmin:1:Nmax
x = rand(1,n);
t = clock;
fft(x);
fft_time(n-Nmin+1) = etime(clock,t);
end %for
n = [Nmin:1:Nmax];
plot(n,fft_time,'.')
xlabel('N');ylabel('Time in Secs');
title('FFT execution times');
Biến đổi Fourier nhanh

Trong Matlab hàm fft để tính Fourier nhanh:
n=-1:3;
x=1:5;
k=0:500;
X=x*exp(-j*2*pi*(n'*k)/500);
X1=fft(x,501);
subplot(221);plot(2*k/500,abs(X));
subplot(222);plot(2*k/500,abs(X1));
subplot(223);plot(2*k/500,angle(X));
subplot(224);plot(2*k/500,angle(X1));
Biểu diễn hệ thống trong miền tần số liên 
tục

Đáp ứng tần số:
[ ] njjnjnjkj
knjnj
eeHenhFTeekh
ekhenhny
ωωωωω
ωω
)()()(
)()()( )(
==





=
=⊗=
∑
∑
∞+
∞−
−
+∞
∞−
−
( )[ ] ∑
+∞
−∞=
−
==
n
njj enhnhFTeH ωω )()(
[ ] ∫
+
−
==
pi
pi
ωωω ω
pi
deeHeHIFTnh njjj )(
2
1)()(
Biểu diễn H(ejω)
 H(ejω) là hàm biến số phức:
( ) Re ( ) Im ( )j j jH e H e j H eω ω ω   = +   
[ ])(arg)()( ωωω jeHjjj eeHeH =
Biểu diễn H(ejω)

Ta cũng có quan hệ giữa đáp ứng tần số và đáp
ứng pha của hệ thống với phần thực và phần ảo
của H(ejω)
 Phổ biên độ:
 Phổ pha:
[ ] [ ])(Im)(Re)( 22j ωωω jj eHeHeH +=
[ ] [ ][ ])(Re
)(Im)(arg
ω
ω
ω
j
j
j
eH
eH
arctgeH =
Công thức quan trọng
)(nh)(nx )(*)()( nxnhny =
( ), ( )jH e h nω( )jX e ω ( ) ( ) ( )j j jY e H e X eω ω ω=
Các bộ lọc số lý tưởng

Bộ lọc thông thấp lý tưởng

Bộ lọc thông cao lý tưởng

Bộ lọc thông dải lý tưởng

Bộ lọc chắn dải lý tưởng
Bộ lọc thông thấp lý tưởng

Bộ lọc thông thấp được định nghĩa bằng công thức:
Với –pi ≤ ω ≤ pi
otherwise
-
0
1)( cc ωωωω ≤≤



=
jeH
-ωc ωc pi-pi
1
ω
|H(ejω)|
0
Bộ lọc thông thấp lý tưởng

Xét thông thấp lý tưởng có tần số cắt ωc=pi/3

Áp dụng công thức ta có:
pi
pi
.
)sin()( 3
n
n
nh =
...
...
0,33
h(n)
0,280,28
0,14
-0,07
0,14
-0,07
0,04
0,04 n0,030,03
-0,05
-0,05
Bộ lọc thông cao lý tưởng

Bộ lọc thông cao được định nghĩa bằng công thức:
Với –pi ≤ ω ≤ pi
otherwise
and
eH j
piωωωωpiω ≤≤−≤≤



=
cc-
0
1)(
-ωc ωc pi-pi
1
ω
|H(ejω)|
0
Bộ lọc thông dải lý tưởng

Bộ lọc thông dải được định nghĩa bằng công thức:
Với –pi ≤ ω ≤ pi
otherwise
and
eH j c2c1c1c2
`-
0
1)( ωωωωωωω ≤≤−≤≤



=
-ωc2 ωc2 pi-pi
1
ω
|H(ejω)|
0
-ωc1 ωc1
Biến đổi Fourier của 
tín hiệu rời rạc hai chiều
Khái niệm và công thức

Phép biến đổi Fourier của tín hiệu rời rạc 2 chiều 
(ảnh số) được tính bằng công thức: 
∑∑
−
=
−
=






+pi−
=
1M
0m
1N
0n
N
vn
M
umj2
e)n,m(x)v,u(X

Ánh xạ ngược của phép biến đổi
∑∑
−
=
−
=






+pi
=
1M
0u
1N
0v
N
vn
M
umj2
e)v,u(X
MN
1)n,m(x
Khái niệm và công thức

Các thành phần tần số mang giá trị phức nên 
ta có thể biểu diễn như sau:
),arg(),(),( vujevuXvuX −=
Khi đó |X(u,v)| được gọi là độ lớn hay phổ biên độ, arg(u,v) được gọi 
là phổ pha
Một số tính chất

Tính tuần hoàn

Đối xứng và đơn vị
)Nv,Mu(F)Nv,u(F)v,Mu(F)v,u(F ++=+=+=
 FT(u,v)=F(u,v)
 F-1(u,v)=F*(u,v)
Một số tính chất

Tính chất chuyển đổi
∑∑
− − 




 −
+
−
−
=−−
1 1 )()(2
),(),(
M N
N
nbv
M
mauj
enmxbvauX
pi
∑∑
−
=
−
=






+





+−
= =
=
1
0
1
0
22
0 0
),(
M
m
N
n
N
bn
M
amj
N
vn
M
umj
m n
eenmx
pipi
Tính chất chuyển đổi (tiếp)

Nhân tín hiệu với e2jpi(am/M+bn/N) trong miền không gian 
thực sẽ tương đương với dịch chuyển phổ đi một 
khoảng (a,b). Xét trường hợp đặc biệt khi a=M/2, 
b=N/2
∑∑
− − 



+pi−1M 1N vnumj2NM

Nhân vào ảnh ban đầu giá trị (-1)(m+n) trước khi biến 
đổi, ta sẽ thu được phổ tần số mà điểm tần số F(0,0) 
của nó sẽ nằm giữa mảng 2 chiều. 
= =
+
−=−−
0m 0n
)nm(NM )1(e)n,m(x)
2
v,
2
u(X
Một số tính chất

Tích chập
 Ta có
DFT(x(m,n))=X(u,v)
DFT(h(k,l))=H(u,v)
 Khi đó:
DFT(x(m,n)*h(k,l))=X(u,v)H(u,v)
Phép lọc trên miền tần số
x(m,n) y(m,n)
DFT filter IDFT
F(u,v)
H(u,v)F(u,v)
Ví dụ: phép lọc thông thấp

Bộ lọc thông thấp



>
≤
=
0
0
Dv)D(u, if 0
Dv)D(u, if 1)v,u(H

Khoảng cách tới nguồn
22
2
N
v
2
M
u)v,u(D 





−+





−=
Ví dụ: phép lọc thông thấp
Trong Matlab

Sử dụng hàm fft2
I=imread('cameraman.tif');
F=fft2(I);
imshow(abs(F),[]);

Sử dụng lệnh fftshift để dịch phổ về tâm 
FC=fftshift(abs(F))
imshow(FC,[])

Để hiện thị phổ được rõ hơn sử dụng thêm hàm log
FC2=log(1+FC);
imshow(FC2,[]);
Kết quả
Thực hành chương III

Thực hiện biến đổi Fourier, biển diễn phổ biên 
độ và phổ pha của các tín hiệu:
 x(n)=2(0.8)n[u(n)-u(n-20)]
 x(n)=n(0.9)n[u(n)-u(50)]
 x(n)={4,3,2,2,1,4,6,2}
 x(n)=(n+2)(-0.7)n-1u(n-2)
 x(n)=5(-0.9)ncos(0.1pin)u(n)
Bài tập

Đáp ứng tần số của hệ thống
 Viết chương trình Matlab để biểu diễn đáp ứng tần
số ở dạng phổ biên độ, phổ pha và dạng phần
thực, phần ảo của các hệ thống tuyến tính bất biến
được mô tả bởi phương trình sai phân sau:
( )jH e ω
• a).
• b).
( ) 0,5 ( 1) ( ) 2 ( 1) ( 2)y n y n x n x n x n− − = + − + −
5 1 1( ) ( 1) ( 2) ( 1)
6 6 3
y n y n y n x n− − + − = −
Ví dụ

Hệ thống cho bởi phương trình sai phân: 
y(n)=0.9y(n-1)+x(n)
 Xác định H(z) và biểu diễn các điểm không và 
điểm cực
 Vẽ |H(ejω)| và ∠ H(ejω)
 Xác định đáp ứng xung h(n)
Ví dụ (tiếp)

 Áp dụng công thức, ta có

 Matlab để kiểm tra:
b=[1];
a=[1,-0.9];
19.01
1)(
−
−
=
z
zH
zplane(b,a);

 Trong Matlab muốn tính H(ejω) ta sử dụng hàm freqz.
[H,w]=freqz(b,a,100);