史上最全数字信号处理实验报告完美版

实验一、零极点分布对系统频率响应的影响
Y(n)=x(n)+ay(n-1)
1、调用MATLAB函数freqz计算并绘制的幅频特性和相频特性其中：1 代表a=0.7；2代表a=0.8；3代表a=0.9
a=0.7时的零极点图
A=0.8时的零极点图
a=0.9时的零极点图
观察零极点的分布与相应曲线易知：
小结：系统极点z=a，零点z=0，当B点从w=0逆时针旋转时，在w=0点，由于极点向量长度最短，形成波峰,并且当a越大，极点越接近单位圆，峰值愈高愈尖锐；在w=pi点形成波谷；z=0处零点不影响幅频响应
2、先求出系统传函的封闭表达式，通过直接计算法得出的幅频特性和相频特性曲线。

其中：1代表a=0.7；2代表a=0.8；3代表a=0.9
附录程序如下：（对程序进行部分注释）
>> a=0.7;
w=0:0.01:2*pi;%设定w的范围由0到2π，间隔为0.01
y=1./(1-a*exp(-j*w)); %生成函数
subplot(211);plot(w/2/pi,10*log(abs(y)),'g');%生成图像其中通过调用abs函数计算幅值
hold on;
xlabel('Frequency(Hz)');%定义横坐标名称
ylabel('magnitude(dB)');%定义纵坐标名称
title('a=0.8,直接计算h(ejw)');grid on;%定义图片标题
subplot(212);plot(w/2/pi,unwrap(angle(y)),'g');grid on;%生成图像其中通过调用angle计算相角，‘g’为规定线条颜色
hold on;
>> a=0.8;
w=0:0.01:2*pi;
y=1./(1-a*exp(-j*w));
subplot(211);plot(w/2/pi,10*log(abs(y)),'r');
hold on;
xlabel('Frequency(Hz)');
ylabel('magnitude(dB)');
title('a=0.8,直接计算h(ejw)');grid on;
subplot(212);plot(w/2/pi,unwrap(angle(y)),'r');grid on;
hold on;
>> a=0.9;
w=0:0.01:2*pi;
y=1./(1-a*exp(-j*w));
subplot(211);plot(w/2/pi,10*log(abs(y)),'b');
hold on;
xlabel('Frequency(Hz)');
ylabel('magnitude(dB)');
title('a=0.9,直接计算h(ejw)');grid on;
subplot(212);plot(w/2/pi,unwrap(angle(y)),'b');grid on;
hold on;
2、y(n)=x(n)=ax(n-1)
通过调用freqz函数绘图，其中：1代表a=0.7,；2代表a=0.8；3代表a=0.9
附录程序如下：（因为程序同实验一相同不再进行注释）a=0.7;
A=1;
B=[1,a];
freqz(B,A,256,'whole',1);
title('a=0.7');
hold on;
a=0.8;
A=1;
B=[1,a];
freqz(B,A,256,'whole',1);
title('a=0.8');
hold on;
a=0.9;
A=1;
B=[1,a];
freqz(B,A,256,'whole',1);
title('a=0.9');
以下为a为不同数值时的零极点图
a=0.7
A=0.8
A=0.9
小结：系统极点z=0，零点z=a，当B点从w=0逆时针旋转时，在w=0点，由于零点向量长度最长，形成波峰：在w=pi点形成波谷；z=a处极点不影响相频响应。

3、y(n)=1.273y(n-1)-0.81y(n-2)+x(n)+x(n-1)
通过调用freqz函数绘制响应曲线如下：
附录程序如下：
A=[1,-1.273,0.81];
B=[1,1];
freqz(B,A,256,'whole',1);
峰值频率：
谷值频率：
零极点分布图：
小结：系统极点z1=0.79+j0.62*1.62^(-2)，z2=0.79-j0.62*1.62^(-2)零点z1=-1，z2=0当B点从w=0逆时针旋转时，当旋转到接近极点z1=0.79+j0.62*1.62^(-2)是极点向量长度最短，幅度特性出现峰值。

当转到w=pi点形成波谷；z=a处零点不影响幅频响应。

当旋转到接近极点z2=0.79-j0.62*1.62^(-2)是极点向量长度再次最短，幅度特性再次出现峰值。

实验总结：
当旋转点转到极点附近时，幅度特性出现峰值，并且极点越靠近单位圆，峰值越尖锐。

如果极点在单位圆上，系统不稳定。

当旋转点转到零点附近时，幅度特性出现谷值，并且零点越靠近单位圆，谷值越接近零，零点在单位圆上时，谷值为零。

所以：极点位置主要影响频响的峰值位置及尖锐程度，零点主要影响频响的谷值位置及形状。

实验二、信号截断及补零对频谱的影响
N<30时，取N=20
1）N=20 补零2*30
T=0.001;
t=-0.2:T:0.2;%设定时间间隔
f1=100;f2=120;%设定频率
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);%编写函数
figure(1);subplot(211);plot(t,x);xlabel('t/s');title('连续信号及频谱图');grid on;%绘制第一幅图以及命名。

omg1=f1;omg2=f2;%设定创建函数的变量
X_freq=1/2*(impseq(omg1,-150,150)+impseq(-omg1,-150,150)+impseq(omg2,-150,150)+impseq( -omg2,-150,150));%通过调用创建的函数impseq来求取相应频率对应的幅度值
f_hz=-150:150;
figure(1);subplot(212);plot(f_hz,X_freq);xlabel('f/Hz');ylabel('幅值');grid on;%绘图以及标明X、y 轴的单位
下图为n=20时的没有补零以及补零2*30、7*30、20*30后的频谱图。

附录程序如下：
fs=600;x2=0;%设定采样信号的频率
x2=cos(2*pi*f1*(0:40)/fs)+cos(2*pi*f2*(0:40)/fs);%创建函数，并截断函数其范围为0到40
N=20;L=20;%N表示采样点数,L表示进行傅立叶变换的点数
x3=0;f_axis=0;%初始化工作
x3=x2(1:N);
x_pu=fft(x3);一维傅里叶变换，将时域信号转换为频域信号。

f_axis=k_to_fs(fs,L);%%%%%%%%%%%使横坐标k能用频率fs表示%%%%
%也可以用如下更为简单的函数形式实现
%f_axis=-fs/2+(0:L-1)*fs/l;
figure(2);subplot(411);plot(f_axis,abs(fftshift(x_pu)));xlabel('f/Hz,N=20,没有补零');ylabel('幅值');grid on; hold on;
stem(f_axis,abs(fftshift(x_pu)));hold off;
%采样信号补零后的频谱图
L=20+2*30;
x_pu=fft(x3,L);
f_axis=k_to_fs(fs,L);
figure(2);subplot(412);plot(f_axis,abs(fftshift(x_pu)));xlabel('f/Hz,N=20,补零2*30');ylabel('幅值');grid on; hold on;
stem(f_axis,abs(fftshift(x_pu)));hold off;绘制离散序列的火柴杆图。

其中通过fftshift函数实现对fft傅里叶变换后的坐标处理，使得零频率分量移到坐标中心。

L=20+7*30;
x_pu=fft(x3,L);一维傅里叶变换，并补l个个数个零。

f_axis=k_to_fs(fs,L);
figure(2);subplot(413);plot(f_axis,abs(fftshift(x_pu)));xlabel('f/Hz,N=20,补零7*30');ylabel('幅值');grid on; hold on;
stem(f_axis,abs(fftshift(x_pu)));hold off;
L=20+20*30;
x_pu=fft(x3,L);
f_axis=k_to_fs(fs,L);
figure(2);subplot(414);plot(f_axis,abs(fftshift(x_pu)));xlabel('f/Hz,N=20,补零20*30');ylabel('幅值');grid on; hold on;
stem(f_axis,abs(fftshift(x_pu)));hold off;
以下为实验中不同补零情况的的分图（为了更清楚的观察）
2）N=20 补零7*30
N=20 补零20*30
N>30时，取N=40
1）N=40 补零2*30
T=0.001;
t=-0.2:T:0.2;
f1=100;f2=120;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
figure(1);subplot(211);plot(t,x);xlabel('t/s');title('连续信号及频谱图');grid on;
omg1=f1;omg2=f2;
X_freq=1/2*(impseq(omg1,-150,150)+impseq(-omg1,-150,150)+impseq(omg2,-150,150)+impseq( -omg2,-150,150));
f_hz=-150:150;
figure(1);subplot(212);plot(f_hz,X_freq);xlabel('f/Hz');ylabel('幅值');grid on;
%采样信号的频谱
下图为n=40时，没有补零、补零2*30、7*30、20*30时的频谱图。

附录程序（除了改变参数外，其他均相同，故不再对函数进行注释）
fs=600;x2=0;
x2=cos(2*pi*f1*(0:40)/fs)+cos(2*pi*f2*(0:40)/fs);
N=40;L=40;%N表示采样点数,L表示进行傅立叶变换的点数
x3=0;f_axis=0;%初始化工作
x3=x2(1:N);
x_pu=fft(x3);
f_axis=k_to_fs(fs,L);%%%%%%%%%%%使横坐标k能用频率fs表示%%%%
%也可以用如下更为简单的函数形式实现
%f_axis=-fs/2+(0:L-1)*fs/l;
figure(2);subplot(411);plot(f_axis,abs(fftshift(x_pu)));xlabel('f/Hz,N=40,没有补零');ylabel('幅值');grid on; hold on;
stem(f_axis,abs(fftshift(x_pu)));hold off;
%采样信号补零后的频谱图
L=40+2*30;
x_pu=fft(x3,L);
f_axis=k_to_fs(fs,L);
figure(2);subplot(412);plot(f_axis,abs(fftshift(x_pu)));xlabel('f/Hz,N=40,补零2*30');ylabel('幅值');grid on; hold on;
stem(f_axis,abs(fftshift(x_pu)));hold off;
L=40+7*30;
x_pu=fft(x3,L);
f_axis=k_to_fs(fs,L);
figure(2);subplot(413);plot(f_axis,abs(fftshift(x_pu)));xlabel('f/Hz,N=40,补零7*30');ylabel('幅值');grid on; hold on;
stem(f_axis,abs(fftshift(x_pu)));hold off;
N=40 补零20*30
更改后程序如下：
%采样信号补零后的频谱图
L=40+20*30;
x_pu=fft(x3,L);
f_axis=k_to_fs(fs,L);
figure(2);subplot(414);plot(f_axis,abs(fftshift(x_pu)));xlabel('f/Hz,N=40,补零20*30');ylabel('幅值');grid on; hold on;
stem(f_axis,abs(fftshift(x_pu)));hold off;
以下为分开的图片：
N=40 补零7*30
附录创建的另个函数的程序：
impseq函数
function[f,omg]=impseq(omg0,omg1,omg2)
omg=[omg1:omg2];
f=[(omg-omg0)==0];
k_to_fs函数
function f_axis=k_to_fs(fs,N)
for k=0:N-1
if(k<=N/2-1)
f_axis(k+1+N/2)=fs/N*k;
else
f_axis(k+1-N/2)=fs/N*(k-N);
end
end
实验总结：
信号截断对频谱的影响：
运用计算机实现工程测试信号处理时，不可能对无限长的信号进行测量和运算，而是取其有限的时间片段进行分析。

做法是从信号中截取一个时间片段，从而对频谱产生两方面的影响：一是会影响频谱分辨率；二是使频谱产生泄漏。

截断时采用的窗函数的主瓣宽度N影响频率分辨率（信号的频率分辨率为Fs/N），而边瓣则会影响频谱的泄漏，产生能量的泄漏，使得信号失真。

补零对于频谱的影响：
由以上几幅图进行对比可知：对于n=20和n=40时分辨率有着明显的不同，n=40时可以明显的区分两个频率，而在n=20时则不能。

所以当提高窗函数的宽度的时候可以明显的提高物理分辨率。

而对于补零来说并没有增加有效的数据长度T或M（宽度），因而不可能增加原数据的新的信息，并不能依靠补零将两个靠的很近的谱峰分开，只能提高计算分辨率，不能提高物理分辨率。

但补零一方面可以使数据为2的N次幂，以便于快速傅里叶换算，并且在一定程度上补零可以对截断时产生的泄漏进行克服。

问答题：
1、答：对于连续时间周期信号而言，其Fourier级数就是他的一个周期的截取后的非周期信号的的傅立叶变换采样，连续时间信号采样后所得到的离散信号的DTFT可看成原来连续时间傅立叶变换在横轴做一下模拟——数字频率变换后进行周期延拓而成。

离散傅里叶变换可以看成DTFT在主值区间（0到2*pi）的等间隔采样。

2、答：幅值相同，只是补零后的频谱相当于在补零前的频谱中插入（2^n)-1条谱线。

与补零前的频谱中相重合的谱线，它们的幅值和相位完全一致。

3、答：混叠现象：由于采样信号频谱发生变化，而出现高、低频成分发生混淆的一种现象。

抽样时频率不够高，抽样出来的点既代表了信号中的低频信号的样本值，也同时代表高频信号样本值，在信号重建的时候，高频信号被低频信号代替，两种波形完全重叠在一起，形成严重失真。

泄漏现象：在实际问题中遇到的离散时间序列x(n)通常是无限长序列，因而处理这个序列的时候需要将它截短。

截短相当于将序列乘以窗函数w(n)。

根据频域卷积定理，时域中x(n)和w(n)相乘对应于频域中它们的离散傅立叶变换X(jw)和W(jw)的卷积。

因此，x(n)截矩后的频谱不同于它以前的频谱
栅栏现象：对一函数实行采样，即是抽取采样点上的对应的函数值。

其效果如同透过栅栏的缝隙观看外景一样，只有落在缝隙前的少数景象被看到，其余景象均被栅栏挡住而视为零，这种现象称为栅栏效应，不管是时域采样还是频域采样，都有相应的栅栏效应。

只是当时域采样满足采样定理时，栅栏效应不会有什么影响。

而频域采样的栅栏效应则影响很大。

实验体会：
Matlab在工程中应用非常广泛，但是由于学得不够好，给本次实验带来不小问题，要努力学习使用MATLAB了。