%《MATLAB在電子信息課程中的應用》第六章例6.14程序q614
% 由差分方程遞推求輸出
% 電子工業出版社出版  陳懷琛 吳大正 高西全合著 2001年10月

clear
disp('輸入方程左端的系數向量a=[a(1),...a(na)]')
a = input('a = (書上取[1,0.1,0.15,-0.225])');
disp('輸入方程右端的系數向量b=[b(1),...b(nb)]=  ')
b = input('b=  (書上取[3,7,1])  ');
disp('輸入信號序列u(1:lu)（注意lu也就是計算長度,書上取exp(0.1*[1:20]）)')
u = input('u =  ');
na=length(a);nb=length(b);nu=length(u);
s=['起算點前',int2str(na-1),'點y的值 =[y(',int2str(na-2),'),..,y(0)]= 書上取[0,0,0]'];
ym=zeros(1,na+nu);ym(1:na-1) = input(s); 	% 建立ym序列并賦予初值
um = [zeros(1,na),u];					 % 建立um序列并賦予初值
for n=na:na+nu							% 這個n以ym的起點為準
   ys = ym(n-1:-1:n-na+1); us = um(n:-1:n-nb+1);	% 生成us及ys
   ym(n) = (b*us'-a(2:na)*ys')/a(1);		% 差分方程遞推求ym
end
% 把ym時間坐標右移na位,求出y
y = ym(na+1:na+nu);stem(y),grid
line([0,nu],[0,0])		% 給出起點和終點的x和y數組,畫橫坐標軸