1、% 實驗 1：系統(tǒng)響應及系統(tǒng)穩(wěn)定性close all;clear all%= 內容 1：調用 filter 解差分方程，由系統(tǒng)對 u(n) 的響應判斷穩(wěn)定性 = A=1,-0.9;B=0.05,0.05; %系統(tǒng)差分方程系數向量 B 和 Ax1n=1 1 1 1 1 1 1 1 zeros(1,50); % 產生信號 x1(n)=R8(n) x2n=ones(1,128); % 產生信號 x2(n)=u(n) hn=impz(B,A,58); % 求系統(tǒng)單位脈沖響應 h(n) subplot(2,2,1);y='h(n)'tstem(hn,y); % 調用函數 tstem 繪圖

2、 title('(a) 系統(tǒng)單位脈沖響應 h(n)');box on y1n=filter(B,A,x1n); % 求系統(tǒng)對 x1(n) 的響應 y1(n) subplot(2,2,2);y='y1(n)'tstem(y1n,y);title('(b) 系統(tǒng)對 R8(n)的響應 y1(n)');box ony2n=filter(B,A,x2n);% 求系統(tǒng)對 x2(n)的響應 y2(n)subplot(2,2,4);y='y2(n)'tstem(y2n,y);title('(c)系統(tǒng)對 u(n)的響應 y2(n)'

3、);box on %= 內容 2：調用 conv 函數計算卷積 = x1n=1 1 1 1 1 1 1 1 ;% 產生信號 x1(n)=R8(n)h1n=ones(1,10) zeros(1,10);h2n=1 2.5 2.5 1 zeros(1,10); y21n=conv(h1n,x1n); y22n=conv(h2n,x1n);figure(2) subplot(2,2,1);y='h1(n)'tstem(h1n,y);%調用函數 tstem 繪圖title('(d) 系統(tǒng)單位脈沖響應 h1(n)');box on subplot(2,2,2);y=

4、9;y21(n)'tstem(y21n,y);title('(e) h1(n)與 R8(n)的卷積 y21(n)');box on subplot(2,2,3);y='h2(n)'tstem(h2n,y);% 調用函數 tstem 繪圖title('(f) 系統(tǒng)單位脈沖響應 h2(n)');box on subplot(2,2,4);y='y22(n)'tstem(y22n,y);title('(g) h2(n)與 R8(n)的卷積 y22(n)');box on%= 內容 3：諧振器分析 = un=one

5、s(1,256);% 產生信號 u(n)n=0:255;xsin=sin(0.014*n)+sin(0.4*n); % 產生正弦信號A=1,-1.8237,0.9801;B=1/100.49,0,-1/100.49;% 系統(tǒng)差分方程系數向量B 和 Ay31 n=filter(B,A,un);%諧振器對 u(n)的響應 y31(n)y32n=filter(B,A,xsin);% 諧振器對 u(n)的響應 y31(n)figure(3) subplot(2,1,1);y='y31(n)'tstem(y31n,y);title('(h)諧振器對 u(n)的響應 y31(n)&

6、#39;);box on subplot(2,1,2);y='y32(n)'tstem(y32n,y);title('(i) 諧振器對正弦信號的響應 y32(n)');box on10.2.2 實驗程序清單1 時域采樣理論的驗證程序清單Tp=64/1000;%觀察時間 Tp=64 微秒%產生 M 長采樣序列 x(n)% Fs=1000;T=1/Fs;Fs=1000;T=1/Fs;M=Tp*Fs;n=0:M-1;A=444.128;alph=pi*50*2A0.5;omega=pi*50*2A0.5; xnt=A*exp(-alph*n*T).*sin(omega

7、*n*T);Xk=T*fft(xnt,M); %M 點 FFTxnt) yn='xa(nT)'subplot(3,2,1);tstem(xnt,yn);%調用自編繪圖函數 tstem 繪制序列圖box on;title('(a) Fs=1000Hz');k=0:M-1;fk=k/Tp; subplot(3,2,2);plot(fk,abs(Xk);title('(a) T*FTxa(nT),Fs=1000Hz');xlabel('f(Hz)');ylabel('幅度');axis(0,Fs,0,1.2*max(ab

8、s(Xk) %= % Fs=300Hz和Fs=200Hz的程序與上面 Fs=1000Hz完全相同。2 頻域采樣理論的驗證程序清單M=27;N=32;n=0:M;%產生 M 長三角波序列 x(n)xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=xa,xb;Xk=fft(xn,1024);%1024 點 FFTx(n),用于近似序列 x(n)的 TFX32k=fft(xn,32);%32 點 FFTx(n)x32n=ifft(X32k); %32 點 IFFTX32(k) 得到 x32(n) X16k=X32k(1:2:N);%隔點抽取 X32k 得到 X16(

9、K)x16n=ifft(X16k,N/2);%16 點 IFFTX16(k) 得到 x16(n)subplot(3,2,2);stem(n,xn,'.');box on title('(b) 三角波序列 x(n)');xlabel('n');ylabel('x(n)');axis(0,32,0,20) k=0:1023;wk=2*k/1024;%subplot(3,2,1);plot(wk,abs(Xk);title('(a)FTx(n)');xlabel('omega/pi');ylabel(&#

10、39;|X(eAjAomega)|');axis(0,1,0,200) k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');box ontitle('(c) 16 點頻域采樣 ');xlabel('k');ylabel('|X_1_6(k)|');axis(0,8,0,200) n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,'.');box ontitle('(d) 16 點 IDFTX_1_6(k)');xlabel

11、('n');ylabel('x_1_6(n)');axis(0,32,0,20) k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');box ontitle('(e) 32 點頻域采樣 ');xlabel('k');ylabel('|X_3_2(k)|');axis(0,16,0,200) n1=0:N-1;subplot(3,2,6);stem(n1,x32n,'.');box ontitle('(f) 32 點 IDFTX_3_2(k

12、)');xlabel('n');ylabel('x_3_2(n)');axis(0,32,0,20)10.3.2 實驗程序清單% 用 FFT 對信號作頻譜分析clear all;close all%實驗內容 (1)=x1n=ones(1,4);%產生序列向量 x1(n)=R4(n)M=8;xa=1:(M/2); xb=(M/2):-1:1; x2n=xa,xb; % 產生長度為 8 的三角波序列 x2(n)% 計算 x1n 的 8 點 DFT %計算 x1n 的 16 點 DFT% 計算 x1n 的 8 點 DFT%計算 x1n 的 16 點 DFT%

13、計算 x1n 的 8 點 DFT%計算 x1n 的 16 點 DFTx3n=xb,xa;X1k8=fft(x1n,8);X1k16=fft(x1n,16);X2k8=fft(x2n,8);X2k16=fft(x2n,16);X3k8=fft(x3n,8);X3k16=fft(x3n,16);%以下繪制幅頻特性曲線subplot(2,2,1);mstem(X1k8); % 繪制 8 點 DFT 的幅頻特性圖title('(1a) 8 點 DFTx_1(n)');xlabel(' 3 /n ');ylabel('幅度')；axis(0,2,0,1.2

14、*max(abs(X1k8)subplot(2,2,3);mstem(X1k16); % 繪制 16點 DFT 的幅頻特性圖 title('(1b)16 點 DFTx_1(n)');xlabel(' 3 / n ');ylabel('幅度');axis(0,2,0,1.2*max(abs(X1k16)figure(2) subplot(2,2,1);mstem(X2k8); % 繪制 8 點 DFT 的幅頻特性圖 title('(2a) 8 點 DFTx_2(n)');xlabel(' 3 /n ');ylabel

15、('幅度'); axis(0,2,0,1.2*max(abs(X2k8) subplot(2,2,2);mstem(X2k16); % 繪制 16 點 DFT 的幅頻特性圖 title('(2b)16 點 DFTx_2(n)');xlabel(' 3 / n ');ylabel('幅度'); axis(0,2,0,1.2*max(abs(X2k16) subplot(2,2,3);mstem(X3k8); % 繪制 8 點 DFT 的幅頻特性圖 title('(3a) 8 點 DFTx_3(n)');xlabel(

16、' 3 /n ');ylabel('幅度'); axis(0,2,0,1.2*max(abs(X3k8) subplot(2,2,4);mstem(X3k16); % 繪制 16 點 DFT 的幅頻特性圖 title('(3b)16 點 DFTx_3(n)');xlabel(' 3 / n ');ylabel('幅度');axis(0,2,0,1.2*max(abs(X3k16)%實驗內容 (2) 周期序列譜分析 =N=8;n=0:N-1;%FFT 的變換區(qū)間 N=8x4n=cos(pi*n/4);x5n=cos(

17、pi*n/4)+cos(pi*n/8);X4k8=fft(x4n);% 計算 x4n 的 8 點 DFTX5k8=fft(x5n);% 計算 x5n 的 8 點 DFTN=16;n=0:N-1; %FFT 的變換區(qū)間 N=16x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k16=fft(x4n);% 計算 x4n 的 16 點 DFTX5k16=fft(x5n);% 計算 x5n 的 16 點 DFTfigure(3)subplot(2,2,1);mstem(X4k8); % 繪制 8 點 DFT 的幅頻特性圖title('(4a) 8 點

18、 DFTx_4(n)');xlabel(' 3 /n ');ylabel('幅度')；axis(0,2,0,1.2*max(abs(X4k8)subplot(2,2,3);mstem(X4k16); % 繪制 16 點 DFT 的幅頻特性圖title('(4b)16 點 DFTx_4(n)');xlabel(' 3 / n ');ylabel('幅度');axis(0,2,0,1.2*max(abs(X4k16)subplot(2,2,2);mstem(X5k8); % 繪制 8 點 DFT 的幅頻特性圖t

19、itle('(5a) 8 點 DFTx_5(n)');xlabel(' 3 /n ');ylabel('幅度'); axis(0,2,0,1.2*max(abs(X5k8)subplot(2,2,4);mstem(X5k16); % 繪制 16 點 DFT 的幅頻特性圖 title('(5b)16 點 DFTx_5(n)');xlabel(' 3 / n ');ylabel('幅度'); axis(0,2,0,1.2*max(abs(X5k16)%實驗內容 (3) 模擬周期信號譜分析 =figure

20、(4)Fs=64;T=1/Fs;N=16;n=0:N-1;%FFT 的變換區(qū)間 N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%對 x6(t)16 點采樣X6k16=fft(x6nT);%計算 x6nT 的 16點 DFTX6k16=fftshift(X6k16);%將零頻率移到頻譜中心Tp=N*T;F=1/Tp;%頻率分辨率 Fk=-N/2:N/2-1;fk=k*F;%產生16點DFT對應的采樣點頻率以零頻率為中心subplot(3,1,1);stem(fk,abs(X6k16),'.');box on % 繪制 8

21、 點 DFT 的幅頻特性圖 title('(6a) 16 點 |DFTx_6(nT)|');xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16)N=32;n=0:N-1;%FFT 的變換區(qū)間 N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%對 x6(t)32 點采樣X6k32=fft(x6nT);%計算 x6nT 的 32 點 DFTX6k32=fftshift(X6k32);%將零頻率移到頻譜中心

22、Tp=N*T;F=1/Tp;%頻率分辨率 Fk=-N/2:N/2-1;fk=k*F;%產生16點DFT對應的采樣點頻率以零頻率為中心subplot(3,1,2);stem(fk,abs(X6k32),'.');box on % 繪制 8 點 DFT 的幅頻特性圖 title('(6b) 32 點 |DFTx_6(nT)|');xlabel('f(Hz)');ylabel('幅度');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32)N=64;n=0:N-1; %FFT 的變換區(qū)間 N=16 x6n

23、T=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); % 對 x6(t)64 點采樣 X6k64=fft(x6nT);%計算 x6nT 的 64 點 DFTX6k64=fftshift(X6k64);%將零頻率移到頻譜中心Tp=N*T;F=1/Tp;%頻率分辨率 Fk=-N/2:N/2-1;fk=k*F;%產生16點DFT對應的采樣點頻率以零頻率為中心subplot(3,1,3);stem(fk,abs(X6k64),'.'); box on% 繪制 8點 DFT 的幅頻特性圖 title('(6a) 64 點 |DFTx_6(n

24、T)|');xlabel('f(Hz)');ylabel(' 幅度 ');axis(-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64)% IIR 數字濾波器設計及軟件實現(xiàn)clear all;close allFs=10000;T=1/Fs; % 采樣頻率%調用信號產生函數 mstg 產生由三路抑制載波調幅信號相加構成的復合信號 st st=mstg;%低通濾波器設計與實現(xiàn) = fp=280;fs=450;wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF 指標低通濾波器的通、阻帶邊界頻 N,wp=ell

25、ipord(wp,ws,rp,rs); %調用ellipord計算橢圓 DF階數N和通帶截止頻率 wpB,A=ellip(N,rp,rs,wp); %調用ellip 計算橢圓帶通DF系統(tǒng)函數系數向量 B和Ay1t=filter(B,A,st); %濾波器軟件實現(xiàn)% 低通濾波器設計與實現(xiàn)繪圖局部figure(2);subplot(3,1,1);myplot(B,A); % 調用繪圖函數 myplot 繪制損耗函數曲線yt='y_1(t)'subplot(3,1,2);tplot(y1t,T,yt); %調用繪圖函數 tplot 繪制濾波器輸出波形%帶通濾波器設計與實現(xiàn) = fpl

26、=440;fpu=560;fsl=275;fsu=900;wp=2*fpl/Fs,2*fpu/Fs;ws=2*fsl/Fs,2*fsu/Fs;rp=0.1;rs=60;N,wp=ellipord(wp,ws,rp,rs); %調用ellipord計算橢圓DF階數N和通帶截止頻率 wpB,A=ellip(N,rp,rs,wp); % 調用ellip 計算橢圓帶通DF系統(tǒng)函數系數向量 B和A y2t=filter(B,A,st); %濾波器軟件實現(xiàn)% 帶通濾波器設計與實現(xiàn)繪圖局部省略%高通濾波器設計與實現(xiàn) = fp=890;fs=600;wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;r

27、s=60; %DF 指標低通濾波器的通、阻帶邊界頻 N,wp=ellipord(wp,ws,rp,rs); %調用ellipord計算橢圓DF階數N和通帶截止頻率 wpB,A=ellip(N,rp,rs,wp,'high'); %調用ellip計算橢圓帶通DF系統(tǒng)函數系數向量B和Ay3t=filter(B,A,st); %濾波器軟件實現(xiàn)% 上下通濾波器設計與實現(xiàn)繪圖局部省略2、實驗程序清單% FIR 數字濾波器設計及軟件實現(xiàn)clear all;close all;%=調用xtg產生信號xt, xt長度N=1000,并顯示xt及其頻譜，= N=1000;xt=xtg(N);fp=

28、120; fs=150;Rp=0.2;As=60;Fs=1000; % 輸入給定指標% (1) 用窗函數法設計濾波器wc=(fp+fs)/Fs;%理想低通濾波器截止頻率 (關于 pi 歸一化B=2*pi*(fs-fp)/Fs;%過渡帶寬度指標Nb=ceil(11*pi/B); %blackman 窗的長度 N hn=fir1(Nb-1,wc,blackman(Nb);Hw=abs(fft(hn,1024);% 求設計的濾波器頻率特性ywt=fftfilt(hn,xt,N);%調用函數 fftfilt 對 xt 濾波%以下為用窗函數法設計法的繪圖局部濾波器損耗函數，濾波器輸出信號波形%省略% (2) 用等波紋最正確逼近法設計濾波器fb=fp,fs;m=1,0;% 確定 remezord 函數所需參數 f,m,devdev=(10A(Rp/20)-1(10A(Rp/20)+1),10A(-As/20);Ne,fo,mo,W=remezord(fb,m,dev,Fs); % 確定 remez 函數所需參數 hn=remez(Ne,fo,mo,W);% 調用 remez 函數進行設計Hw=abs(fft(hn,1024);% 求設計的濾波器頻率特性yet=fftfilt(hn,xt,N);% 調用函數 fftfilt 對 x