1、1,第一講 信號(hào)及其運(yùn)算的MATLAB表示,1.1連續(xù)信號(hào),1、指數(shù)信號(hào) 形式：y=A*exp（a*t）,2、正弦信號(hào) 形式：y=A*cos（w0*t+f） 或：y=A*sin（w0*t+f）,2,3、抽樣函數(shù) 用sinc函數(shù)表示，其定義為：sinc(t)=sin(t)/( t) MATLAB語言調(diào)用形式: y=sinc(t) t=-3*pi:pi/100:3*pi; y=sinc(t/pi); plot(t,y), grid on,3,4、矩形脈沖信號(hào) MATLAB語言調(diào)用形式: y=rectpuls(t,width) 產(chǎn)生一個(gè)幅值為1,寬度為width、相對(duì)于t=0點(diǎn)左右對(duì)稱的矩形波信號(hào)

2、t=0:0.001:4;T=1; y=rectpuls(t-2*T,2*T); plot(t,y), grid on axis(0,4,-0.5,1.5) 以t-2*T=0,即t=2*T為 對(duì)稱中心,4,周期性矩形脈沖信號(hào) MATLAB語言調(diào)用形式: y=square(t,duty) 產(chǎn)生一個(gè)周期為2 ,幅值為1的周期性方波, duty表示占空比(duty cycle) t=-0.0625:0.0001:0.0625; y=square(2*pi*30*t,75); plot(t,y),grid on axis(-0.1,0.1,-1.5,1.5) 占空比75%,5,5、三角脈沖信號(hào) MATL

3、AB語言調(diào)用形式: y=tripuls(t,width,skew) 產(chǎn)生一個(gè)最大幅值為1,寬度為width、斜度為skew的三角波信號(hào). t=-3:0.0001:3; y=tripuls(t,4,0.5); plot(t,y),grid on axis(-3,3,-1.5,1.5),6,周期性三角脈沖信號(hào) MATLAB語言調(diào)用形式: y=sawtooth(t,width) 產(chǎn)生一個(gè)周期為2 , 最大幅值為1,最小幅值為-1周期性三角波信號(hào). t=-5*pi:0.0001:5*pi; y=sawtooth(t,0.5); plot(t,y),grid on axis(-16,16,-1.5,1.

4、5),7,6、一般周期性脈沖信號(hào) MATLAB語言調(diào)用形式: y=pulstran(t,d,func) 依據(jù)名為func的連續(xù)函數(shù)并以之為一個(gè)周期,從而產(chǎn)生一串周期性的連續(xù)函數(shù)。用于指定周期性的 偏移量，這個(gè)func函數(shù)會(huì)被計(jì)算length(d)次,最后值相當(dāng)于:y=func(t-d(1)+ func(t-d(2)+.,更一般的調(diào)用形式為:y=pulstran(t,d,func,p1,p2) 。其中p1,p2為需要傳送給func函數(shù)的額外輸入?yún)?shù)值,最后值相當(dāng)于: y=func(t-d(1),p1,p2,.)+ func(t-d(2),p1,p2,.)+.,8,t=0:0.0001:1;d=0

5、:1/3:1; y=pulstran(t,d,rectpuls,0.1); figure(1),plot(t,y),grid on axis(0,1,-0.1,1.5) yy=pulstran(t,d,tripuls,0.1,-1);figure(2),plot(t,yy),grid on axis(0,1,-0.1,1.5),9,1.2 離散信號(hào),1、指數(shù)序列 a.k,2、正弦序列,離散的指數(shù)序列與正弦序列用MATLAB的表示與連續(xù)信號(hào)類似,只不過用stem函數(shù)而不是用plot函數(shù)來畫出序列的波形.,10,3、單位沖激序列,單位沖激序列,k=-50:50; y=zeros(1,50),1,z

6、eros(1,50); stem(k,y),function y,k=impseq(k0,k1,k2); k=k1:k2; y=k-k0=0,寫成MATLAB函數(shù)文件形式若 k-k0在k1k k2范圍內(nèi),11,4、單位階躍序列,單位階躍序列,k=-50:50; y=zeros(1,50), ones(1,51); stem(k,y), axis(-60,60,-0.5,1.5),寫成函數(shù)文件形式若 uk-k0 在k1k k2,function y,k=stepseq(k0,k1,k2); k=k1:k2; y=k-k0=0,12,1.3 信號(hào)運(yùn)算的MATLAB實(shí)現(xiàn),2、離散序列的差分與求和,1

7、、信號(hào)的尺度變換、反轉(zhuǎn)和平移,差分f(k)=f(k)-f(k-1)用diff函數(shù) 求和用sum函數(shù),13,3、連續(xù)信號(hào)的微分與積分,微分可用diff函數(shù) 積分用quad函數(shù), 調(diào)用形式: quad(function_name,a,b),例:對(duì)于三角波，畫出df(t)/dt和t- f(t)dt的波形,function y=funtrip(t) y=tripuls(t,4,0.5);,14,h=0.001,t=-3:h:3; y1=diff(funtrip(t)*1/h; figure(1),plot(t(1:length(t)-1),y1); title(df/dt);grid on,xlabe

8、l(t),15,for x=1:length(t) y2(x)=quad(funtrip,-3,t(x); end figure(2),plot(t,y2); title(Intergral of df);grid on,xlabel(t),16,第二講 時(shí)域分析的MATLAB表示,2.1連續(xù)時(shí)間系統(tǒng)的零狀態(tài)響應(yīng),LTI連續(xù)時(shí)間系統(tǒng)以常系數(shù)微分方程描述，系統(tǒng)的零狀態(tài)響應(yīng)可通過求解初始狀態(tài)為零的微分方程得到。 形式：y=lsim（sys,f,t） sys是LTI系統(tǒng)模型，f是系統(tǒng)輸入信號(hào)向量，是計(jì)算系統(tǒng)響應(yīng)的抽樣點(diǎn)向量。 LTI系統(tǒng)模型sys要通過MATLAB中的tf函數(shù)來獲得: sys=tf(

9、b,a), 其中b、a分別為微分方程右端和左端各項(xiàng)的系數(shù)向量 。,17,若系統(tǒng)的微分方程為：,系統(tǒng)的輸入信號(hào)為f(t)=10sin(t)，計(jì)算物體的位移y(t)。,ts=0;te=5;dt=0.01; sys=tf(1,1, 2, 100); t=ts:dt:te; f=10*sin(2*pi*t); y=lsim(sys,f,t); plot(t,y); xlabel(t(sec); ylabel(y(t);grid on,18,2.2連續(xù)時(shí)間系統(tǒng)的沖擊響應(yīng)和階躍響應(yīng),1、沖擊響應(yīng) 形式：y=impulse（sys,t） sys是LTI系統(tǒng)模型, sys=tf(b,a), 其中b、a分別為微

10、分方程右端和左端各項(xiàng)的系數(shù)向量 。,2、階躍響應(yīng) 形式：y=step(sys,t) sys是LTI系統(tǒng)模型, sys=tf(b,a), 其中b、a分別為微分方程右端和左端各項(xiàng)的系數(shù)向量 。,19,例：在一水平面上物體與彈簧組成的力學(xué)系統(tǒng)中，彈簧一端固定，物體的質(zhì)量m=1kg，彈簧的彈性系數(shù)ks=100N/m，物體與地面的摩擦系數(shù)fd=2Ns/m，在物體的水平方向施加一外力f(t)，外力是強(qiáng)度為10的沖擊信號(hào)，若系統(tǒng)的微分方程為：,求物體的位移y(t)。,20,ts=0;te=5;dt=0.01; sys=tf(10,1 2 100); t=ts:dt:te; y=impulse(sys,t);

11、 plot(t,y); xlabel(t(sec); ylabel(y(t); grid on,21,2.3 離散時(shí)間系統(tǒng)的零狀態(tài)響應(yīng),LTI離散時(shí)間系統(tǒng)常用線性常系數(shù)差分方程描述：,已知差分方程的n個(gè)初始狀態(tài)及輸入fk，可由下式迭代算出系統(tǒng)的輸出：,fk和yk 分別表示系統(tǒng)的輸入和輸出, n是差分方程的階數(shù)。,22,在零初始狀態(tài)下， MATLAB調(diào)用的形式為： y=filter（b，a，f）,23,R=51;d=rand(1,51)-0.5; k=0:R-1; s=2*k.*(0.9.k); f=s+d; subplot(1,2,1) plot(k,d,r-.,k,s,b:,k,f,k-);

12、 xlabel(Time k);legend(dk,sk,fk); M=5;b=ones(M,1)/M;a=1; y=filter(b,a,f); subplot(1,2,2) plot(k,s,b:,k,y,r-); xlabel(Time k);legend(sk,yk);,24,25,2.4 離散時(shí)間系統(tǒng)的沖擊響應(yīng)和階躍響應(yīng),1、單位沖擊響應(yīng) 形式：h=impz（b,a,k） 其中b、a分別為差分方程右端和左端的系數(shù)向量 ，k是輸出序列的取值范圍，h是系統(tǒng)的單位沖擊響應(yīng),2、單位階躍響應(yīng) 形式：h=stepz（b,a,k）,26,例：用impz函數(shù)求下列離散時(shí)間系統(tǒng)的單位沖擊響應(yīng)hk，并

13、與理論值hk=-(-1)k+2(-2)-k，k0進(jìn)行比較： yk+3yk-1+2yk-2=fk,k=0:10; a=1 3 2;b=1; h=impz(b,a,k); subplot(2,1,1) stem(k,h);title(單位沖擊相應(yīng)的近似值);grid on; hk=-(-1).k+2*(-2).k; subplot(2,1,2) stem(k,hk);title(單位沖擊相應(yīng)的理論值);grid on;,27,28,2.5 卷積的計(jì)算,兩個(gè)離散序列卷積和：c=conv(a,b),29,第三講 頻域分析的MATLAB表示,3.1周期信號(hào)的頻域分析,1、周期性離散時(shí)間信號(hào) 離散傅里葉級(jí)

14、數(shù)(DFS)和逆離散傅里葉級(jí)數(shù)（IDFS）表示為：,其中 WN=e -j2/N,30,MATLAB表示為：F=fft(f) f=ifft(F),例：利用MATLAB計(jì)算一周期矩形波序列的DFS系數(shù)。,N=32;M=4;%定義周期矩形波序列的參數(shù) %產(chǎn)生該周期矩形波序列 f=ones(1,M+1),zeros(1,N-2*M-1),ones(1,M); F=fft(f);%計(jì)算DFS系數(shù) m=0:N-1; figure(1),stem(m,real(F);title(Fm的實(shí)部);xlabel(m); figure(2),stem(m,imag(F);title(Fm的虛部);xlabel(m)

15、; fr=ifft(F);%計(jì)算由DFS系數(shù)Fm重建的fk figure(3)stem(m,real(fr);title(重建的fk);xlabel(k);,31,32,其中func是被積函數(shù)的函數(shù)文件名，a、b分別為定積分的下上限，Tol表示指定允許的相對(duì)或絕對(duì)積分誤差，非零的Trace表示以被積函數(shù)的點(diǎn)繪圖形式來跟蹤函數(shù)生成的返回值，若Tol和Trace均賦予“”，則兩者均自動(dòng)使用缺省值；”p1,p2,.”表示被積函數(shù)所需的多個(gè)額外輸入?yún)?shù),y=quad8(func,a,b) y=quad8(func,a,b,Tol,Trace,p1,p2,.),3.2 非周期信號(hào)的頻域分析,1、非周期性

16、連續(xù)時(shí)間信號(hào)的頻譜 quad8和 quadl可以計(jì)算非周期性連續(xù)時(shí)間信號(hào)的頻譜（高階算法）。,33,例：利用MATLAB采用數(shù)值方法近似計(jì)算 三角波信號(hào) 的 頻譜。,function y=sfl(t,w) %對(duì)于不同的，計(jì)算出傅里葉變換中被積函數(shù)的值 y=(abs(t)=1).*(1-abs(t).*exp(-j*w*t);,34,w=linspace(-6*pi,6*pi,512); N=length(w);F=zeros(1,N); for k=1:N F(k)=quad8(sfl,-1,1,w(k); end figure(1),plot(w,real(F); xlabel(omega)

17、;ylabel(F(jomega); title(三角波信號(hào)近似頻譜)； figure(2),plot(w,real(F)-sinc(w/2/pi).2); % sinc函數(shù)值即為理論計(jì)算結(jié)果 xlabel(omega);ylabel(error); title(計(jì)算誤差),35,36,2、非周期性離散時(shí)間信號(hào)的頻譜 當(dāng)序列的DTFT（離散時(shí)間傅里葉變換）可寫成ej的有理多項(xiàng)式時(shí)，MATLAB的函數(shù)freqz可用來計(jì)算其DTFT的值。設(shè)DTFT的有理多項(xiàng)式為：,函數(shù)freqz的一般調(diào)用形式：h=freqz(b,a,w),W為頻率抽樣點(diǎn)向量,37,例：利用MATLAB畫出a=0.9時(shí) 的幅度頻譜

18、，,b=1;a1=1,-0.9;a2=1,0.9; w=linspace(0,2*pi,512); h1=freqz(b,a1,w);h2=freqz(b,a2,w); plot(w/pi,abs(h1),w/pi,abs(h2),k:) xlabel(omega/pi); legend(alpha=0.9,alpha=-0.9);,38,39,1、連續(xù)時(shí)間系統(tǒng)的頻率特性,3.3 系統(tǒng)的頻率特性分析,當(dāng)頻率響應(yīng)H(jw)是jw的有理多項(xiàng)式時(shí)，即：,或?qū)憺镠(s)的形式：,40,MATLAB提供freqs函數(shù)計(jì)算連續(xù)時(shí)間系統(tǒng)的頻率響應(yīng): H=freqs(b,a,w) H,w=freqs(b,a,

19、N) 其中，W為頻率抽樣點(diǎn)向量, MATLAB自動(dòng)選擇N個(gè)頻率點(diǎn) 計(jì)算其頻率響應(yīng)， N的缺省值為200,41,w=linspace(0,6*pi,200);b=1;a=1,2,2,1; H=freqs(b,a,w); subplot(2,1,1)plot(w,abs(H) set(gca,xtick,0,1,2,3,4,5);set(gca,ytick,0,0.4,0.707,1); xlabel(omega(rad/s);ylabel(|H(jomega)|);grid on,42,subplot(2,1,2)plot(w,angle(H) set(gca,xtick,0,1,2,3,4,5

20、);set(gca,ytick,-4:2:4); xlabel(omega(rad/s);ylabel(phi(rad);grid on,43,2、離散時(shí)間系統(tǒng)的頻率特性,當(dāng)頻率響應(yīng)H(ej)是ej的有理多項(xiàng)式時(shí)，即：,或?qū)憺镠(z)的形式，即：,44,MATLAB提供freqz函數(shù)計(jì)算離散時(shí)間系統(tǒng)的頻率響應(yīng): H=freqz(b,a,w) H,w=freqz(b,a,N) 其中，W為頻率抽樣點(diǎn)向量, MATLAB在單位圓的上半周上等間距分布的N個(gè)頻率點(diǎn)計(jì)算其頻率響應(yīng)， N的缺省值為512,45,MATLAB音頻支持的函數(shù),1、sound 重放矢量成為聲音， 格式：sound(y,fs,bit

21、s) 2、soundsc 聲音自動(dòng)定標(biāo)和重放矢量， 格式：soundsc(y,fs,bits,slim), slim = slow, shigh 3、wavplay 使用視窗音頻的輸出裝置重放聲音 格式： wavplay(y,fs),46,4、wavread 讀*.wav聲音文件 x,fs,bits=wavread(file) 5、wavrecord 使用視窗音頻的輸入裝置記錄聲音 格式： wavrecord(n,fs,ch) 6、wavwrite 寫*.wav聲音文件 格式： wavwrite(y,fs,nbits,wavefile),47,程序設(shè)計(jì)實(shí)例,1 、 音頻信號(hào)的分析與合成,采用M

22、ATLAB分析WAV文件。,步驟一：選擇一個(gè)WAV文件作為分析的對(duì)象 。（*.wav）,步驟二：讀WAV文件數(shù)據(jù)并畫時(shí)域圖形 。,步驟三：進(jìn)行 FFT變換并畫頻域圖形 。,步驟四：進(jìn)行該聲波主要頻譜的分析。,步驟五：根據(jù)該聲音的頻譜，反演時(shí)域圖形 。（有失真）,常用函數(shù)：,x，fs，bits=wavread(filename),d= FFT (w，l),sound(w，fs，bits),步驟六：進(jìn)行付立葉逆變換IFFT并畫頻域圖形 。,48,音頻信號(hào)的分析與合成程序,返回,%= clear; if nargin1;action=initialized;end; %顯示打開文件的對(duì)話框 fnam

23、e,pname=uigetfile(SOUND1.wav,Open Wave File); file=pname,fname; x,fs,bits=wavread(file); % 讀入聲音文件（*.wav） sound(x,fs,bits); % 數(shù)據(jù)通過聲卡轉(zhuǎn)換為聲音 %=選擇聲道 %y=x(:); % 雙聲道信號(hào)波形數(shù)據(jù) y=x(:,1); % 左聲道信號(hào)波形數(shù)據(jù) %y=x(:,2); % 右聲道信號(hào)波形數(shù)據(jù),49,%=disp(按任意鍵顯示左聲道聲音波形！); disp(-); pause plot(y); set(gca,FontName,Arial,FontSize,10); title(左聲道聲音波形); xlabel(采樣點(diǎn)數(shù)); disp(按任意鍵聽左聲道信號(hào)三次！); disp(-); pause,sound(y,fs,bits); % 數(shù)據(jù)通過聲卡轉(zhuǎn)換為聲音 pause,sound(y,fs,bits); pause,sound(y,fs,bits); k=fft(y,32768); disp(按任意鍵顯示左聲道聲音頻域的幅值！); disp(-