实验一  时域离散信号与系统变换域分析
一、实验目的
1．了解时域离散信号的产生及基本运算实现。
2．掌握离散时间傅里叶变换实现及系统分析方法。
3. 熟悉离散时间傅里叶变换性质。
4. 掌握系统Z域分析方法。
5. 培养学生运用软件分析、处理数字信号的能力。
二、实验内容 
1. 序列的基本运算
1.1 产生余弦信号及带噪信号  0<=n<=50（噪声采用randn函数）
1.2 已知 ， ,求两个序列的和、乘积、序列x1的移位序列（右移2位）,序列x2的翻褶序列，画出原序列及运算结果图。
2. 序列的傅里叶变换
2.1 已知序列。试求它的傅里叶变换，并且画出其幅度、相角、实部和虚部的波形，并分析其含有的频率分量主要位于高频区还是低频区。
2.2 令，求其傅立叶变换。分别用和对其进行采样，求出离散时间傅立叶变换，写出程序，并画出相应频谱，分析结果的不同及原因。
3. 序列的傅里叶变换性质分析
3.1 已知序列，，求其傅里叶变换，并讨论其傅里叶变换的周期性和对称性。
3.2 已知序列，，求其傅里叶变换，并讨论其傅里叶变换的周期性和对称性。
为了方便，考虑在两个周期，例如[]中2M+1个均匀频率点上计算FT，并且观察其周期性和对称性。为此给出function文件如下，求解FT变换：
function[X,w]=ft1(x,n,k)
w=(pi/abs(max(k)/2))*k
X=x*(exp(-j*pi/abs(max(k)/2))).^(n'*k)
3.3 编写程序验证序列傅里叶变换频移性质，时域卷积定理（时域卷积后的频域特性）。（所需信号自行选择）
4. 时域差分方程的求解
4.1求解差分方程y(n)＋a1y(n-1)＋a2y(n-2)=b0x(n)＋b1x(n-1)的零状态响应和全响应。已知X(n)为单位取样序列，y(-1)=1,y(-2)=2,a1=0.5,a2=0.06,b0=2,b1＝3。
5. 离散系统的Z域分析
5.1 利用系统函数分析系统的稳定性。假设系统函数如下式：
，试判断系统是否稳定。
5.2 已知线性时不变系统的系统函数，编写程序求其单位取样响应，频率响应及系统零极点，并画出相应图形。
6. 创新训练拓展内容
6.1 利用Matlab自带的录音功能，或利用Goldwave等音频编辑软件，对语音或其他音频信号进行采集并保存为*.wav文件。
要求：（1）采用不同的采样频率（2000Hz，4000Hz，8000Hz，16000Hz等）。
（2）对采集得到的信号进行播放。
（3）分析在不同采样频率下得到的信号有何不同。
6.2 设定一个连续时间信号，进行抽样和恢复，要求分析不同采样频率对恢复结果的影响，给出实验程序及各关键步骤图形结果。
6.3 设计内容：
设计一个离散系统，给定系统函数或差分方程，设定激励及初始条件。要求：
（1）绘制系统函数零极点图，判断稳定性；
（2）求单位序列响应h（n）；
（3）求系统零输入响应及零状态响应，要求零状态响应采样三种方法求解（卷积的方法、迭代解法、差分方程求解函数方法），激励自定；
（4）分析系统频响特性，画出频响函数幅频曲线和相频曲线。
三、试验要求
第一部分：验证试验内容
根据给定的试验内容，部分试验给出了参考程序段，见下面各段程序。请基于Matlab环境进行验证试验。
第二部分：编程试验内容
对于给定的试验内容中，没有参考程序段的部分，进行编程，并给出试验结果，进行相应的分析。
第三部分：创新训练拓展内容
此部分内容，要求根据个人能力，进行选作。
1. 序列的基本运算
%1.单位取样序列  x(n)=delta(n-n0)    要求n1<=n0<=n2 
function[x,n]=impseq(n0,n1,n2)
n=[n1:n2];  x=[(n-n0)==0];    == 是逻辑判断
%2.单位阶跃序列  x(n)=u(n-n0)      要求n1<=n0<=n2
function[x,n]=stepseq(n0,n1,n2)
n=[n1:n2];    x=[(n-n0)>=0];
%3.信号加        y(n)=x1(n)+x2(n)
%find函数：出非零元素的索引号
%x1:第一个序列的值，n1:序列x1的索引号
%x2:第二个序列的值，n2:序列x2的索引号
function[y,n]=sigadd(x1,n1,x2,n2)
n=min(min(n1),min(n2)):max(max(n1),max(n2));
y1=zeros(1,length(n));  y2=y1;
y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;     
y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;
y=y1+y2;
%4.信号乘      y(n)=x1(n)*x2(n)
function[y,n]=sigmult(x1,n1,x2,n2)
n=min(min(n1),min(n2)):max(max(n1),max(n2));
y1=zeros(1,length(n));  y2=y1;
y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;
y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;
y=y1.*y2;
%5.移位    y(n)=x(n-n0)
function[y,n]=sigshift(x,m,n0)
n=m+n0;  y=x;
%6.翻褶    y(n)=x(-n)
function[y,n]=sigfold(x,n)
y=fliplr(x);  n=-fliplr(n);
2. 序列的傅里叶变换
%7. 求序列的傅里叶变换
w=[0:1:500]*pi/500
X=exp(j*w)./(exp(j*w)-0.5*ones(1,501))
magX=abs(X)
angX=angle(X)
realX=real(X)
imagX=imag(X)
subplot(2,2,1)
plot(w/pi,magX)
grid
xlabel('frequency in pi units')
title('Magnitude Part')
ylabel('Magnitude')
subplot(2,2,3)
plot(w/pi,angX)
grid
xlabel('frequency in pi units')
title('Angle Part')
ylabel('Radians')
subplot(2,2,2)
plot(w/pi,realX)
grid
xlabel('frequency in pi units')
title('Real Part')
ylabel('Real')
subplot(2,2,4)
plot(w/pi,imagX)
grid
xlabel('frequency in pi units')
title('Imaginary Part')
ylabel('Imaginary')
  程序执行结果：
matlab求傅里叶变换%8 令，绘制其傅立叶变换。用不同频率对其进行采样，分别画出。
Dt=0.00005;    %步长为0.00005s
t=-0.005:Dt:0.005;   
xa=exp(-1000*abs(t));  %取时间从-0.005s到0.005s这段模拟信号
Wmax=2*pi*2000;  %信号最高频率为2*2000
K=500;        %频域正半轴取500个点进行计算
k=0:1:K;
W=k*Wmax/K;    %  求模拟角频率
Xa=xa*exp(-j*t'*W)*Dt;      %计算连续时间傅立叶变换（利用矩阵运算实现）
Xa=real(Xa);  %取实部
W=[-fliplr(W),W(2:501)];    %将角频率范围扩展为从-到+
Xa=[fliplr(Xa),Xa(2:501)]; 
subplot(2,2,1);