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

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);