Giáo trình Xử lý tín hiệu số I - Đỗ Huy Khôi

hư các tiện ích của nó
dành cho xử lý tín hiệu số. Để chương trình đơn giản và dễ dàng thấy được thuật
toán của nó, ta sẽ không thực hiện giao diện cho người dùng và chương trình được
viết theo cách đối thoại trực tiếp trên cửa sổ lệnh (Command Window) của
MATLAB, bằng cách dùng các lệnh disp và input. Hầu hết các chương trình sau
đây được viết dưới dạng Script và lưu vào các M-file cùng tên của chương trình.
Sau khi nhập vào một thư mục nào đó của MATLAB và tạo đường dẫn (nếu thư
mục này chưa có sẳn đường dẫn), để chạy chương trình, ta chỉ cần nhập tên
chương trình vào, trên Command Window, và gõ Enter. Nếu chương trình được
viết dưới dạng Function, người sử dụng cần nắm được các thông số vào, ra, để
nhập lệnh đúng cú pháp.
1. dsp13
% Nhập vào vector biến thời gian và biểu thức của tín hiệu, vẽ các loại tín hiệu: tương
tự, rời rạc, số.
%---------------------------------------------------
t=input('Nhap khoang thoi gian, VD:0:0.1:40, t= ');
y=input('Nhap ham so muon ve co bien t, VD:sin(t/4+1), y= ');
loai=input('(analog,type=1;discrete,type=2;digital,type=3)Type = ');
duong=input('(___,style=1;...,style=2;-.,style=3) stype = ');
if loai==1
 DS1=figure('Name','Type of signal','Color','w',...
 'NumberTitle','off','Position',[50 50 400 300]);
 if duong = =1
 plot(t,y,'r-');
198
 elseif duong = =2
 plot(t,y,'r:');
 elseif duong = =3
 plot(t,y,'r-.');
 end;
elseif loai= =2
 cham=input('(cham den,cham=1;cham trang,cham=2;...
 khong,cham=3) cham= ');
 DS1=figure('Name','Type of signal', 'Color','w',...
 'NumberTitle','off','Position',[50 50 400 300]);
 if cham= =1
 stem(t,y,'fulled');
 elseif cham= =2
 stem(t,y);
 elseif cham= =3
 stem1(t,y);
 end;
elseif loai= =3
 [x,z]=stairs(t,y);
 xt(1)=x(1);zt(1)=z(1);
 for n=1:length(x)/2-1
 ni=2*n+1;
 xt(n)=x(ni);zt(n)=z(ni);
 end;
 cham=input('(cham den,cham=1;cham trang,cham=2;...
 khong,cham=3) cham= ');
 plot(x,z,'g:');hold on;
 if cham= =1
 stem(xt,zt,'fulled');
 elseif cham= =2
 stem(xt,zt);
 elseif cham= =3
 stem1(xt,zt); 
 end;
end
axis off;
2. function dsphinh3_26(N,L)
%Ve bien do va pha cua DFT n diem cua day co do dai L.
% Doan hoa minh 2001
%--------------------------------------------------------
function dsphinh3_26(N,L)
xn=ones(1,L);
199
X=fft(xn,N);
X1=abs(X);
theta1=angle(X);
DS2=figure('Name','DFT N diem’,'Color','w',...
 'NumberTitle','off','Position',[50 50 580 300]);
stem(X1,'filled')
DS2=figure('Name','Type of signal','Color','w',...
 'NumberTitle','off','Position',[50 50 580 300]);
stem(theta1,'filled')
3. dsphinh5_16
% Ve dac tuyen cua mach loc thiet ke bang cua so co chieu dai bang 9 va bang 61.
% Doan Hoa Minh
syms w v;
y=sin((w-v)*9/2)/sin((w-v)/2);
z=int(y,v,-pi/4,pi/4);
z=simple(z)
w=0:0.01:pi;
for n=1:length(w)
Ht(n)=subs(z,'w',w(n));
end
H=exp(-j*4.*w)./(2*pi).*Ht;
tHt=abs(H);
Hdb=20*log10(tHt);
DS1=figure('Name','Type of signal','Color','w',...
 'NumberTitle','off','Position',[50 50 500 200]);
plot(w,tHt)
grid on
DS1=figure('Name','Type of signal','Color','w',...
 'NumberTitle','off','Position',[50 50 500 200]);
plot(w,Hdb,'k')
grid on
syms w v;
y1=sin((w-v)*61/2)/sin((w-v)/2);
z1=int(y1,v,-pi/4,pi/4);
z1=simple(z1)
w=0:0.01:pi;
for n=1:length(w)
Ht1(n)=subs(z1,'w',w(n));
end
H1=exp(-j*4.*w)./(2*pi).*Ht1;
tHt1=abs(H1);
200
Hdb1=20*log10(tHt1);
DS3=figure('Name','Type of signal',...
 'Color','w','NumberTitle','off','Position',[50 50 500 200]);
plot(w,tHt1)
grid on
DS4=figure('Name','Type of signal',...
 'Color','w','NumberTitle','off','Position',[50 50 500 200]);
plot(w,Hdb1,'k')
grid on
4.firequiripple
% Thiet ke bo loc FIR thong thap pha tuyen tinh dung thuat toan Remez exchange.
% Doan Hoa Minh
M=input('Nhap chieu dai cua dap ung xung, M = ');
dx=11; pdx=12;
disp('Chon dieu kien doi xung, neu doi xung thi nhap: dx')
disp(' , neu phan doi xung thi nhap: pdx')
dk=input('Dieu kien doi xung : ');
W=input('Nhap vector trong so,so phan tu bang so dai bang,...
 Vd: W=[1.2 1],W= ');
disp('Nhap vector cac tan so c../Anh bang tan,...
 mot cap tan so cho moi ')
disp('bang tan, cac tan so nay nam giua 0 va 1,...
 Vd F=[0 .1 .15 1]')
F=input('F = ');
disp('Nhap vector gia tri dap ung tan so mong muon A (gia tri thuc),')
disp('tai cac diem tan so bang c../Anh, A co kich thuoc bang F')
disp ('Vi du: A=[1 1 0 0]')
A=input('A = ');
N=M-1;
if dk= =11
 [hn,err]=remez(N,F,A,W)
elseif dk= =12
 [hn,err]=remez(N,F,A,W,'Hilbert')
end
w=0:0.001:pi;
f=w./pi;
H= freqz(hn,1,w); 
H1=20*log10(abs(H));
 DS1=figure('Name','Impulse Response','Color','w',...
 'NumberTitle','off','Position',[50 50 500 300]);
201
 n=0:1:M-1;
 stem(n,hn,'filled','k')
 axis off
 DS1=figure('Name','Frequency Response', 'Color','w',...
 'NumberTitle','off','Position',[50 50 500 300]);
 plot(f,abs(H),'k')
 grid on
 DS1=figure('Name','Frequency Response (dB)','Color','w',...
 'NumberTitle','off','Position',[50 50 500 300]);
 plot(f,H1,'k')
 ylim([-100 10])
 grid on
5. firsample
% Thiet ke bo loc FIR thong thap pha tuyen tinh bang phuong phap lay may tan so.
% Doan hoa minh 2001.
%---------------------------------------------------------
M=input('Nhap chieu dai cua dap ung xung, M = ');
dx=11; pdx=12;
disp('Chon dieu kien doi xung, neu doi xung thi nhap: dx')
disp(' , neu phan doi xung thi nhap: pdx')
dk=input('Dieu kien doi xung : ');
alpha=input('Chon he so alpha, alpha= ');
disp('Voi h(n) dx k=[0:(M-1)/2] neu M le,...
 k=[0:(M/2)-1] neu M chan')
disp('Voi h(n) pdx k=[0:(M-3)/2] neu M le,...
 k=[1:(M/2)] neu M chan')
disp('Nhap dac tuyen tan so mong muon,...
 tai cac diem tan so wk=2*pi*k/M')
if mod(M,2)= =0
 U=M/2-1;
else
 U=(M-1)/2;
end
for ii=1:U+1
 %kk=int2str(ii);
 %disp('k = 'kk);
 Hrk(ii)=input('Hr(k) = ');
end
G=zeros(U+1,1);
hn=zeros(M,1);
for k=1:U+1
 G(k)=((-1)^(k-1))*Hrk(k);
end
202
if alpha= =0
 if dk= =11
 for n=1:M
 for k=2:U+1
 hn(n)=hn(n)+G(k)*cos(pi*(k-1)*(2*(n-1)+1)/M);
 end
 hn(n)=(2*hn(n)+G(1))/M;
 end
 elseif dk = =12
 if mod(M,2)= =1
 for n=1:M
 for k=1:U+1
 hn(n)=hn(n)-2*G(k)*sin(2*pi*(k-1)*((n-1)+0.5)/M)/M;
 end
 end
 else
 for n=1:M
 for k=1:U
 hn(n)=hn(n)-2*G(k)*sin(2*pi*k*((n-1)+0.5)/M)/M;
 end
 hn(n)=hn(n)+((-1)^(n))*G(U+1)/M;
 end
 end
 end
elseif alpha= = 0.5
 if dk= =11
 for n=1:M
 for k=1:U+1
 hn(n)=hn(n)+2*G(k)*sin(2*pi*(k-1+1/2)*((n-1)+0.5)/M)/M;
 end
 end
 elseif dk= =12
 for n=1:M
 for k=1:U+1 
 hn(n)=hn(n)+2*G(k)*cos(2*pi*(k-1+1/2)*((n-1)+0.5)/M)/M;
 end
 end
 end
end
hn
om=0:0.01:pi;
if mod(M,2)= =0
 Hr=hn(1).*cos(om.*((M-1)/2));
 n=1;
203
 while n<=U
 n=n+1;
 Hr=Hr+hn(n).*cos(om.*((M-1)/2-n+1));
 end
 Hr=2.*Hr;
else
 Hr=hn(1).*cos(om.*((M-1)/2));
 n=1;
 while n<=(M-3)/2
 n=n+1;
 Hr=Hr+hn(n).*cos(om.*((M-1)/2-n+1));
 end
 Hr=2.*Hr;
 Hr=Hr+hn(U+1);
end
modunH=abs(Hr);
DS1=figure('Name', 'Dap ung bien do', 'Color','w',...
 'NumberTitle','off','Position',[50 50 400 300]);
plot(om,modunH,'k');
grid on
modunHdb=20.*log10(modunH);
DS2=figure('Name','Type of signal', 'Color','w',...
 'NumberTitle','off','Position',[50 50 400 300]);
plot(om,modunHdb,'k');
grid on 
teta=-om.*(M-1)/2+angle(Hr);
DS3=figure('Name','Dap ung pha','Color','w',...
 'NumberTitle','off','Position',[50 50 400 300]);
plot(om,teta,'k');
grid on
DS4=figure('Name','Dap ung xung', 'Color','w',...
 'NumberTitle','off','Position',[50 50 400 300]);
stem(hn,'filled','k');
grid on
6. dsphinh 5_15
% Ve dap ung tan so cua cua so chu nhat co chieu dai bang M=9,
% M=51 va m=101.
% Doan hoa minh 2001
om=0:0.001:pi;
204
M=9;
W1=20*log10(abs(sin(om.*M/2)./sin(om./2)));
DS1=figure('Name','Dap ung tan so cua cua so chu nhat M=9',
 'Color','w','NumberTitle','off','Position',[50 50 500 200]);
plot(om,W1,'k')
title('M = 9');
xlabel('w (rad)');
ylabel('|W(w)|(dB)');
axis on
grid on
M=51;
W2=20*log10(abs(sin(om.*M/2)./sin(om./2)));
DS2=figure('Name',' Dap ung tan so cua cua so chu nhat M=51',
 'Color','w','NumberTitle','off','Position',[50 50 500 200]);
plot(om,W2,'k')
title('M = 51');
xlabel('w (rad)');
ylabel('|W(w)|(dB)');
axis on
grid on
M=101;
W3=20*log10(abs(sin(om.*M/2)./sin(om./2)));
DS2=figure('Name','Dap ung tan so cua cua so chu nhat M=101',
 'Color','w','NumberTitle','off','Position',[50 50 500 200]);
plot(om,W3,'k')
title('M = 101');
xlabel('w (rad)');
ylabel('|W(w)|(dB)');
axis on
grid on 
function hh = stem1(varargin)
% Hàm này được cải biên từ hàm stem của MATLAB, vẽ dãy rời rạc không có chấm
trêm đầu.
%STEM1 Discrete sequence or "stem" plot.
% STEM1(Y) plots the data sequence Y as stems from the x axis
%
% STEM1(X,Y) plots the data sequence Y at the values specfied
% in X.
% STEM1(...,'LINESPEC') uses the linetype specifed for the stems and
% markers. See PLOT for possibilities.
%
% H = STEM(...) returns a vector of line handles.
%
205
% See also PLOT, BAR, STAIRS.
% Copyright (c) by Doan Hoa Minh.
% Date: 2000/6/4.
nin = nargin;
fill = 0;
ls = '-';
ms = 'o';
col = '';
% Parse the string inputs
while isstr(varargin{nin}),
 v = varargin{nin};
 if ~isempty(v) & strcmp(lower(v(1)),'f')
 fill = 1;
 nin = nin-1;
 else
 [l,c,m,msg] = colstyle(v);
 if ~isempty(msg), 
 error(sprintf('Unknown option "%s".',v));
 end
 if ~isempty(l), ls = l; end
 if ~isempty(c), col = c; end
 if ~isempty(m), ms = m; end
 nin = nin-1;
 end
end
error(nargchk(1,2,nin));
[msg,x,y] = xychk(varargin{1:nin},'plot');
if ~isempty(msg), error(msg); end
if min(size(x))= =1, x = x(:); end
if min(size(y))= =1, y = y(:); end
% Set up data using fancing ../indexing
[m,n] = size(x);
xx = zeros(3*m,n);
xx(1:3:3*m,:) = x;
xx(2:3:3*m,:) = x;
xx(3:3:3*m,:) = NaN;
[m,n] = size(y);
yy = zeros(3*m,n);
yy(2:3:3*m,:) = y;
yy(3:3:3*m,:) = NaN;
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;
206
h2 = plot(xx,yy,[col,ls],'parent',cax);
if nargout>0, hh = h2; end