北京交通大学《信号与系统》研究性学习实验报告

1、信号与系统课程研究性学习手册姓名 学号 同组成员 14* 14* 14* 指导教师 时间 1. 信号的时域分析专题研讨【目的】(1) 掌握基本信号及其特性，了解实际信号的建模。(2) 掌握基本信号的运算，加深对信号时域分析基本原理和方法的理解，并建立时频之间的感性认识。(3) 学会仿真软件MATLAB的初步使用方法，掌握利用MATLAB进行信号表示和信号运算。【研讨内容】题目1：基本信号的产生，语音的读取与播放1） 生成一个正弦信号，改变正弦信号的角频率和初始相位，观察波形变化，并听其声音的变化。2） 生成一个幅度为1、基频为2Hz、占空比为50%的周期方波。3） 观察一定时期内的股票上证指数

2、变化，生成模拟其变化的指数信号。4） 分别录制一段男声、女声信号，进行音频信号的读取与播放，画出其时域波形。【温馨提示】(1)利用MATLAB函数 wavread(file)读取.wav格式文件。(2)利用MATLAB函数 sound(x, fs)播放正弦信号和声音信号。【题目分析】【仿真程序】1）生成一个正弦信号t=0:0.001:8;y=sin(2*pi*t+pi/6);plot(t,y) 改变其角频率和初始相位t=0:0.001:8;y=sin(pi*t+pi/2);plot(t,y)2）生成一个幅度为1、基频为2Hz、占空比为50%的周期方波t=0:0.001:10;y=square(

3、2*t,50);plot(t,y);axis(0,10,-1.2,1.2)3）观察一定时期内的股票上证指数变化，生成模拟其变化的指数信号。x1=0:0.0015;y1=2630+1.75*exp(x1);x2=5:0.001:10;y2=2895-1.54*exp(0.8*x2);x3=10:0.001:15;y3=2811+152*exp(-0.08*x3);x4=15:0.001:20;y4=2600-151*exp(-0.08*x4);x=x1,x2,x3,x4;y=y1,y2,y3,y4;plot(x,y);4）分别录制一段男声、女声信号，进行音频信号的读取与播放，画出其时域波形。y,

4、x,nbits=wavread('C:UsersiDesktop信号研讨图片文件日野聪.wav');sound(y,x); plot(y)y,x,nbits=wavread('C:UsersiDesktop信号研讨图片文件钉宫理惠.wav');sound(y,x); plot(y)【仿真结果】1）生成一个正弦信号 改变其角频率和初始相位2）生成一个幅度为1、基频为2Hz、占空比为50%的周期方波3）观察一定时期内的股票上证指数变化，生成模拟其变化的指数信号。4）分别录制一段男声、女声信号，进行音频信号的读取与播放，画出其时域波形。【结果分析】(1) 随着正弦信号

5、角频率的变化，其波形变得更紧密，声音逐渐变得尖细而高。(2) 男声多低沉粗犷，主要是因为男声中低频分量更多；女声多高亢清脆，主要是因为女声中高频分量更多。【自主学习内容】wavread的读取与sound函数的使用，音频文件的放置。【阅读文献】MATLAB语言与实践教程（第2版）、信号与系统【发现问题】 (专题研讨或相关知识点学习中发现的问题)：根据声音信号的什么特征能有效区分出男声和女声？【问题探究】 高频低频成分的多少，以及声音的尖细程度【研讨内容】题目2：信号的基本运算（语音信号的翻转、展缩）1） 将原始音频信号在时域上进行延展、压缩，2） 将原始音频信号在时域上进行幅度放大与缩小，3）

6、将原始音频信号在时域上进行翻转，【题目分析】【仿真程序】 原始音频信号y,x,nbits=wavread('D:音乐动画歌曲（207曲）灼眼的夏娜无损 17曲Light My Fire.wav');plot(y)1） 将原始音频信号在时域上进行延展、压缩 将原始音频信号在时域上延展3倍y,x,nbits=wavread('D:音乐动画歌曲（207曲）灼眼的夏娜无损 17曲Light My Fire.wav');sigLength=length(y);t=3*(0:sigLength-1)/x;figure;plot(t,y);xlabel('Time(s

7、)'); 将原始音频信号在时域上压缩为1/5y,x,nbits=wavread('D:音乐动画歌曲（207曲）灼眼的夏娜无损 17曲Light My Fire.wav');sigLength=length(y);t=1/5*(0:sigLength-1)/x;figure;plot(t,y);xlabel('Time(s)');2）将原始音频信号在时域上进行幅度放大与缩小 将原始音频信号的幅度扩大为2倍y,x,nbits=wavread('D:音乐动画歌曲（207曲）灼眼的夏娜无损 17曲Light My Fire.wav');sigLe

8、ngth=length(y);m=fft(y,sigLength);Pyy=m.*conj(m)/sigLength;halflength=floor(sigLength/2);f=x*(0:halflength)/sigLength;t=(0:sigLength-1)/x;figure;plot(t,2*y);xlabel('Time(s)'); 将原始音频信号的幅度缩小为1/2y,x,nbits=wavread('D:音乐动画歌曲（207曲）灼眼的夏娜无损 17曲Light My Fire.wav');sigLength=length(y);m=fft(y,

9、sigLength);Pyy=m.*conj(m)/sigLength;halflength=floor(sigLength/2);f=x*(0:halflength)/sigLength;t=(0:sigLength-1)/x;figure;plot(t,1/2*y);xlabel('Time(s)');3）将原始音频信号在时域上进行翻转y,x,nbits=wavread('D:音乐动画歌曲（207曲）灼眼的夏娜无损 17曲Light My Fire.wav');sigLength=length(y);t=-1*(0:sigLength-1)/x;figure

10、;plot(t,y);xlabel('Time(s)');【仿真结果】 原始音频信号1）将原始音频信号在时域上进行延展、压缩 将原始音频信号在时域上延展3倍 将原始音频信号在时域上压缩为1/52）将原始音频信号在时域上进行幅度放大与缩小 将原始音频信号的幅度扩大为2倍 将原始音频信号的幅度缩小为1/23）将原始音频信号在时域上进行翻转【结果分析】以时间为横轴的图比以频率为横轴的图看起来复杂得多，时域上的分析在现阶段很困难。对于原始音频信号，时域上的分析太复杂，而如果在频域上分析，则简化很多。【自主学习内容】信号的获取以及如何使用；MATLAB的一些基本语法；如何根据MATLAB

11、语法来编写自己想表达的语句。【阅读文献】MATLAB语言与实践教程（第2版）、信号与系统【发现问题】 (专题研讨或相关知识点学习中发现的问题)：在查找资料时，发现音频单声道和双声道的编写此程序会有不同，同时语音信号的录入需要我们的注意。【问题探究】 2.系统的时域分析专题研讨【目的】(1) 掌握系统响应的时域求解，加深对系统时域分析基本原理和方法的理解。(2) 掌握连续系统零状态响应（卷积积分）数值计算的方法。(3) 学会仿真软件MATLAB的初步使用方法，掌握利用MATLAB求解连续系统和离散系统的零状态响应。(4) 培养学生发现问题、分析问题和解决问题的能力。【研讨内容】题目1：系统响应时

12、域求解1） 求一个RLC电路的零输入响应和零状态响应，2） 将原始音频信号中混入噪声，然后用M点滑动平均系统对受噪声干扰的信号去噪，改变M点数，比较不同点数下的去噪效果，【题目分析】（1） 一个RLC电路，若,L=1H,C=1F,电容上的初始储能为Vc=2V，电感初始储能为i=1A，试求输入激励为X（t）时的零输入响应和零状态响应。（2） 题目要求采用M点滑动平均系统进行去噪。M点滑动平均系统可以看成是N=0的差分方程。调用filter函数时，调用参数a-1=1,b为有M个元素的向量，b中每个元素的值为1/M。即M点的滑动平均系统输入输出关系为：，同时我们将噪声设为n,函数为n=rand(n,

13、1);原始信号为s。通过调整M值，观察和比较去噪效果，从而得出结论。【仿真程序】（1）零输入响应 num=6;den=1 3 6;R0=2 1; sys=tf(num,den); sys1=ss(sys); t=0:0.01:8; u=0*t; lsim(sys1,u,t,R0); axis(0 8 -0.5 2.5）零状态响应num=6;den=1 3 6; sys=tf(num,den); t=0:0.01:8; xt=10*sin(2*pi*t); yt=lsim(sys,xt,t); plot(t,yt,r); axis(0 8 -2 4)（2）fs=44100;bits=16;R=1

14、00000y,fs,bits=wavread('yuyin.wav',R);k=0:R-1;wavplay(y,fs);d=(rand(R,2)-0.5)*0.3;x=y+d;wavplay(x,fs);figure(1);plot(k,d, 'r-.', k,s, 'b-', k,x, 'g-'); xlabel('k'); legend('dk', 'sk', 'xk');M=1;b=ones(M,1)/M;a=1;y=filter(b,a,x);wavplay(

15、y,fs);figure(2);plot(k,s, 'b-', k,y, 'r-'); xlabel('k');legend('sk', 'yk');【仿真结果】零输入响应零状态响应M=1时M=10时【结果分析】1. 首先，我们取了M=5，去噪后的信号与原信号相比有一定的相近，即去噪功能相等。2. 其次，我们取了M=200，此时的操作已经不能称之为去噪，信号有一定的失真。3. 最后，我们取了M=1，去噪后的信号加上噪音的信号基本无区别。可得结论：M太大时，信号会失去失真度；M太小，平均值范围太小与加噪信号区别不大。

16、【自主学习内容】Lsim函数的运用，filter函数的运用，对正常信号的加噪处理（即randn函数的运用）还有运用M点滑动平均系统对噪声信号进行去噪【阅读文献】MATLAB释义与实现(第二版)【发现问题】 (专题研讨或相关知识点学习中发现的问题)：axis(0 8 -0.5 2.5错误: 表达式或语句不正确-可能 (、 或 不对称。是不是想输入:>> axis(0 8 -0.5 2.5)当开始时使用rand函数加噪时，令 d=0.01*rand(1,R)*0.3-0.5;结果提示? Error using => plusMatrix dimensions must agree

17、.意思d和x的维数没有保持一致。改用randn函数时，令d=0.01*randn(size(x)*0.3-0.5就成功地把噪声加入了进去。【问题探究】 在对两个函数进行加法计算时要保持维数相同。【研讨内容】题目2：连续信号卷积的近似计算两个连续信号的卷积定义为为了进行数值计算，需对连续信号进行抽样。记xk=x(kD), hk=h(kD), D为进行数值计算的抽样间隔。则连续信号卷积可近似的写为(1)这就可以利用conv函数可近似计算连续信号的卷积。设x(t)=u(t)-u(t-1)，h(t)=x(t)*x(t)，(a)为了与近似计算的结果作比较，用解析法求出y(t)=x(t)*h(t)；(b)

18、用不同的D计算出卷积的数值近似值，并和(a)中的结果作比较；(c)证明(1)式成立；(d)若x(t)和h(t)不是时限信号，则用上面的方法进行近似计算会遇到什么问题？给出一种解决问题的方案；(e) 若将x(t)和h(t)近似表示为推导近似计算卷积的算法。取相同的抽样间隔，比较两种方法的计算卷积误差。【题目分析】（a），求出y(t)的函数表达式，并画出它的图像，便于和后面通过算出的信号作比较。（b），Matlab中不可以直接进行连续的信号的卷积，必须得先对连续信号采样，得到离散的信号，然后求得两个离散信号的卷积，得到另外一个离散信号yk，最后再把离散的信号连续化，得到要得到的y(k)。通过改变的

19、值来与（a）中所得结果进行对比。（c），由于积分可以用求和的极限来表示，所以，当足够小时，成立。（d），Matlab不能表示出非时限信号。我们可以用符号运算或者（e），两个相同的宽度的矩形波的卷积为一个三角波。根据卷积性质，推导出卷积的近似算法。【仿真程序】(a)为了与近似计算的结果作比较，用解析法求出y(t)=x(t)* h(t)；dt=0.000001;t=-1:dt:4;y=(0.5*t.2).*t>0&t<=1+(-t.2+3*t-1.5).* t>1&t<=2 +(0.5*t.2-3*t+4.5).*t>2&t<=3;plo

20、t(t,y,'r');title('解析法求出y（t）的波形');clear all;dt=0.00001;t=-1:dt:4;x=rectpuls(t-0.5,1);axis(-1,3,-2,2);hold;plot(t,x,'b');title('x(t)的波形');h=tripuls(t-1,2);plot(t,h,'g');title('h(t)的波形');clear all;(b)用不同的D计算出卷积的数值近似值，并和(a)中的结果作比较；clear all;c1=0.01,0.1,0.2

21、5,0.3,0.5;for n=1:length(c1) dt=c1(n);t=-0.5:dt:2;x=rectpuls(t-0.5,1);h=tripuls(t-1,2);y=conv(x,h)*dt;t1=-1:dt:4;plot(t1,y)str2=num2str(dt);str3=strcat('当dt=',str2);str1=strcat(str3,'的波形');title(str1);pause;end；(e)若将x(t)和h(t)近似表示为 推导近似计算卷积的算法。取相同的抽样间隔，比较两种方法的计算卷积误差。clear all;dt=0.001

22、;c1=0.01,0.1,0.25,0.3,0.5;for n=1:length(c1)Dt=c1(n);N=Dt/dt;t=0:dt:4;m=length(t);n1=1/Dt;n2=2/Dt;x=rectpuls(t-0.5,1);x1=x(1:N:m);h=tripuls(t-1,2,0);h1=h(1:N:m);sum=0;for m=1:n2 for n=1:n1 y=t>(m+n)*Dt&t<(m+n+1)*Dt.*(t-(m+n)*Dt)+t>(m+n+1)*Dt&t<(m+n+2)*Dt.*(m+n+2)*Dt-t); sum=sum+x

23、1(n)*h1(m)*y; end;end; sum=sum(find (sum=0) y=sum; t=linspace(0,3,length(y);plot(t,y,'r');str2=num2str(Dt);str3=strcat('当Dt=',str2);str1=strcat(str3,'的波形');title(str1);pause;end；【仿真结果】'解析法求出y（t）的波形'x(t)波形：h（t）波形：(b)用不同的D计算出卷积的数值近似值，并和(a)中的结果作比较；用不同的D计算出卷积：(e)若将x(t)和h(

24、t)近似表示为推导近似计算卷积的算法。取相同的抽样间隔，比较两种方法的计算卷积误差。【结果分析】（a）为了与近似计算的结果作比较，用解析法求出y(t)=x(t)*h(t)；x(t)=u(t)-u(t-1)，y(t)=x(t)*h(t)= u(t)-u(t-1)* r(t)-2r(t-1)+r(t-2)当dt越小的时候，图像越平滑，越接近于利用解析法求得的结果。dt越大是，误差越大。(d)若x(t)和h(t)不是时限信号：方法一： matlab软件不能取到无穷大，函数值在到达某一值时会变成0。解决办法是我们可以根据自己的需要，设定x(t)和h(t)在某一范围内的函数值，其他值均为0。方法二：我们

25、可以采用matlab中的符号运算。比如：x1(t)=exp(-abs(t); x2(t)=sin(t);那么x1(t)与x2(t)卷积可以通过以下程序实现：syms t x ;x1(t)=exp(-abs(t); x2(t)=sin(t);%c= x1(t)* x2(t-x);c= exp(-abs(x)* (sin(t-x);y=int(c,x,-inf,inf)求得; x1(t)与 x2(t)的卷积为y=sin(t)(e)经过取值从0.01,0.1,0.25,0.3,0.5两种算法的不同波形比较，我们可以得到（b）的方法比较好，更接近于理论值。【自主学习内容】卷积函数conv的调用格式【阅

26、读文献】1陈后金 胡健 薛健 信号与系统. M 北京：高等教育出版社2肖燕彩 邱成 齐红元 MATLAB语言及实践教程. M 北京：北京交通大学出版社【发现问题】 (专题研讨或相关知识点学习中发现的问题)：暂无。3.信号的频域分析专题研讨【目的】(1) 建立工程应用中有效带宽的概念，了解有限次谐波合成信号及吉伯斯现象。(2) 掌握带限信号，带通信号、未知信号等不同特性的连续时间信号的抽样，以及抽样过程中的参数选择与确定。认识混叠误差，以及减小混叠误差的措施。(3) 加深对信号频域分析基本原理和方法的理解。(4) 锻炼学生综合利用所学理论和技术，分析与解决实际问题的能力。【研讨内容】基础题题目1

27、：吉伯斯现象(1)以定义信号的有效带宽，试确定下图所示信号的有效带宽，取A=1，T=2。(2)画出有效带宽内有限项谐波合成的近似波形，并对结果加以讨论和比较。(3)增加谐波的项数，观察其合成的近似波形，并对结果加以讨论和比较。(a) 周期矩形信号 (b) 周期三角波信号【知识点】连续周期信号的频域分析，有效带宽，吉伯斯现象【信号频谱及有效带宽计算】(a)周期矩形信号： 由于此周期矩形信号的周期T=2，所以 周期信号在区间-1,1的表达式为由于x(t)是奇对称信号，因此有  C0=0根据傅里叶系数计算公式，有 计算可得周期矩形信号的频谱Cn为，由

28、功率计算公式：得：P=1/4。根据有效带宽的定义：，解得。(b)周期三角波信号的频谱及有效带宽的计算： 由于此周期矩形信号的周期T=2，所以 W0=2/T= 周期信号在区间0,2的表达式为由于这个周期三角波信号的偶对称信号，所以根据傅里叶系数的计算 公式，在时有计算可得周期矩形信号的频谱为 由功率计算公式得P=1/3根据有效带宽的定义，得【仿真程序】（1） 周期矩形信号： t=-2:0.001:2; N=input('N='); c0=0; xN=c0*ones(1,length(t);f

29、or n=1:2:N xN=xN+sin(n*pi*t)/(n*pi)-cos(n*pi)*sin(n*pi*t)/(n*pi);end plot(t,xN);Grid on; 周期三角波信号： t=-2:0.001:2;N=input('N='); c0=0.5; xN=c0*ones(1,length(t); for n=1:2:N xN=xN+2*cos(n*pi)*(1-cos(n*pi)*cos(n*pi*t)/(n*n*pi*pi); end pl

30、ot(t,xN);grid on 【仿真结果】（a）  周期矩形信号的仿真结果是：N=3时N=7时N=31时(b) 周期三角波的信号的仿真结果：N=3时N=7时N=31时【结果分析】提示：应从以下几方面对结果进行分析：（1） 随着N值的增加，图(a) 和图(b)信号有效带宽内有限项谐波合成波形与原波形越来越接近。到一定的程度，两个图形基本达到一致。（2） 时域特性与有效带宽内谐波次数的关系。谐波会随着其次数的增加，而减弱其对信号时域的影响。（3） 谐波次数增加，（a）在不连续点附近部分和 x(t)所呈现的起伏,这个起伏的峰值大小似乎不

31、随N 增大而下降，（b）图中也是一样。也就是说，一个不连续信号x(t)的傅里叶级数的截断近似xN(t),一般来说,在接近不连续点处将呈现高频起伏和超量,而且,若在实际情况下利用这样一个近似式的话,就应该选择足够大的N ,以保证这些起伏拥有的总能量可以忽略.当然,在极限情况下,近似误差的能量是零,而且一个不连续的信号(如方波)的傅里叶级数表示是收敛的。【自主学习内容】吉伯斯现象与傅里叶级数的联系，matlab中信号的合成，连续时间周期信号的傅里叶级数的物理意义和分析方法。 【阅读文献】1陈后金 胡健 薛健 信号与系统. M 北京：高等教育出版社2肖燕彩 邱成 齐红元

32、 MATLAB语言及实践教程. M 北京：北京交通大学出版社【发现问题】信号的频谱不会用简单的函数表示，如Sa（t）函数的表示，使得其在matlab中的实现显得简单，容易解得，这样程序也不会显得繁琐。【问题探究】  在一些信号的出来中要尽量的利用简单函数的线性组合来表示，然后利用一些频移性质，来简单的处理这些信号。【研讨内容】中等题题目2：分析音阶的频谱(1) 录制你所喜欢乐器（如钢琴、小提琴等）演奏的音阶，并存为wav格式(可以从网上下载)。(2) 画出各音阶的时域波形，并进行比较。(3) 对所采集的音阶信号进行频谱分析，比较各音阶的频谱。【知识点】连续时间信号的频域分析【温馨提示

33、】利用MATLAB提供的函数fft计算频谱。【题目分析】【仿真程序】y,fs,bits=wavread('D:/qin.wav'); sound(y,fs,bits); m=length(y); Y=fft(y,m); subplot(2,1,1); plot(y) title('y') subplot(2,1,2); plot(abs(Y); title('abs');【仿真结果】【结果分析】提示：应从以下几方面对结果进行分析：(1) 你所选择乐器演奏的音阶，其时域波形的包络有何特点？一开始幅度比较小，渐渐进入主题后振幅有一定加大，其幅度的涨落

34、很有规律(2) 你所选择乐器演奏的音阶，其频谱有何特点？基波是多少？谐波是多少？x,fs,bits=wavread('D:/qin.wav'); sound(x,fs,bits); N=length(x); % x 是待分析的数据 n=1:N;%1-FFT X=fft(x); % FFT X=X(1:N/2); Xabs=abs(X); Xabs(1) = 0; %直流分量置0 for i= 1 : m Amax,index=max(Xabs);if(Xabs(index-1) > Xabs(index+1)a1 = Xabs(index-1) / Xabs(index)

35、; r1 = 1/(1+a1);k01 = index -1; else a1 = Xabs(index) / Xabs(index+1); r1 = 1/(1+a1); k01 = index; end Fn = (k01+r1-1)*fs/N; %基波频率 An = 2*pi*r1*Xabs(k01)/(N*sin(r1*pi); %基波幅值 Pn = phase(X(k01)-pi*r1; %基波相角 单位弧度 Pn = mod(Pn(1),pi); end 【自主学习内容】格式转换；基波分析、谐波分析。【阅读文献】1陈后金 胡健 薛健 信号与系统. M 北京：高等教育出版社2肖燕彩 邱

36、成 齐红元 MATLAB语言及实践教程. M 北京：北京交通大学出版社【发现问题】(1) 改变音阶的包络，相应音阶听起来会有什么变化？ 音色和响度都发生了变化。(2) 音阶频谱中的谐波分量有什么作用？ 影响着音色。(3) 你所分析的乐器各音阶对应的频率是多少，之间存在什么关系？音高和频率是指数关系。f=440×2（p-69)/12 ，其中P是音高，f是频率【研讨内容】拓展题题目3：连续时间信号的抽样(1) 对带限信号（如,等），确定合适的抽样间隔T，分析的频谱和抽样所得到离散信号的频谱X(ejW)，并将两者进行比较。 (2) 将正弦信号按抽样频率fs=8kHz进行1 秒钟抽样，得离散

37、正弦序列xk为比较f0=2kHz, 2.2 kHz, 2.4 kHz, 2.6 kHz和 f0=7.2 kHz, 7.4 kHz, 7.6 kHz, 7.8 kHz 两组信号抽样所得离散序列的声音，解释所出现的现象。(3) 对于许多具有带通特性的信号，举例验证可否不需要满足?【知识点】连续非周期信号的频谱，离散非周期信号的频谱，时域抽样，频域抽样【温馨提示】 (1) 利用MATLAB提供的函数fft计算抽样所得序列xk的频谱。(2) 利用MATLAB函数 sound(x, fs)播放正弦信号和声音信号。(3) 可以利用仪器或仿真软件产生具有带通特性的信号。(1)【题目分析】 抽样函数Sa(t)

38、的频谱应该为非周期的矩形脉冲，但由于计算机处理时的信号均为离散信号，因此我们只关心其频谱在一个周期内的形状。 因此，一方面，我们决定定义出Sa(t)函数后，采用Matlab的函数fft()来计算频谱。 另一方面，通过信号与系统课程的学习，我们可以从理论上推导出Sa(t)函数的频谱并用计算机画图。 之后将两种方法得到的频谱进行比较，观察其中的差异。 最后，再增大Sa(t)函数时域的抽样间隔T，即减小抽样频率fs,观察时域与频域波形的变化。【仿真程序】fx=1;fs=10;N=150*fs;n=(0:N-1)-N/2;t=n/fs;y=sinc(fx*t/pi);subplot(111);plot

39、(t,y);grid onxlim(-5*pi/fx 5*pi/fx);xlabel('t');ylabel('x(t)');XA=abs(fftshift(fft(y,N)/fs;subplot(111);plot(2*pi*n*fs/N,XA);xlabel('w');ylabel('X(jw)');xlim(-2*fx 2*fx);grid ontemp=(pi/fx).*ones(1,fx*20);wth=linspace(-fx,fx,length(temp);subplot(111);plot(wth,temp);ho

40、ld on;plot(-fx,-fx,0,pi/fx)plot(fx,fx,0,pi/fx)xlabel('w');ylabel('X(jw)');xlim(-2*fx 2*fx);grid onfx=1;fs=0.3;N=1500;n=(0:N-1)-N/2;t=n/fs;y=sinc(fx*t/pi);figure(2);subplot(111);plot(t,y);xlim(-5*pi/fx 5*pi/fx);xlabel('t');ylabel('x(t)');grid onXA=abs(fftshift(fft(y,N)

41、/fs;subplot(111);plot(2*pi*n*fs/N,XA);xlabel('w');ylabel('X(jw)');axis(-2*fx 2*fx 0 8);grid on【仿真结果】抽样信号Sa(t)抽样信号Sa(t)的频谱抽样信号Sa(t)的频谱（理论）增大抽样间隔T的抽样信号Sa(t)增大抽样间隔T的抽样信号Sa(t)的频谱【结果分析】可以看出，当增大采样间隔T到一定值后，信号频谱将发生混叠，时域波形也发生较大失真。频域混叠后的幅度值为原来不混叠时的两倍，可知是相邻周期间的频谱发生了叠加。经过我们试验，得到当采样间隔约为1/pi时，为临界采

42、样，及临界采样频率为pi。另外，我们发现Sa(t)函数的频域波形有非常明显的关于稳定值的震荡，可知其为吉伯斯现象。增大取样点数N后，函数在跳变点处仍然有震荡，即吉伯斯现象不消失。（2）【题目分析】由Nyquist抽样定理，将正弦信号以大于其频率二倍的抽样频率抽样时，信号可以被较好还原，而当抽样频率小于正弦频率二倍抽样时，信号将发生失真。【仿真程序】%产生指定频率的正弦声音信号并按固定的抽样频率播放f0=2000 2200 2400 2600 7200 7400 7600 7800;fs=8000;N=2*fs;n=(0:N-1);t=n/fs;for i=1:8y=sin(2*pi*f0(i)

43、*t);%定义声音信号XA=abs(fftshift(fft(y,N)/N;%计算幅度频谱f=n*fs/N;f=f-max(f)/2;sound(y,fs);subplot(211);plot(t,y);%画出时域波形xlabel('t');ylabel('x(t)');grid on;title('频率f=',num2str(f0(i),'Hz时域波形')xlim(0 0.02);subplot(212);plot(f,XA);%画出幅度频谱xlabel('f');ylabel('X(jw)');

44、grid on;title('频率f=',num2str(f0(i),'Hz频域波形')pause(1);end【仿真结果】我们选取2000Hz和7200Hz的两个正弦信号以抽样频率fs=8000Hz抽样后的时域波形和频谱进行对比，仿真后画出的图形如下所示：【结果分析】如我们所料，7200Hz的正弦信号在以抽样频率fs=8000Hz抽样后，发生了严重的失真，由于其周期性而使得频率变为800Hz。这说明，抽样频率低于正弦信号最大频率的两倍抽样时，信号会发生失真。这种声音信号实际频率随正弦信号频率与抽样频率关系而变化的现象引起了我们的兴趣，于是我们进行了以下仿真来探

45、索这一现象。%用来生成一个抽样频率固定，频率随时间变化的正弦声音信号fs=8000;N=0.5*fs;n=0:N-1;t=n/fs;y=sin(2*pi*0*t);for f=50:50:8000 %频率从50Hz到8000Hz间隔50Hz变化y=y sin(2*pi*f*t);endsound(y,fs)运行以上仿真程序，我们听到实际声音的频率，首先随定义的正弦信号的频率升高而升高，到4000Hz时不可闻，此时时域波形的幅度全为零。之后，实际声音的频率，随定义的正弦信号的频率升高而降低，到接近8000Hz时，音高最低。（3）【题目分析】对于许多具有带通特性的信号，其带宽远小于2fH,一般处理

46、带通信号时常遵循带通信号采样定理，即采样频率fs>2(fH-fL)即可。（其中，fH为信号的最大频率，fL为信号的最大频率）经讨论，我们决定用matlab生成一个有一定频谱特点的带限声音信号，然后对这个信号进行抽样。通过观察抽样前后的频谱变化，来验证以上分析的正确性。【仿真程序】fs=44100;N=0.5*fs;n=0:N-1;t=n/fs;y=sin(2*pi*1950*t);for f=2002:2:2100;y=y sin(2*pi*f*t);endfor f=2150:2:2250;y=y sin(2*pi*f*t);endwavwrite(y,fs,'C:Usersi

47、Desktop信号研讨图片文件钉宫理惠.wav');cleary,fs=wavread('C:UsersiDesktop信号研讨图片文件钉宫理惠.wav');y=y(:,1);N=length(y);n=0:N-1;f=n*fs/N;Xabs=abs(fft(y,N)/N;Xabs=Xabs/max(Xabs);subplot(211);plot(f(1:N/2),Xabs(1:N/2);xlabel('f/Hz');ylabel('X(jw)');grid on;axis(0 3000 0 1.1);fsc=646;Nc=length(

48、y)/fs*fsc;B=N/Nc;for i=1:Nc yc(i,:)=y(round(i*B),:);endNc=length(yc);nc=0:Nc-1;fc=nc*fsc/Nc;Xabsc=abs(fft(yc,Nc)/Nc;Xabsc=Xabsc/max(Xabsc);subplot(212);plot(fc(1:Nc/2),Xabsc(1:Nc/2);xlabel('f/Hz');ylabel('X(jw)');grid on;axis(0 800 0 1.1);【仿真结果】【结果分析】 我们生成的信号的频率从1950Hz到2250Hz之间变化，按我们

49、学过的低通信号的Nyquist抽样定理，此时采样频率fsc必须至少大于2250*2=4500Hz，才能保证信号失真较小。但是，我们发现，当我们取采样频率为646Hz时，画出的抽样后声音信号频域波形，在形状上与原信号频域波形相差不大，这说明，我们我们可以从646Hz抽样后的信号中重构得到我们需要的原始信息中的大部分，并且失真较小。这样，便验证了带通信号采样定理的正确性。【自主学习内容】连续信号的采样及重构傅里叶快速算法fft的功能及使用方法带通信号采样定理【阅读文献】【发现问题】 若连续时间信号的最高频率未知，该如何确定对信号进行抽样的最大间隔？ 一般通过查询资料，确定自己所处理类型的信号的最大

50、频率fm，然后，让信号通过一个截止频率为fm的低通滤波器，再依据抽样定理进行抽样。【问题探究】带通信号抽样频率确定的理论分析。【研讨内容】扩展题题目4：男生女生声音信号的转换(1) 采集wav格式的男女生语音信号。(2) 对所采集的男女生信号进行频谱分析，比较男女生频谱的特点。(3) 实现男生女生声音信号的转换。【知识点】连续非周期信号的频谱，离散非周期信号的频谱，时域抽样，频域抽样【温馨提示】可以根据傅里叶变换公式，利用数值积分计算；也可以利用MATLAB提供的函数fft计算。【男生女生声音信号转换的方法】(1)更改基频。不同人的基频不同统计如下：正常成年男声：0200Hz；正常成年女声：2

51、00450Hz；小孩声音的基频要比女声的高，老年人的基频要比男声的低。过整理统计可知女声基频=男声基频*1.5。 本程序使用的是通过抽样与插值的方式来达到基频的改变。以女变男为例：用整数D 对语音信号X（n）进行抽取 Xd=X（Dn）；然后X（n）的抽样频率提高到 I（整数）倍，即为对X（n）的插值。(2)时长规整通过抽样插值来改变基频 也使播放速 度，播放时间发生改变，因此 通过时 长规整的方式来使播放速度和时间恢 复到原来。 本程序使用的是用重叠叠加算法来达到 时长规整。 重叠叠加

52、算法原理： 它分为两个阶段分解和合成 将原始信号以帧长 N，帧间距sa 进行分解，然后以帧间距 ss 进行合成。 sa 与ss 的的比值决定了时长规整因子 F=sa/ss。【仿真程序】（2）x,fs,bits=wavread('a.wav')  x=x(:,1);  sigLength=length(x) t=(0:sigLength-1)/fs; y=fft(x); y=fftshift(y); plot(t,y)&#

53、160; xlabel('Time(s)'); sound(x,fs)    x,fs,bits=wavread('b.wav') x=x(:,1);  sigLength=length(x) t=(0:sigLength-1)/fs; y=fft(x); y=fftshift(y); plot(t,y)  xlabel('Time(s)'); sound(x,fs)  其中，a为女生音频，b为男生音频。 （

54、3）女变男： function Y=voical(x)  x,fs,bits=wavread('a.wav') d=resample(x,3,2) W=400; Wov=W/2; Kmax=2*W; Wsim=Wov; xdecim=8; kdecim=2; X=d' F=1.5;  Ss=W-Wov; xpts=size(X,2); ypts=round(xpts/F); Y=zeros(1,ypts);

55、60; xfwin=(1:Wov)/(Wov+1); ovix=(1-Wov):0; newix=1:(W-Wov);  simix=(1:xdecim:Wsim)-Wsim;  padX=zeros(1,Wsim),X,zeros(1,Kmax+W-Wov); Y(1:Wsim)=X(1:Wsim); xabs=0; lastxpos=0; km=0;  for ypos=Wsim:Ss:(ypts-W);  xpos=F*ypos;  kmpred=km+(xp

56、os-lastxpos); lastxpos=xpos; if(kmpred<=Kmax) km=kmpred; else  ysim=Y(ypos+simix); rxy=zeros(1,Kmax+1); rxx=zeros(1,Kmax+1); Kmin=0;  for k=Kmin:kdecim:Kmax  xsim=padX(Wsim+xpos+k+simix); rxx(k+1)=norm(xsim); rxy=(ysim*xsim'); end  Rxy=(rxx=0).*rxy./(rxx+(rxx=0); km=min(find(Rxy=max(Rxy)-1);  end  xabs=xpos+km;  Y(ypos+ovix)=(1-xfwin).*Y(ypos+ovix)+(xfwin.*padX(Wsim+xabs+ovix); Y(ypos+newix)=padX(Wsim+xab