时域离散时间信号与系统

简单序列波形绘制

首先是利用matlab实现最基本的序列，包括单位采样序列、单位阶跃序列、矩形序列、三角波、方波、锯齿波、非周期方波、非周期三角波和Sinc函数，并绘制出响应波形。

过程并无太大不同，但总会有一个异曲同工的参照流程。首先，界定一个取值范围，表示你所要绘制的波形区域。
简单的举例
N=5;
n=0:N-1;
然后是范围内的赋值操作，除了直接赋值，亦可以利用zeros与ones函数
x=zeros(1,N);%表示由第一位到第N位均为0，
x(1)=1;%起始位为1，对应图里的0
值得注意的是，matlab中的序列起始位1所指代的坐标为0

当赋值完成之后，利用subplot、stem、plot、axis、title、grid等函数进行图像绘制。
值得注意的是stem与plot的差别，前者绘制序列，后者绘制连续波形
subplot(3,3,1);%框图布置，放置于3行3列第一张图
stem(n,x,’.’);%序列的点线画法，端点为‘.’，‘filled’为实心点
axis([-2 2 -1 2]);
%matlab中坐标轴的设定，依次为横坐标min，max，纵坐标min，max
grid;%grid作用为显示图片中的网线

提一手矩形序列的波形绘制采用的比较判断方式，以及*与.*的差别，如下
x=[(n-n1)>=0];%当n-n1大于等于0时判断赋给x=1，否则x=0
x1=[(n-n2+1)<=0];%同上
y=x.*x1;
%*与.在数乘中没有差别，但在矩阵运算时意义不同，ab表示矩阵a和b进行矩阵
%相乘，需要a行匹配b列，a.*b则表示a和b中元素按位置依次相乘，需要行列互相匹配。
%此处即为两者元素依次相乘，从而保留下两者均取值为1的点

三角波，等区间的划分，以及sawtooth（T,m）的应用，解释如下
fs=5000;
t=0:1/fs:0.1;%在0到0.1的范围内按1/fs的间隔划分区间
y=sawtooth(2pi50*t,0.5);
%生成三角波，前一变量规定周期，默认为2pi，后一变量表示在0-1T的范围内何时达到峰值，不填默认1，当其为0.5时，为三角波
%此处周期为（2pi）/（100pit）=0.02

同理的方波函数square（T,n），前者为周期，后者为占空比（0-100）
y=square(2pi60*t,50);%生成方波，前一变量同三角波，后一变量表示占空比，以百分号为单位，50表示占空比为50%

值得一提的是，sawtooth与square函数两者的周期默认为2pi，其实际周期则为（2pi）/T

%sinc函数 即Sa（t）函数=sint/t，基本照搬，与上述并无差别之处
t=-5:0.01:5;
y=sinc(t);%此处即为生成标准正弦函数
subplot(3,3,6);
plot(t,y);
axis([-5 5 -1 2]);
grid;

锯齿波的生成与三角波并无不同，只是sawtooth后一参数由0.5改至1而已

%非周期三角波，函数linspace，tripuls
t=linspace(-2,2,1000);%用于产生x1，x2之间的N点行矢量，其中x1为起点，x2为终点，N为元素个数，默认100
y1=tripuls(t,2);%生成单个三角波，变量依次为t序列长度，w三角波长度，缺省一个s变量，表示斜率，s默认为0
%s=0时，三角波顶点对应横轴0，不为0表示顶点平移了w/2*s个单位长度，底边不动。

%非周期方波，函数linspace以及rectpuls
t=linspace(-2,2,1000);
y2=rectpuls(t,1);
%同理tripuls，不过没有后面的s变量，依次表示T序列长度，w方波宽度

简单输入经简单系统的输出结果

设某LSI系统的单位冲激响应为h(n)=0.8^nu(n),当输入为三角形脉冲x(n)=tripuls([1:30]-15,20)时，求该系统的输出y(n)。试用MATLAB实现，并画出图形。

%三角波与幂函数的线性卷积的matlab运用
%利用tripuls生成单个三角波，利用conv计算卷积

x=tripuls([1:30]-15,20);%生成三角波，序列范围为1-30位，15表示生成的三角波的中心点坐标，20表示三角波底长，s默认为0
x1=[x,zeros(1,20)];%生成矩阵/数组？，前面为x长度为30，后面再加上20个连0,
N1=length(x);%计算x长度
n1=0:N1-1;%以0起始，所以后面长度要-1
N2=50;
n2=0:N2-1;
h=0.8.^n2;
y=conv(x,h);%表示y=x与h线性卷积
N=N1+N2-1;
n=0:N-1;%计算总长度
subplot(3,1,1);
stem(n2,x1,'filled');%端点为实心点，若没有默认端点为空心点
subplot(3,1,2);
stem(n2,h,'filled');
subplot(3,1,3);
stem(n(1:50),y(1:50),'filled');

%改进的卷积程序，若x(n),h(n)的起点不为0，可用conv_m计算卷积

function [y,ny]=conv_m(x,nx,h,nh)
nyb=nx(1)+nh(1);%计算卷积起点
nye=length(x)+length(h);
nye=nye+nyb-2;%计算卷积终点=两个函数的长度之和-2（去除初始点）
ny=[nyb:nye];
y=conv(x,h);

已知系统表达式，求h（n），s（n）及其他参量

例 2.9

%n∈[0,49]的系统单位冲激响应h（n）
%利用filter解方程，zeros，ones
B=0.866;A=[1,-0.8,0.64];%差分方程各项参数，A表示y（n）的系数，越往后阶次越低，B表示x（n）的系数，阶次同y
xn=[1,zeros(1,49)];%生成首项（0处）为1，后49项均为0的序列，即单位脉冲序列
%zeros（1,49）生成1行49列的零矩阵，ones同理
hn=filter(B,A,xn);%调用filter解差分方程，求系统输出信号h（n），由于输入x（n）为单位脉冲序列，故输出为单位冲激响应h（n）
n=0:length(hn)-1;%长度匹配
subplot(3,1,1);stem(n,hn,'.');grid;%绘图，subplot（311）为简略表示
title('(a)系统的单位冲激响应');%设置标题
xlabel('n');ylabel('h(n)');%xlabel与ylabel分别给x轴y轴作批注
%n∈[0,100]的系统的单位阶跃响应
xn=ones(1,100);%生成长度为100的单位阶跃序列，由0开始？
sn=filter(B,A,xn);%解差分方程，输出为单位阶跃响应
n=0:length(sn)-1;
subplot(3,1,2);stem(n,sn,'.');axis([0,30,0,2]);grid;
title('(b)系统的单位阶跃响应');
xlabel('n');ylabel('s(n)');
%利用一部分的h（n）前部后加0补位作为单位冲激响应生成单位阶跃响应，用FIR来代替或者说逼近IIR
for m=1:15
    hnfir(m)=hn(m);%去h（n）的前十五位
end
sn=filter(B,A,xn);%同理生成对应的单位阶跃序列
n=0:length(sn)-1;
subplot(3,1,3);stem(n,sn,'.');axis([0,30,0,2]);grid;
title('(c)FIR系统的单位阶跃响应');
xlabel('n');ylabel('h(n)');
%比较生成的（b）（c）波形，可以发现两者基本相同，说明某些IIR数字滤波器可以用FIR数字滤波器来逼近
%此时FIR滤波器的单位冲激响应可以通过截取IIR滤波器单位冲激响应的一段序列值来获得
%显然，当截取长度越长，逼近的误差越小。