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

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

2、變化，生成模擬其變化的指數(shù)信號。4） 分別錄制一段男聲、女聲信號，進(jìn)行音頻信號的讀取與播放，畫出其時域波形?！緶剀疤崾尽?1)利用MATLAB函數(shù) wavread(file)讀取.wav格式文件。(2)利用MATLAB函數(shù) 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）觀察一定時期內(nèi)的股票上證指數(shù)變化，生成模擬其變化的指數(shù)信號。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）分別錄制一段男聲、女聲信號，進(jìn)行音頻信號的讀取與播放，畫出其時域波形。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)【仿真結(jié)果】1）生成一個正弦信號 改變其角頻率和初始相位2）生成一個幅度為1、基頻為2Hz、占空比為50%的周期方波3）觀察一定時期內(nèi)的股票上證指數(shù)變化，生成模擬其變化的指數(shù)信號。4）分別錄制一段男聲、女聲信號，進(jìn)行音頻信號的讀取與播放，畫出其時域波形?！窘Y(jié)果分析】(1) 隨著正弦信號

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

6、將原始音頻信號在時域上進(jìn)行翻轉(zhuǎn)，【題目分析】【仿真程序】 原始音頻信號y,x,nbits=wavread('D:音樂動畫歌曲（207曲）灼眼的夏娜無損 17曲Light My Fire.wav');plot(y)1） 將原始音頻信號在時域上進(jìn)行延展、壓縮 將原始音頻信號在時域上延展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）將原始音頻信號在時域上進(jìn)行幅度放大與縮小 將原始音頻信號的幅度擴(kuò)大為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）將原始音頻信號在時域上進(jìn)行翻轉(zhuǎn)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)');【仿真結(jié)果】 原始音頻信號1）將原始音頻信號在時域上進(jìn)行延展、壓縮 將原始音頻信號在時域上延展3倍 將原始音頻信號在時域上壓縮為1/52）將原始音頻信號在時域上進(jìn)行幅度放大與縮小 將原始音頻信號的幅度擴(kuò)大為2倍 將原始音頻信號的幅度縮小為1/23）將原始音頻信號在時域上進(jìn)行翻轉(zhuǎn)【結(jié)果分析】以時間為橫軸的圖比以頻率為橫軸的圖看起來復(fù)雜得多，時域上的分析在現(xiàn)階段很困難。對于原始音頻信號，時域上的分析太復(fù)雜，而如果在頻域上分析，則簡化很多?！咀灾鲗W(xué)習(xí)內(nèi)容】信號的獲取以及如何使用；MATLAB的一些基本語法；如何根據(jù)MATLAB

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

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

13、1);原始信號為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時M=10時【結(jié)果分析】1. 首先，我們?nèi)×薓=5，去噪后的信號與原信號相比有一定的相近，即去噪功能相等。2. 其次，我們?nèi)×薓=200，此時的操作已經(jīng)不能稱之為去噪，信號有一定的失真。3. 最后，我們?nèi)×薓=1，去噪后的信號加上噪音的信號基本無區(qū)別?？傻媒Y(jié)論：M太大時，信號會失去失真度；M太小，平均值范圍太小與加噪信號區(qū)別不大。

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

17、.意思d和x的維數(shù)沒有保持一致。改用randn函數(shù)時，令d=0.01*randn(size(x)*0.3-0.5就成功地把噪聲加入了進(jìn)去?！締栴}探究】 在對兩個函數(shù)進(jìn)行加法計(jì)算時要保持維數(shù)相同?！狙杏憙?nèi)容】題目2：連續(xù)信號卷積的近似計(jì)算兩個連續(xù)信號的卷積定義為為了進(jìn)行數(shù)值計(jì)算，需對連續(xù)信號進(jìn)行抽樣。記xk=x(kD), hk=h(kD), D為進(jìn)行數(shù)值計(jì)算的抽樣間隔。則連續(xù)信號卷積可近似的寫為(1)這就可以利用conv函數(shù)可近似計(jì)算連續(xù)信號的卷積。設(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)不是時限信號，則用上面的方法進(jìn)行近似計(jì)算會遇到什么問題？給出一種解決問題的方案；(e) 若將x(t)和h(t)近似表示為推導(dǎo)近似計(jì)算卷積的算法。取相同的抽樣間隔，比較兩種方法的計(jì)算卷積誤差。【題目分析】（a），求出y(t)的函數(shù)表達(dá)式，并畫出它的圖像，便于和后面通過算出的信號作比較。（b），Matlab中不可以直接進(jìn)行連續(xù)的信號的卷積，必須得先對連續(xù)信號采樣，得到離散的信號，然后求得兩個離散信號的卷積，得到另外一個離散信號yk，最后再把離散的信號連續(xù)化，得到要得到的y(k)。通過改變的

19、值來與（a）中所得結(jié)果進(jìn)行對比。（c），由于積分可以用求和的極限來表示，所以，當(dāng)足夠小時，成立。（d），Matlab不能表示出非時限信號。我們可以用符號運(yùn)算或者（e），兩個相同的寬度的矩形波的卷積為一個三角波。根據(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ì)算卷積誤差。【結(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越小的時候，圖像越平滑，越接近于利用解析法求得的結(jié)果。dt越大是，誤差越大。(d)若x(t)和h(t)不是時限信號：方法一： matlab軟件不能取到無窮大，函數(shù)值在到達(dá)某一值時會變成0。解決辦法是我們可以根據(jù)自己的需要，設(shè)定x(t)和h(t)在某一范圍內(nèi)的函數(shù)值，其他值均為0。方法二：我們

25、可以采用matlab中的符號運(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）的方法比較好，更接近于理論值?！咀灾鲗W(xué)習(xí)內(nèi)容】卷積函數(shù)conv的調(diào)用格式【閱

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

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

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

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

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

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)從以下幾方面對結(jié)果進(jìn)行分析：(1) 你所選擇樂器演奏的音階，其時域波形的包絡(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)換；基波分析、諧波分析。【閱讀文獻(xiàn)】1陳后金 胡健 薛健 信號與系統(tǒng). M 北京：高等教育出版社2肖燕彩 邱

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

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

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抽樣定理，將正弦信號以大于其頻率二倍的抽樣頻率抽樣時，信號可以被較好還原，而當(dāng)抽樣頻率小于正弦頻率二倍抽樣時，信號將發(fā)生失真?！痉抡娉绦颉?產(chǎn)生指定頻率的正弦聲音信號并按固定的抽樣頻率播放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;%計(jì)算幅度頻譜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【仿真結(jié)果】我們選取2000Hz和7200Hz的兩個正弦信號以抽樣頻率fs=8000Hz抽樣后的時域波形和頻譜進(jìn)行對比，仿真后畫出的圖形如下所示：【結(jié)果分析】如我們所料，7200Hz的正弦信號在以抽樣頻率fs=8000Hz抽樣后，發(fā)生了嚴(yán)重的失真，由于其周期性而使得頻率變?yōu)?00Hz。這說明，抽樣頻率低于正弦信號最大頻率的兩倍抽樣時，信號會發(fā)生失真。這種聲音信號實(shí)際頻率隨正弦信號頻率與抽樣頻率關(guān)系而變化的現(xiàn)象引起了我們的興趣，于是我們進(jìn)行了以下仿真來探

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

46、帶通信號時常遵循帶通信號采樣定理，即采樣頻率fs>2(fH-fL)即可。（其中，fH為信號的最大頻率，fL為信號的最大頻率）經(jīng)討論，我們決定用matlab生成一個有一定頻譜特點(diǎn)的帶限聲音信號，然后對這個信號進(jìn)行抽樣。通過觀察抽樣前后的頻譜變化，來驗(yàn)證以上分析的正確性?！痉抡娉绦颉縡s=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);【仿真結(jié)果】【結(jié)果分析】 我們生成的信號的頻率從1950Hz到2250Hz之間變化，按我們

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

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

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

52、算法原理： 它分為兩個階段分解和合成 將原始信號以幀長 N，幀間距sa 進(jìn)行分解，然后以幀間距 ss 進(jìn)行合成。 sa 與ss 的的比值決定了時長規(guī)整因子 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