北京交通大學(xué)《信號(hào)與系統(tǒng)》研究性學(xué)習(xí)實(shí)驗(yàn)報(bào)告

1、信號(hào)與系統(tǒng)課程研究性學(xué)習(xí)手冊(cè)姓名 學(xué)號(hào) 同組成員 14* 14* 14* 指導(dǎo)教師 時(shí)間 1. 信號(hào)的時(shí)域分析專題研討【目的】(1) 掌握基本信號(hào)及其特性，了解實(shí)際信號(hào)的建模。(2) 掌握基本信號(hào)的運(yùn)算，加深對(duì)信號(hào)時(shí)域分析基本原理和方法的理解，并建立時(shí)頻之間的感性認(rèn)識(shí)。(3) 學(xué)會(huì)仿真軟件MATLAB的初步使用方法，掌握利用MATLAB進(jìn)行信號(hào)表示和信號(hào)運(yùn)算?！狙杏憙?nèi)容】題目1：基本信號(hào)的產(chǎn)生，語(yǔ)音的讀取與播放1） 生成一個(gè)正弦信號(hào)，改變正弦信號(hào)的角頻率和初始相位，觀察波形變化，并聽其聲音的變化。2） 生成一個(gè)幅度為1、基頻為2Hz、占空比為50%的周期方波。3） 觀察一定時(shí)期內(nèi)的股票上證指數(shù)

2、變化，生成模擬其變化的指數(shù)信號(hào)。4） 分別錄制一段男聲、女聲信號(hào)，進(jìn)行音頻信號(hào)的讀取與播放，畫出其時(shí)域波形?！緶剀疤崾尽?1)利用MATLAB函數(shù) wavread(file)讀取.wav格式文件。(2)利用MATLAB函數(shù) sound(x, fs)播放正弦信號(hào)和聲音信號(hào)。【題目分析】【仿真程序】1）生成一個(gè)正弦信號(hào)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）生成一個(gè)幅度為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）觀察一定時(shí)期內(nèi)的股票上證指數(shù)變化，生成模擬其變化的指數(shù)信號(hào)。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）分別錄制一段男聲、女聲信號(hào)，進(jìn)行音頻信號(hào)的讀取與播放，畫出其時(shí)域波形。y,

4、x,nbits=wavread('C:UsersiDesktop信號(hào)研討圖片文件日野聰.wav');sound(y,x); plot(y)y,x,nbits=wavread('C:UsersiDesktop信號(hào)研討圖片文件釘宮理惠.wav');sound(y,x); plot(y)【仿真結(jié)果】1）生成一個(gè)正弦信號(hào) 改變其角頻率和初始相位2）生成一個(gè)幅度為1、基頻為2Hz、占空比為50%的周期方波3）觀察一定時(shí)期內(nèi)的股票上證指數(shù)變化，生成模擬其變化的指數(shù)信號(hào)。4）分別錄制一段男聲、女聲信號(hào)，進(jìn)行音頻信號(hào)的讀取與播放，畫出其時(shí)域波形?！窘Y(jié)果分析】(1) 隨著正弦信號(hào)

5、角頻率的變化，其波形變得更緊密，聲音逐漸變得尖細(xì)而高。(2) 男聲多低沉粗獷，主要是因?yàn)槟新曋械皖l分量更多；女聲多高亢清脆，主要是因?yàn)榕曋懈哳l分量更多。【自主學(xué)習(xí)內(nèi)容】wavread的讀取與sound函數(shù)的使用，音頻文件的放置?！鹃喿x文獻(xiàn)】MATLAB語(yǔ)言與實(shí)踐教程（第2版）、信號(hào)與系統(tǒng)【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識(shí)點(diǎn)學(xué)習(xí)中發(fā)現(xiàn)的問題)：根據(jù)聲音信號(hào)的什么特征能有效區(qū)分出男聲和女聲？【問題探究】 高頻低頻成分的多少，以及聲音的尖細(xì)程度【研討內(nèi)容】題目2：信號(hào)的基本運(yùn)算（語(yǔ)音信號(hào)的翻轉(zhuǎn)、展縮）1） 將原始音頻信號(hào)在時(shí)域上進(jìn)行延展、壓縮，2） 將原始音頻信號(hào)在時(shí)域上進(jìn)行幅度放大與縮小，3）

6、將原始音頻信號(hào)在時(shí)域上進(jìn)行翻轉(zhuǎn)，【題目分析】【仿真程序】 原始音頻信號(hào)y,x,nbits=wavread('D:音樂動(dòng)畫歌曲（207曲）灼眼的夏娜無(wú)損 17曲Light My Fire.wav');plot(y)1） 將原始音頻信號(hào)在時(shí)域上進(jìn)行延展、壓縮 將原始音頻信號(hào)在時(shí)域上延展3倍y,x,nbits=wavread('D:音樂動(dòng)畫歌曲（207曲）灼眼的夏娜無(wú)損 17曲Light My Fire.wav');sigLength=length(y);t=3*(0:sigLength-1)/x;figure;plot(t,y);xlabel('Time(s

7、)'); 將原始音頻信號(hào)在時(shí)域上壓縮為1/5y,x,nbits=wavread('D:音樂動(dòng)畫歌曲（207曲）灼眼的夏娜無(wú)損 17曲Light My Fire.wav');sigLength=length(y);t=1/5*(0:sigLength-1)/x;figure;plot(t,y);xlabel('Time(s)');2）將原始音頻信號(hào)在時(shí)域上進(jìn)行幅度放大與縮小 將原始音頻信號(hào)的幅度擴(kuò)大為2倍y,x,nbits=wavread('D:音樂動(dòng)畫歌曲（207曲）灼眼的夏娜無(wú)損 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)'); 將原始音頻信號(hào)的幅度縮小為1/2y,x,nbits=wavread('D:音樂動(dòng)畫歌曲（207曲）灼眼的夏娜無(wú)損 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）將原始音頻信號(hào)在時(shí)域上進(jìn)行翻轉(zhuǎn)y,x,nbits=wavread('D:音樂動(dòng)畫歌曲（207曲）灼眼的夏娜無(wú)損 17曲Light My Fire.wav');sigLength=length(y);t=-1*(0:sigLength-1)/x;figure

10、;plot(t,y);xlabel('Time(s)');【仿真結(jié)果】 原始音頻信號(hào)1）將原始音頻信號(hào)在時(shí)域上進(jìn)行延展、壓縮 將原始音頻信號(hào)在時(shí)域上延展3倍 將原始音頻信號(hào)在時(shí)域上壓縮為1/52）將原始音頻信號(hào)在時(shí)域上進(jìn)行幅度放大與縮小 將原始音頻信號(hào)的幅度擴(kuò)大為2倍 將原始音頻信號(hào)的幅度縮小為1/23）將原始音頻信號(hào)在時(shí)域上進(jìn)行翻轉(zhuǎn)【結(jié)果分析】以時(shí)間為橫軸的圖比以頻率為橫軸的圖看起來復(fù)雜得多，時(shí)域上的分析在現(xiàn)階段很困難。對(duì)于原始音頻信號(hào)，時(shí)域上的分析太復(fù)雜，而如果在頻域上分析，則簡(jiǎn)化很多?！咀灾鲗W(xué)習(xí)內(nèi)容】信號(hào)的獲取以及如何使用；MATLAB的一些基本語(yǔ)法；如何根據(jù)MATLAB

11、語(yǔ)法來編寫自己想表達(dá)的語(yǔ)句?！鹃喿x文獻(xiàn)】MATLAB語(yǔ)言與實(shí)踐教程（第2版）、信號(hào)與系統(tǒng)【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識(shí)點(diǎn)學(xué)習(xí)中發(fā)現(xiàn)的問題)：在查找資料時(shí)，發(fā)現(xiàn)音頻單聲道和雙聲道的編寫此程序會(huì)有不同，同時(shí)語(yǔ)音信號(hào)的錄入需要我們的注意。【問題探究】 2.系統(tǒng)的時(shí)域分析專題研討【目的】(1) 掌握系統(tǒng)響應(yīng)的時(shí)域求解，加深對(duì)系統(tǒng)時(shí)域分析基本原理和方法的理解。(2) 掌握連續(xù)系統(tǒng)零狀態(tài)響應(yīng)（卷積積分）數(shù)值計(jì)算的方法。(3) 學(xué)會(huì)仿真軟件MATLAB的初步使用方法，掌握利用MATLAB求解連續(xù)系統(tǒng)和離散系統(tǒng)的零狀態(tài)響應(yīng)。(4) 培養(yǎng)學(xué)生發(fā)現(xiàn)問題、分析問題和解決問題的能力?！狙杏憙?nèi)容】題目1：系統(tǒng)響應(yīng)時(shí)

12、域求解1） 求一個(gè)RLC電路的零輸入響應(yīng)和零狀態(tài)響應(yīng)，2） 將原始音頻信號(hào)中混入噪聲，然后用M點(diǎn)滑動(dòng)平均系統(tǒng)對(duì)受噪聲干擾的信號(hào)去噪，改變M點(diǎn)數(shù)，比較不同點(diǎn)數(shù)下的去噪效果，【題目分析】（1） 一個(gè)RLC電路，若,L=1H,C=1F,電容上的初始儲(chǔ)能為Vc=2V，電感初始儲(chǔ)能為i=1A，試求輸入激勵(lì)為X（t）時(shí)的零輸入響應(yīng)和零狀態(tài)響應(yīng)。（2） 題目要求采用M點(diǎn)滑動(dòng)平均系統(tǒng)進(jìn)行去噪。M點(diǎn)滑動(dòng)平均系統(tǒng)可以看成是N=0的差分方程。調(diào)用filter函數(shù)時(shí)，調(diào)用參數(shù)a-1=1,b為有M個(gè)元素的向量，b中每個(gè)元素的值為1/M。即M點(diǎn)的滑動(dòng)平均系統(tǒng)輸入輸出關(guān)系為：，同時(shí)我們將噪聲設(shè)為n,函數(shù)為n=rand(n,

13、1);原始信號(hào)為s。通過調(diào)整M值，觀察和比較去噪效果，從而得出結(jié)論?！痉抡娉绦颉浚?）零輸入響應(yīng) 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）零狀態(tài)響應(yīng)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');【仿真結(jié)果】零輸入響應(yīng)零狀態(tài)響應(yīng)M=1時(shí)M=10時(shí)【結(jié)果分析】1. 首先，我們?nèi)×薓=5，去噪后的信號(hào)與原信號(hào)相比有一定的相近，即去噪功能相等。2. 其次，我們?nèi)×薓=200，此時(shí)的操作已經(jīng)不能稱之為去噪，信號(hào)有一定的失真。3. 最后，我們?nèi)×薓=1，去噪后的信號(hào)加上噪音的信號(hào)基本無(wú)區(qū)別。可得結(jié)論：M太大時(shí)，信號(hào)會(huì)失去失真度；M太小，平均值范圍太小與加噪信號(hào)區(qū)別不大。

16、【自主學(xué)習(xí)內(nèi)容】Lsim函數(shù)的運(yùn)用，filter函數(shù)的運(yùn)用，對(duì)正常信號(hào)的加噪處理（即randn函數(shù)的運(yùn)用）還有運(yùn)用M點(diǎn)滑動(dòng)平均系統(tǒng)對(duì)噪聲信號(hào)進(jìn)行去噪【閱讀文獻(xiàn)】MATLAB釋義與實(shí)現(xiàn)(第二版)【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識(shí)點(diǎn)學(xué)習(xí)中發(fā)現(xiàn)的問題)：axis(0 8 -0.5 2.5錯(cuò)誤: 表達(dá)式或語(yǔ)句不正確-可能 (、 或 不對(duì)稱。是不是想輸入:>> axis(0 8 -0.5 2.5)當(dāng)開始時(shí)使用rand函數(shù)加噪時(shí)，令 d=0.01*rand(1,R)*0.3-0.5;結(jié)果提示? Error using => plusMatrix dimensions must agree

17、.意思d和x的維數(shù)沒有保持一致。改用randn函數(shù)時(shí)，令d=0.01*randn(size(x)*0.3-0.5就成功地把噪聲加入了進(jìn)去。【問題探究】 在對(duì)兩個(gè)函數(shù)進(jìn)行加法計(jì)算時(shí)要保持維數(shù)相同?！狙杏憙?nèi)容】題目2：連續(xù)信號(hào)卷積的近似計(jì)算兩個(gè)連續(xù)信號(hào)的卷積定義為為了進(jìn)行數(shù)值計(jì)算，需對(duì)連續(xù)信號(hào)進(jìn)行抽樣。記xk=x(kD), hk=h(kD), D為進(jìn)行數(shù)值計(jì)算的抽樣間隔。則連續(xù)信號(hào)卷積可近似的寫為(1)這就可以利用conv函數(shù)可近似計(jì)算連續(xù)信號(hào)的卷積。設(shè)x(t)=u(t)-u(t-1)，h(t)=x(t)*x(t)，(a)為了與近似計(jì)算的結(jié)果作比較，用解析法求出y(t)=x(t)*h(t)；(b)

18、用不同的D計(jì)算出卷積的數(shù)值近似值，并和(a)中的結(jié)果作比較；(c)證明(1)式成立；(d)若x(t)和h(t)不是時(shí)限信號(hào)，則用上面的方法進(jìn)行近似計(jì)算會(huì)遇到什么問題？給出一種解決問題的方案；(e) 若將x(t)和h(t)近似表示為推導(dǎo)近似計(jì)算卷積的算法。取相同的抽樣間隔，比較兩種方法的計(jì)算卷積誤差?！绢}目分析】（a），求出y(t)的函數(shù)表達(dá)式，并畫出它的圖像，便于和后面通過算出的信號(hào)作比較。（b），Matlab中不可以直接進(jìn)行連續(xù)的信號(hào)的卷積，必須得先對(duì)連續(xù)信號(hào)采樣，得到離散的信號(hào)，然后求得兩個(gè)離散信號(hào)的卷積，得到另外一個(gè)離散信號(hào)yk，最后再把離散的信號(hào)連續(xù)化，得到要得到的y(k)。通過改變的

19、值來與（a）中所得結(jié)果進(jìn)行對(duì)比。（c），由于積分可以用求和的極限來表示，所以，當(dāng)足夠小時(shí)，成立。（d），Matlab不能表示出非時(shí)限信號(hào)。我們可以用符號(hào)運(yùn)算或者（e），兩個(gè)相同的寬度的矩形波的卷積為一個(gè)三角波。根據(jù)卷積性質(zhì)，推導(dǎo)出卷積的近似算法。【仿真程序】(a)為了與近似計(jì)算的結(jié)果作比較，用解析法求出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計(jì)算出卷積的數(shù)值近似值，并和(a)中的結(jié)果作比較；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('當(dāng)dt=',str2);str1=strcat(str3,'的波形');title(str1);pause;end；(e)若將x(t)和h(t)近似表示為 推導(dǎo)近似計(jì)算卷積的算法。取相同的抽樣間隔，比較兩種方法的計(jì)算卷積誤差。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('當(dāng)Dt=',str2);str1=strcat(str3,'的波形');title(str1);pause;end；【仿真結(jié)果】'解析法求出y（t）的波形'x(t)波形：h（t）波形：(b)用不同的D計(jì)算出卷積的數(shù)值近似值，并和(a)中的結(jié)果作比較；用不同的D計(jì)算出卷積：(e)若將x(t)和h(

24、t)近似表示為推導(dǎo)近似計(jì)算卷積的算法。取相同的抽樣間隔，比較兩種方法的計(jì)算卷積誤差?！窘Y(jié)果分析】（a）為了與近似計(jì)算的結(jié)果作比較，用解析法求出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)當(dāng)dt越小的時(shí)候，圖像越平滑，越接近于利用解析法求得的結(jié)果。dt越大是，誤差越大。(d)若x(t)和h(t)不是時(shí)限信號(hào)：方法一： matlab軟件不能取到無(wú)窮大，函數(shù)值在到達(dá)某一值時(shí)會(huì)變成0。解決辦法是我們可以根據(jù)自己的需要，設(shè)定x(t)和h(t)在某一范圍內(nèi)的函數(shù)值，其他值均為0。方法二：我們

25、可以采用matlab中的符號(hào)運(yùn)算。比如：x1(t)=exp(-abs(t); x2(t)=sin(t);那么x1(t)與x2(t)卷積可以通過以下程序?qū)崿F(xiàn)：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)經(jīng)過取值從0.01,0.1,0.25,0.3,0.5兩種算法的不同波形比較，我們可以得到（b）的方法比較好，更接近于理論值。【自主學(xué)習(xí)內(nèi)容】卷積函數(shù)conv的調(diào)用格式【閱

26、讀文獻(xiàn)】1陳后金 胡健 薛健 信號(hào)與系統(tǒng). M 北京：高等教育出版社2肖燕彩 邱成 齊紅元 MATLAB語(yǔ)言及實(shí)踐教程. M 北京：北京交通大學(xué)出版社【發(fā)現(xiàn)問題】 (專題研討或相關(guān)知識(shí)點(diǎn)學(xué)習(xí)中發(fā)現(xiàn)的問題)：暫無(wú)。3.信號(hào)的頻域分析專題研討【目的】(1) 建立工程應(yīng)用中有效帶寬的概念，了解有限次諧波合成信號(hào)及吉伯斯現(xiàn)象。(2) 掌握帶限信號(hào)，帶通信號(hào)、未知信號(hào)等不同特性的連續(xù)時(shí)間信號(hào)的抽樣，以及抽樣過程中的參數(shù)選擇與確定。認(rèn)識(shí)混疊誤差，以及減小混疊誤差的措施。(3) 加深對(duì)信號(hào)頻域分析基本原理和方法的理解。(4) 鍛煉學(xué)生綜合利用所學(xué)理論和技術(shù)，分析與解決實(shí)際問題的能力。【研討內(nèi)容】基礎(chǔ)題題目1

27、：吉伯斯現(xiàn)象(1)以定義信號(hào)的有效帶寬，試確定下圖所示信號(hào)的有效帶寬，取A=1，T=2。(2)畫出有效帶寬內(nèi)有限項(xiàng)諧波合成的近似波形，并對(duì)結(jié)果加以討論和比較。(3)增加諧波的項(xiàng)數(shù)，觀察其合成的近似波形，并對(duì)結(jié)果加以討論和比較。(a) 周期矩形信號(hào) (b) 周期三角波信號(hào)【知識(shí)點(diǎn)】連續(xù)周期信號(hào)的頻域分析，有效帶寬，吉伯斯現(xiàn)象【信號(hào)頻譜及有效帶寬計(jì)算】(a)周期矩形信號(hào)： 由于此周期矩形信號(hào)的周期T=2，所以 周期信號(hào)在區(qū)間-1,1的表達(dá)式為由于x(t)是奇對(duì)稱信號(hào)，因此有  C0=0根據(jù)傅里葉系數(shù)計(jì)算公式，有 計(jì)算可得周期矩形信號(hào)的頻譜Cn為，由

28、功率計(jì)算公式：得：P=1/4。根據(jù)有效帶寬的定義：，解得。(b)周期三角波信號(hào)的頻譜及有效帶寬的計(jì)算： 由于此周期矩形信號(hào)的周期T=2，所以 W0=2/T= 周期信號(hào)在區(qū)間0,2的表達(dá)式為由于這個(gè)周期三角波信號(hào)的偶對(duì)稱信號(hào)，所以根據(jù)傅里葉系數(shù)的計(jì)算 公式，在時(shí)有計(jì)算可得周期矩形信號(hào)的頻譜為 由功率計(jì)算公式得P=1/3根據(jù)有效帶寬的定義，得【仿真程序】（1） 周期矩形信號(hào)： 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; 周期三角波信號(hào)： 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 【仿真結(jié)果】（a）  周期矩形信號(hào)的仿真結(jié)果是：N=3時(shí)N=7時(shí)N=31時(shí)(b) 周期三角波的信號(hào)的仿真結(jié)果：N=3時(shí)N=7時(shí)N=31時(shí)【結(jié)果分析】提示：應(yīng)從以下幾方面對(duì)結(jié)果進(jìn)行分析：（1） 隨著N值的增加，圖(a) 和圖(b)信號(hào)有效帶寬內(nèi)有限項(xiàng)諧波合成波形與原波形越來越接近。到一定的程度，兩個(gè)圖形基本達(dá)到一致。（2） 時(shí)域特性與有效帶寬內(nèi)諧波次數(shù)的關(guān)系。諧波會(huì)隨著其次數(shù)的增加，而減弱其對(duì)信號(hào)時(shí)域的影響。（3） 諧波次數(shù)增加，（a）在不連續(xù)點(diǎn)附近部分和 x(t)所呈現(xiàn)的起伏,這個(gè)起伏的峰值大小似乎不

31、隨N 增大而下降，（b）圖中也是一樣。也就是說，一個(gè)不連續(xù)信號(hào)x(t)的傅里葉級(jí)數(shù)的截?cái)嘟苮N(t),一般來說,在接近不連續(xù)點(diǎn)處將呈現(xiàn)高頻起伏和超量,而且,若在實(shí)際情況下利用這樣一個(gè)近似式的話,就應(yīng)該選擇足夠大的N ,以保證這些起伏擁有的總能量可以忽略.當(dāng)然,在極限情況下,近似誤差的能量是零,而且一個(gè)不連續(xù)的信號(hào)(如方波)的傅里葉級(jí)數(shù)表示是收斂的?！咀灾鲗W(xué)習(xí)內(nèi)容】吉伯斯現(xiàn)象與傅里葉級(jí)數(shù)的聯(lián)系，matlab中信號(hào)的合成，連續(xù)時(shí)間周期信號(hào)的傅里葉級(jí)數(shù)的物理意義和分析方法。 【閱讀文獻(xiàn)】1陳后金 胡健 薛健 信號(hào)與系統(tǒng). M 北京：高等教育出版社2肖燕彩 邱成 齊紅元

32、 MATLAB語(yǔ)言及實(shí)踐教程. M 北京：北京交通大學(xué)出版社【發(fā)現(xiàn)問題】信號(hào)的頻譜不會(huì)用簡(jiǎn)單的函數(shù)表示，如Sa（t）函數(shù)的表示，使得其在matlab中的實(shí)現(xiàn)顯得簡(jiǎn)單，容易解得，這樣程序也不會(huì)顯得繁瑣?！締栴}探究】  在一些信號(hào)的出來中要盡量的利用簡(jiǎn)單函數(shù)的線性組合來表示，然后利用一些頻移性質(zhì)，來簡(jiǎn)單的處理這些信號(hào)?！狙杏憙?nèi)容】中等題題目2：分析音階的頻譜(1) 錄制你所喜歡樂器（如鋼琴、小提琴等）演奏的音階，并存為wav格式(可以從網(wǎng)上下載)。(2) 畫出各音階的時(shí)域波形，并進(jìn)行比較。(3) 對(duì)所采集的音階信號(hào)進(jìn)行頻譜分析，比較各音階的頻譜。【知識(shí)點(diǎn)】連續(xù)時(shí)間信號(hào)的頻域分析【溫馨提示

33、】利用MATLAB提供的函數(shù)fft計(jì)算頻譜?！绢}目分析】【仿真程序】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');【仿真結(jié)果】【結(jié)果分析】提示：應(yīng)從以下幾方面對(duì)結(jié)果進(jìn)行分析：(1) 你所選擇樂器演奏的音階，其時(shí)域波形的包絡(luò)有何特點(diǎn)？一開始幅度比較小，漸漸進(jìn)入主題后振幅有一定加大，其幅度的漲落

34、很有規(guī)律(2) 你所選擇樂器演奏的音階，其頻譜有何特點(diǎn)？基波是多少？諧波是多少？x,fs,bits=wavread('D:/qin.wav'); sound(x,fs,bits); N=length(x); % x 是待分析的數(shù)據(jù) 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 【自主學(xué)習(xí)內(nèi)容】格式轉(zhuǎn)換；基波分析、諧波分析?！鹃喿x文獻(xiàn)】1陳后金 胡健 薛健 信號(hào)與系統(tǒng). M 北京：高等教育出版社2肖燕彩 邱

36、成 齊紅元 MATLAB語(yǔ)言及實(shí)踐教程. M 北京：北京交通大學(xué)出版社【發(fā)現(xiàn)問題】(1) 改變音階的包絡(luò)，相應(yīng)音階聽起來會(huì)有什么變化？ 音色和響度都發(fā)生了變化。(2) 音階頻譜中的諧波分量有什么作用？ 影響著音色。(3) 你所分析的樂器各音階對(duì)應(yīng)的頻率是多少，之間存在什么關(guān)系？音高和頻率是指數(shù)關(guān)系。f=440×2（p-69)/12 ，其中P是音高，f是頻率【研討內(nèi)容】拓展題題目3：連續(xù)時(shí)間信號(hào)的抽樣(1) 對(duì)帶限信號(hào)（如,等），確定合適的抽樣間隔T，分析的頻譜和抽樣所得到離散信號(hào)的頻譜X(ejW)，并將兩者進(jìn)行比較。 (2) 將正弦信號(hào)按抽樣頻率fs=8kHz進(jìn)行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 兩組信號(hào)抽樣所得離散序列的聲音，解釋所出現(xiàn)的現(xiàn)象。(3) 對(duì)于許多具有帶通特性的信號(hào)，舉例驗(yàn)證可否不需要滿足?【知識(shí)點(diǎn)】連續(xù)非周期信號(hào)的頻譜，離散非周期信號(hào)的頻譜，時(shí)域抽樣，頻域抽樣【溫馨提示】 (1) 利用MATLAB提供的函數(shù)fft計(jì)算抽樣所得序列xk的頻譜。(2) 利用MATLAB函數(shù) sound(x, fs)播放正弦信號(hào)和聲音信號(hào)。(3) 可以利用儀器或仿真軟件產(chǎn)生具有帶通特性的信號(hào)。(1)【題目分析】 抽樣函數(shù)Sa(t)

38、的頻譜應(yīng)該為非周期的矩形脈沖，但由于計(jì)算機(jī)處理時(shí)的信號(hào)均為離散信號(hào)，因此我們只關(guān)心其頻譜在一個(gè)周期內(nèi)的形狀。 因此，一方面，我們決定定義出Sa(t)函數(shù)后，采用Matlab的函數(shù)fft()來計(jì)算頻譜。 另一方面，通過信號(hào)與系統(tǒng)課程的學(xué)習(xí)，我們可以從理論上推導(dǎo)出Sa(t)函數(shù)的頻譜并用計(jì)算機(jī)畫圖。 之后將兩種方法得到的頻譜進(jìn)行比較，觀察其中的差異。 最后，再增大Sa(t)函數(shù)時(shí)域的抽樣間隔T，即減小抽樣頻率fs,觀察時(shí)域與頻域波形的變化?！痉抡娉绦颉縡x=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【仿真結(jié)果】抽樣信號(hào)Sa(t)抽樣信號(hào)Sa(t)的頻譜抽樣信號(hào)Sa(t)的頻譜（理論）增大抽樣間隔T的抽樣信號(hào)Sa(t)增大抽樣間隔T的抽樣信號(hào)Sa(t)的頻譜【結(jié)果分析】可以看出，當(dāng)增大采樣間隔T到一定值后，信號(hào)頻譜將發(fā)生混疊，時(shí)域波形也發(fā)生較大失真。頻域混疊后的幅度值為原來不混疊時(shí)的兩倍，可知是相鄰周期間的頻譜發(fā)生了疊加。經(jīng)過我們?cè)囼?yàn)，得到當(dāng)采樣間隔約為1/pi時(shí)，為臨界采

42、樣，及臨界采樣頻率為pi。另外，我們發(fā)現(xiàn)Sa(t)函數(shù)的頻域波形有非常明顯的關(guān)于穩(wěn)定值的震蕩，可知其為吉伯斯現(xiàn)象。增大取樣點(diǎn)數(shù)N后，函數(shù)在跳變點(diǎn)處仍然有震蕩，即吉伯斯現(xiàn)象不消失。（2）【題目分析】由Nyquist抽樣定理，將正弦信號(hào)以大于其頻率二倍的抽樣頻率抽樣時(shí)，信號(hào)可以被較好還原，而當(dāng)抽樣頻率小于正弦頻率二倍抽樣時(shí)，信號(hào)將發(fā)生失真?！痉抡娉绦颉?產(chǎn)生指定頻率的正弦聲音信號(hào)并按固定的抽樣頻率播放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);%定義聲音信號(hào)XA=abs(fftshift(fft(y,N)/N;%計(jì)算幅度頻譜f=n*fs/N;f=f-max(f)/2;sound(y,fs);subplot(211);plot(t,y);%畫出時(shí)域波形xlabel('t');ylabel('x(t)');grid on;title('頻率f=',num2str(f0(i),'Hz時(shí)域波形')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【仿真結(jié)果】我們選取2000Hz和7200Hz的兩個(gè)正弦信號(hào)以抽樣頻率fs=8000Hz抽樣后的時(shí)域波形和頻譜進(jìn)行對(duì)比，仿真后畫出的圖形如下所示：【結(jié)果分析】如我們所料，7200Hz的正弦信號(hào)在以抽樣頻率fs=8000Hz抽樣后，發(fā)生了嚴(yán)重的失真，由于其周期性而使得頻率變?yōu)?00Hz。這說明，抽樣頻率低于正弦信號(hào)最大頻率的兩倍抽樣時(shí)，信號(hào)會(huì)發(fā)生失真。這種聲音信號(hào)實(shí)際頻率隨正弦信號(hào)頻率與抽樣頻率關(guān)系而變化的現(xiàn)象引起了我們的興趣，于是我們進(jìn)行了以下仿真來探

45、索這一現(xiàn)象。%用來生成一個(gè)抽樣頻率固定，頻率隨時(shí)間變化的正弦聲音信號(hào)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)運(yùn)行以上仿真程序，我們聽到實(shí)際聲音的頻率，首先隨定義的正弦信號(hào)的頻率升高而升高，到4000Hz時(shí)不可聞，此時(shí)時(shí)域波形的幅度全為零。之后，實(shí)際聲音的頻率，隨定義的正弦信號(hào)的頻率升高而降低，到接近8000Hz時(shí)，音高最低。（3）【題目分析】對(duì)于許多具有帶通特性的信號(hào)，其帶寬遠(yuǎn)小于2fH,一般處理

46、帶通信號(hào)時(shí)常遵循帶通信號(hào)采樣定理，即采樣頻率fs>2(fH-fL)即可。（其中，fH為信號(hào)的最大頻率，fL為信號(hào)的最大頻率）經(jīng)討論，我們決定用matlab生成一個(gè)有一定頻譜特點(diǎn)的帶限聲音信號(hào)，然后對(duì)這個(gè)信號(hào)進(jìn)行抽樣。通過觀察抽樣前后的頻譜變化，來驗(yàn)證以上分析的正確性。【仿真程序】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信號(hào)研討圖片文件釘宮理惠.wav');cleary,fs=wavread('C:UsersiDesktop信號(hào)研討圖片文件釘宮理惠.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);【仿真結(jié)果】【結(jié)果分析】 我們生成的信號(hào)的頻率從1950Hz到2250Hz之間變化，按我們

49、學(xué)過的低通信號(hào)的Nyquist抽樣定理，此時(shí)采樣頻率fsc必須至少大于2250*2=4500Hz，才能保證信號(hào)失真較小。但是，我們發(fā)現(xiàn)，當(dāng)我們?nèi)〔蓸宇l率為646Hz時(shí)，畫出的抽樣后聲音信號(hào)頻域波形，在形狀上與原信號(hào)頻域波形相差不大，這說明，我們我們可以從646Hz抽樣后的信號(hào)中重構(gòu)得到我們需要的原始信息中的大部分，并且失真較小。這樣，便驗(yàn)證了帶通信號(hào)采樣定理的正確性。【自主學(xué)習(xí)內(nèi)容】連續(xù)信號(hào)的采樣及重構(gòu)傅里葉快速算法fft的功能及使用方法帶通信號(hào)采樣定理【閱讀文獻(xiàn)】【發(fā)現(xiàn)問題】 若連續(xù)時(shí)間信號(hào)的最高頻率未知，該如何確定對(duì)信號(hào)進(jìn)行抽樣的最大間隔？ 一般通過查詢資料，確定自己所處理類型的信號(hào)的最大

50、頻率fm，然后，讓信號(hào)通過一個(gè)截止頻率為fm的低通濾波器，再依據(jù)抽樣定理進(jìn)行抽樣?！締栴}探究】帶通信號(hào)抽樣頻率確定的理論分析?！狙杏憙?nèi)容】擴(kuò)展題題目4：男生女生聲音信號(hào)的轉(zhuǎn)換(1) 采集wav格式的男女生語(yǔ)音信號(hào)。(2) 對(duì)所采集的男女生信號(hào)進(jìn)行頻譜分析，比較男女生頻譜的特點(diǎn)。(3) 實(shí)現(xiàn)男生女生聲音信號(hào)的轉(zhuǎn)換?！局R(shí)點(diǎn)】連續(xù)非周期信號(hào)的頻譜，離散非周期信號(hào)的頻譜，時(shí)域抽樣，頻域抽樣【溫馨提示】可以根據(jù)傅里葉變換公式，利用數(shù)值積分計(jì)算；也可以利用MATLAB提供的函數(shù)fft計(jì)算?！灸猩曇粜盘?hào)轉(zhuǎn)換的方法】(1)更改基頻。不同人的基頻不同統(tǒng)計(jì)如下：正常成年男聲：0200Hz；正常成年女聲：2

51、00450Hz；小孩聲音的基頻要比女聲的高，老年人的基頻要比男聲的低。過整理統(tǒng)計(jì)可知女聲基頻=男聲基頻*1.5。 本程序使用的是通過抽樣與插值的方式來達(dá)到基頻的改變。以女變男為例：用整數(shù)D 對(duì)語(yǔ)音信號(hào)X（n）進(jìn)行抽取 Xd=X（Dn）；然后X（n）的抽樣頻率提高到 I（整數(shù)）倍，即為對(duì)X（n）的插值。(2)時(shí)長(zhǎng)規(guī)整通過抽樣插值來改變基頻 也使播放速 度，播放時(shí)間發(fā)生改變，因此 通過時(shí) 長(zhǎng)規(guī)整的方式來使播放速度和時(shí)間恢 復(fù)到原來。 本程序使用的是用重疊疊加算法來達(dá)到 時(shí)長(zhǎng)規(guī)整。 重疊疊加

52、算法原理： 它分為兩個(gè)階段分解和合成 將原始信號(hào)以幀長(zhǎng) N，幀間距sa 進(jìn)行分解，然后以幀間距 ss 進(jìn)行合成。 sa 與ss 的的比值決定了時(shí)長(zhǎng)規(guī)整因子 F=sa/ss?！痉抡娉绦颉浚?）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