《数字信号处理》上机全部源代码调试通过完整版.docx
《《数字信号处理》上机全部源代码调试通过完整版.docx》由会员分享，可在线阅读，更多相关《《数字信号处理》上机全部源代码调试通过完整版.docx（26页珍藏版）》请在冰点文库上搜索。
《数字信号处理》上机全部源代码调试通过完整版
《数字信号处理》上机全部源代码调试通过，完整版
（高西全，第四版）
实验一
%实验1：
系统响应及系统稳定性
cｌｏse all;ｃlear all
%调用fliter解差分方程，由系统对ｕｎ的响应判断稳定性
％内容1：
调用filter解差分方程， 由系统对u（n）的响应判断稳定性
A=[1,-0．9];B=[0.05，0.0５];
ｘ1n=［1 111１111zeroｓ（1，5０）］；
x2n=ｏnes（1,１28）；
ｈｎ=imｐz（B,Ａ，５8）；
subplｏt（2,2，１）；ｙ='h（n）';tsｔｅm（ｈn,y）;
ｔitle（'（a） 系统单位脉冲响应h（n）'）
y1n=filter（B,A,x1n）;
sｕbpｌｏｔ（２,2，２）;ｙ＝'y1（n）';tstem（y1n,y）;
ｔitle（'（b）系统对R8（ｎ）的响应y1（n）'）
ｙ2n=ｆｉｌter（B,A,x2ｎ）；
subplｏt（2,2,４）;y=＇y２（n）';tｓtem（y2n， y）;
title（'（c）系统对ｕ（n）的响应y2（ｎ）＇）
y1n=fiｌｔｅr（B,A,x１n）;
ｓｕbpｌot（2，2,2）;y=＇y１（n）';tstem（y1n，y）;
tiｔlｅ（'（b）系统对R8（n）的响应y1（n）'）
y２n=fiｌter（B,A,ｘ2n）；
suｂｐlｏt（２,2,4）;y＝＇y２（n）';ｔstｅm（y2ｎ,y）;
title（'（c）系统对u（n）的响应ｙ2（ｎ）'）
%内容2:
调用ｃonv函数计算卷积
x１ｎ=[11１ 1 1111];％产生信号x1ｎ＝R8n
h1n=[ones（1,10） zeroｓ（1,１0）］;
h２n=[１２.52．５1 ｚｅrｏs（1,１0）]
y２１n=ｃonｖ（h1n，x1n）;
y22n=conv（h２n，x１ｎ）；
fｉgure
（2）
subｐlot（2,2,1）;y='h1（n）';tstｅm（h1n,y）;
％调用函数tsｔeｍ绘图
titｌe（'（ｄ）系统单位脉冲响应h1（n）＇）
subｐlot（2,2,2）;y=＇y２1（ｎ）';tｓｔeｍ（y２1ｎ,y）;
titｌe（'（e） ｈ1（n）与R８（ｎ）的卷积y21（n）'）
subｐｌot（2,２，3）;y＝'h2（n）＇;tstem（ｈ2n,y）;%调用函数ｔsｔeｍ绘图
tｉtｌe（＇（ｆ）系统单位脉冲响应h2（n）'）
suｂploｔ（2， 2,4）;y＝＇y22（n）';ｔｓｔeｍ（y22n,ｙ）;
title（'（g） ｈ2（n）与R８（n）的卷积ｙ22（n）'）
％=====＝=＝=＝=＝＝==＝==＝===＝===＝＝========＝
%内容3：
谐振器分析
uｎ=ones（1,256）; %产生信号ｕn
ｎ=0：
255;
ｘsｉn＝sｉn（0.0１4*n）+siｎ（0.4*n）； %产生正弦信号
A=[１,-1．８23７,0.９801]；
Ｂ=[1/100．４９,0,-1/100.４9];
%系统差分方程系数向量B和A
ｙ31n=filtｅr（B，A,ｕn）; %谐振器对uｎ的响应y31n
ｙ３2n=filｔｅr（Ｂ,Ａ，xｓiｎ）;
%谐振器对正弦信号的响应y3２ｎ
figurｅ（3）
subploｔ（２,1,1）;y='y3１（ｎ）'；tstem（y31n,y）
titlｅ（＇（h）谐振器对ｕ（n）的响应ｙ３1（n）'）
suｂｐlｏt（2,1,２）;y='y3２（n）＇;ｔsｔｅｍ（ｙ３２ｎ,y）;
ｔｉtlｅ（＇（i）谐振器对正弦信号的响应y32（n）'）
functiontstem（xn,ｙｎ）
ｎ=0：
lengｔｈ（xn）－１;
stｅｍ（n,xn,＇.'）；
xlabeｌ（'n'）;ylａｂｅl（'yn'）;
%xlabel（'n'）;ylabeｌ（yn）;
ａｘis（[0,n（eｎd）,min（xｎ）,1.2*ｍａx（xn）］）；
实验二
％时域采样理论验证程序exp2a.ｍ
Tp＝64／1０00;％观察时间Tp=64微秒
%产生Ｍ长采样序列x（ｎ）
%Fs=１０00;T=1／Fs;
Fs=1000;Ｔ=1／Ｆs;
M=Tp*Fs；n=0：
M-1;
A=444.1２8; alph=pi＊５0*2＾0.5； omega=ｐi*５０*2^0．5;
xnｔ=A＊ｅxp（-alph*ｎ*T）．*ｓin（omeｇａ＊n*T）;
Ｘk=T*fft（xnt,M）;％M点FFT［xnt）］
yn＝'xa（ｎT）';sｕbplｏt（３,2,１）;
tsｔem（xnt,ｙn）;%调用自编绘图函数tsｔem绘制序列图
boｘon；titlｅ（'（a）Fｓ=10０0Hz'）；
k＝0:
M-１; fk=k/Tp;
subpｌoｔ（3,2,２）； plｏt（fk,abs（Xk））；
tiｔlｅ（'（a） T＊FT[xa（nT）]，Fs=1000Hｚ'）;
ｘlabel（'f（Ｈz）'）;yｌabeｌ（'幅度'）;
ａxis（［0,Fs,0，1.２*maｘ（abs（Xｋ））]）
％＝=========＝===========＝===＝======
％Fs=300Ｈz和Fs=200Hz的程序与上面Ｆs=10０0Ｈz的程序完全相同。
%％％%%%%%%%%％ｆs＝30０Hz
Ｔp＝64/1000;%观察时间Ｔp=64微秒
%产生M长采样序列x（n）
％Ｆs=1００0;T=1/Fs;
Ｆs=30０；T＝1/Fs;
M=Tp*Fs;n=0:
M-1;
A=４４４.128;ａlph=pi*５0*2^0．5;ｏmega=pi*50＊2^0.5；
xnｔ=Ａ*ｅxｐ（－alph*ｎ*T）．＊sin（omega*n*T）;
Xk=T*fft（ｘnt,M）;％M点FFT[xｎt）］
yｎ='ｘa（nT）＇;fｉgｕre（２）;sｕbpｌot（３,2,1）;
tstem（xｎt,yn）; %调用自编绘图函数tsｔeｍ绘制序列图
ｂoxon;tｉtｌe（＇（b）Fｓ=300Hｚ'）;
k=0:
Ｍ-1； fk＝k／Tp;
subｐloｔ（3,2,2）;pｌot（fk,abs（Xｋ））;
tｉtｌe（'（b）T*FT[ｘa（ｎＴ）],Fs=1000Hｚ'）；
xlabｅｌ（'f（Hｚ）'）;ｙｌａbｅｌ（'幅度'）;
axis（［0,Fs,0,1.2*ｍax（aｂs（Xｋ））]）
%%%%%%％%%%%%％%%％%%％fs=200Hz
Tp=64/１000;%观察时间Ｔp=６4微秒
%产生M长采样序列ｘ（n）
%Ｆs=100０；Ｔ=1/Fs;
Ｆs=２00; Ｔ=1／Fs；
M＝Tｐ*Fｓ；ｎ=０:
M－1；
A=44４．１２8;aｌｐh=pi＊５0*2^0.5;omegａ=pi*50＊2＾0．5；
xnｔ＝A*exp（-alph*n＊T）.*sｉn（ｏｍｅｇａ*n*Ｔ）;
Xk=T＊fft（xnt,Ｍ）;%M点FＦT[xnｔ）］
ｙn='xａ（nＴ）'；ｆiｇure（3）;sｕｂplot（3,2,1）;
tstem（ｘnｔ，yｎ）;%调用自编绘图函数ｔｓｔem绘制序列图
boxｏn;titｌe（'（c）Fｓ=200Hz'）;
ｋ＝0:
Ｍ-１；ｆk=k/Tp;
sｕｂpｌot（3,2,2）；plot（fｋ,ａbs（Xk））;
title（＇（ｃ）T＊FT[xa（ｎＴ）]，Fs＝100０Hz'）;
xlaｂeｌ（'f（Ｈｚ）'）;ylabel（＇幅度'）;
aｘis（[0,Fs,0,1．2＊max（abs（Ｘk））]）
％频域采样理论验证程序exp2ｂ.m
M=２7;N=32;ｎ=0:
M;
％产生M长三角波序列x（n）
xa=0:
floｏｒ（M/2）; xｂ=ｃｅil（M/2）-１:
-１:
０;ｘn=[ｘa,ｘb］;
Xk＝fｆt（xｎ，1024）;%1024点ＦFT[x（n）],用于近似序列x（n）的TF
X３2k＝ffｔ（ｘn,32）;％32点FFT[x（n）]
ｘ3２n=ｉｆft（X32k）;%32点ＩＦFT[X３2（k）］得到x3２（n）
X16k=X32k（1:
２：
N）;%隔点抽取X３2k得到Ｘ1６（K）
ｘ16n=ifft（X16k,Ｎ/２）;％1６点IＦＦT[X16（ｋ）]得到x16（n）
ｓubｐlｏt（3，2,2）;ｓteｍ（n,xn,'.'）;bｏｘon
ｔitle（'（b） 三角波序列x（n）＇）;xlaｂel（＇n'）;yｌaｂel（＇ｘ（n）'）；ａxiｓ（[０,３２,０,２0]）
k=0:
1023;wk=2*k／1024；ﻩ%
subpｌｏt（3，２,1）;ｐlot（wｋ,abs（Xk））；title（'（a）ＦT［x（n）]＇）;
ｘlabｅｌ（'\oｍｅgａ／＼pi＇）;ylａbｅl（'|X（e＾j＾\omegａ）|'）;axｉｓ（[０,1,0，20０］）
k＝0:
Ｎ/2-1；
ｓubpｌot（3,2,3）;steｍ（k,ａbs（X16ｋ），'．＇）;boxoｎ
titlｅ（'（ｃ）16点频域采样'）;xlaｂｅｌ（'k＇）;ylaｂｅl（＇|X_1_6（k）|＇）;ａxis（［0,８，0,200］）
n１=0:
N/2-1;
sｕbｐlot（3,2,4）;stｅm（ｎ1,ｘ16n,＇.＇）;boｘon
tiｔle（'（ｄ） 16点IDＦT［X_１_6（k）]＇）;xｌａｂel（'n'）;ylabel（'x_1_６（n）'）；axｉｓ（［0,32,0,２0]）
k＝0:
N-1;
sｕbｐｌｏｔ（3,2,５）；sｔem（k,ａbs（Ｘ32k）,'.'）;boｘ on
tiｔｌe（'（e）32点频域采样'）;xlabｅl（'k'）；ｙｌaｂel（'|Ｘ_３_２（k）|＇）;axis（[０,１6,０,２0０］）
n1=0:
Ｎ-1;
suｂｐloｔ（3,2,６）;stｅm（n1，x32n,'.'）;boxon
title（'（f） 32点IDＦT［Ｘ_３_2（k）]'）;xlabeｌ（＇ｎ'）;yｌabel（'x_3_２（ｎ）＇）;axｉs（[０,３2,0,2０]）
ｆunctioｎ tｓtem（xn,ｙn）
ｎ=0:
lengtｈ（xn）－1;
stem（n,xｎ，'.＇）;
ｘlａｂel（'n'）;yｌabｅl（'ｙn'）;
%xlaｂｅｌ（'n' ）;ylaｂel（yｎ）;
axis（[0,ｎ（eｎd）,min（xn）,１.2＊ｍａx（xn）]）;
实验三
%实验三程序exp3.ｍ
%用FFT对信号作频谱分析
％clｅar aｌｌ;cｌｏseall
％实验内容（１）===＝===＝======＝＝=======＝=======＝
x1ｎ=[ones（1,4）];％产生序列向量x１（n）=R４（n）
M=8;xa＝1:
（Ｍ/2）; xb=（Ｍ/2）：
－1:
１;
x2n＝[xa，xb]; %产生长度为8的三角波序列x2（ｎ）
x３ｎ=[xb,xa］；
Ｘ1k8=fft（x1n,８）; %计算x1n的8点DFT
X1k１6＝ｆft（ｘ１n,16）; ％计算ｘ１n的１６点DＦT
Ｘ2ｋ8=ｆft（x2n,8）; ％计算x１ｎ的8点DFT
X2ｋ16＝fｆt（x2n，16）; ％计算x1n的１6点DFＴ
Ｘ3k8＝fft（x3n，８）； %计算x1n的８点DFT
Ｘ３k１6=fｆt（ｘ3n，1６）; ％计算x1n的１６点DFT
％以下绘制幅频特性曲线
subpｌot（2,2，１）;mｓtem（X１k8）;
%绘制8点DFT的幅频特性图
ｔｉtle（＇（1a）８点DFT［x_1（n）]'）;xlabel（'ω/π＇）;
ylabel（'幅度'）;
aｘis（[０,２,0,1．2*max（ａｂs（X１ｋ８））]）
subpｌot（2,２,3）；mstem（Ｘ1k16）;
％绘制１6点ＤFT的幅频特性图
tｉtle（'（1b）1６点DFT[x＿1（n）]＇）;xlaｂel（'ω/π'）；
ylabel（'幅度'）;
axｉs（[０,2,0,１.２*mａｘ（abｓ（X1k1６））]）
figure
（2）
ｓｕｂｐｌoｔ（2，2，1）;msteｍ（X2k８）;
%绘制8点ＤFT的幅频特性图
titlｅ（＇（2a）８点DＦT［x_2（ｎ）]＇）；xlａbeｌ（'ω/π'）;
ｙｌabｅｌ（'幅度＇）；
aｘis（[０,2,０,1.２*max（aｂｓ（X2ｋ8））]）
suｂploｔ（２，2,２）;ｍstem（X2k1６）;
％绘制16点DFT的幅频特性图
titｌe（'（２b）16点DFT[x_2（n）]'）;xlabel（＇ω/π'）;
ylabｅl（＇幅度'）;
aｘis（[0,2，0,1．2*max（aｂs（X2k１6））］）
ｓubpｌot（２，2，3）;mstｅｍ（X3k8）;
%绘制8点DFT的幅频特性图
tｉtle（'（３ａ）8点ＤFＴ[ｘ_3（ｎ）]＇）;ｘlabeｌ（'ω／π'）;
ylabel（'幅度'）;
axis（[０，２，0,１．2*ｍａx（aｂs（X3k8））］）
subplot（2，２，４）;mstem（X3k1６）;
%绘制16点DＦT的幅频特性图
title（'（3ｂ）1６点DＦT[x_３（n）]'）;xlaｂｅl（'ω/π'）;
ylａbeｌ（'幅度＇）;
ａxiｓ（[０,２,0,1．2*ｍａｘ（abs（Ｘ3ｋ1６））]）
％实验内容
（2）周期序列谱分析=============＝=====
N=8;n=0:
N-1; %FＦＴ的变换区间N＝8
x４ｎ=cos（ｐｉ＊n/4）;
x５ｎ＝ｃos（ｐi*ｎ／4）+cos（pi＊ｎ／8）;
X4k8=fft（x4n）; %计算x4n的8点DＦT
X5ｋ８＝fft（ｘ５ｎ）; %计算x5n的8点ＤFT
N=１６;n＝0:
N-1;%ＦFT的变换区间Ｎ＝１6
x4ｎ=ｃｏs（ｐｉ*n/4）;
x5n=cos（ｐi*n/4）+cos（pi*ｎ/8）;
X４k16=ｆｆt（x４ｎ）； %计算x４ｎ的1６点ＤＦＴ
X5ｋ１6=fft（x5n）; %计算x5ｎ的1６点DFT
fｉｇuｒe（3）
subｐloｔ（2，2,1）;msｔem（X4k８）;
%绘制8点DＦT的幅频特性图
ｔitle（＇（4a）8点DFＴ［x_４（n）]'）；
xlaｂｅl（＇ω／π'）;ylabel（'幅度＇）;
axis（[0，2，0，1.２＊maｘ（abs（X4k8））]）
subplot（2,2,3）；ｍstｅm（X4k1６）;
％绘制１6点ＤFT的幅频特性图
tiｔｌe（'（4b）16点DＦT[x＿4（ｎ）］'）;
ｘlａbeｌ（＇ω/π'）;ylaｂｅl（'幅度'）;
aｘｉs（[0,2，0,１．2*maｘ（abs（X4ｋ１6））]）
ｓｕｂplot（2,２,2）;mstem（X5k８）；
%绘制8点DＦＴ的幅频特性图
title（'（5ａ） 8点ＤFＴ[x_5（n）］'）;ｘlabel（＇ω/π＇）;
yｌabel（＇幅度'）;
axiｓ（［0,2,0,1.2*max（abs（Ｘ5k8））]）
subpｌot（2，2,4）;mｓtem（Ｘ5ｋ16）；
%绘制16点DＦT的幅频特性图
tｉtle（'（５b）16点DFT［ｘ＿５（ｎ）]'）;xlabeｌ（'ω／π'）;
ylaｂel（'幅度'）；
aｘis（[０,2,0,1.2＊ｍax（ａbｓ（X５ｋ１６））]）
%实验内容（3）模拟周期信号谱分析===＝==＝===＝==＝===
figurｅ（4）
Fs＝６4;T=1／Fs；
N=16;n=0:
N-1; %FFT的变换区间N=16
x６nＴ＝ｃoｓ（8*pi*n*Ｔ）+cos（１6*pi*ｎ*Ｔ）+cos（20＊pｉ*n＊Ｔ）;
%对x6（t）16点采样
X6ｋ１6＝fft（x６nT）； %计算x６nT的１6点DFT
Ｘ6ｋ16=fftｓｈift（X6k16）; %将零频率移到频谱中心
Tｐ＝N*T；F=１/Tp； ％频率分辨率F
k=-Ｎ/2：
N/2-1；fk=k*Ｆ;
％产生1６点DFＴ对应的采样点频率（以零频率为中心）
sｕbｐlot（3,１,1）；steｍ（fk,aｂs（X6k１6）,'.＇）;
boxon %绘制8点DFT的幅频特性图
tｉtｌｅ（＇（６a）16点|DＦT[x_6（nT）]|'）;xlabel（'f（Hz）＇）;
％%%%%%%%%
％%%%%％%％％
ylabｅｌ（'幅度'）;
aｘis（［－N*Ｆ/2-１,N*F／2－１,0,1.2＊max（abｓ（X６k1６））］）
N=3２;ｎ=0:
Ｎ-１; %FFT的变换区间N=1６
x６nT=ｃos（8＊pi＊ｎ*Ｔ）＋cｏs（16＊pi＊n*Ｔ）+ｃos（20＊pi*ｎ*T）；
％对x6（t）３2点采样
Ｘ6ｋ3２＝fft（x６nT）； %计算x６nT的3２点DFT
X6ｋ32＝fftshifｔ（Ｘ6ｋ３2）; %将零频率移到频谱中心
Tp=N*Ｔ;Ｆ=１/Ｔp;%频率分辨率F
k＝－N／2:
N/2－1;
fk＝k*F;
％产生16点ＤFＴ对应的采样点频率（以零频率为中心）
subplot（3,1,2）;sｔeｍ（ｆk，abs（X6k３2）, '.＇）;
ｂox oｎ %绘制8点DFＴ的幅频特性图
ｔitle（'（6b）32点｜DFT［x_６（ｎＴ）］|＇）；
xlabel（'f（Ｈz）'）;ｙｌabeｌ（'幅度'）;
axｉs（[－N*F/2－1,Ｎ*Ｆ/２-1，０,1.2*maｘ（abｓ（Ｘ6k３2））]）
N=64;n=0：
Ｎ-１; ％ＦFＴ的变换区间N=16
x6nT＝cos（8*ｐｉ*n*T）+cos（16*pi*n*Ｔ）+ｃos（2０*pi*n*Ｔ）;
％对ｘ６（t）６4点采样
X6k64=fft（x6nT）; %计算x6nT的６4点DFT
X6k６4=fftｓhift（X6k6４）; %将零频率移到频谱中心
Ｔｐ=N*T；Ｆ＝1／Tp; %频率分辨率Ｆ
k=-N／2：
N／2-1;fｋ=k＊F;
%产生１6点DFT对应的采样点频率（以零频率为中心）
suｂｐlｏｔ（3,1,3）；stｅm（ｆｋ,abｓ（X6ｋ64）,'.'）;
boｘon％绘制８点DFT的幅频特性图
title（'（６a）64点|DFT[x_6（ｎＴ）]｜'）;
ｘlabel（＇f（Hz）'）;ylaｂｅl（'幅度'）;
aｘiｓ（[-N*F/2-１,N*F/2-1，０,１.2*max（ａbs（X6k6４））］）
fｕnctｉoｎmｓtem（Xk）
%ｍsｔem（Xk）绘制频域采样序列向量Xk的幅频特性图
Ｍ＝lｅｎgtｈ（Xk）;
k=0:
Ｍ-1;wk=2＊k/M;%产生M点DFT对应的采样点频率（关于ｐi归一化值）
stem（wｋ,abs（Xk）,'.'）；ｂoxon;%绘制Ｍ点ＤFT的幅频特性图?
xｌａｂel（'w／\ｐi＇）；ylabel（'幅度'）;
ａxiｓ（[0,２，0,1.2＊mａx（abｓ（Xk））］）;
ﻬ
实验四
%实验四程序eｘp4．ｍ
%IＩＲ数字滤波器设计及软件实现
%ｃlearall;closeall
Fｓ=10000;Ｔ=1/Fs； %采样频率
%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st
ｓt＝msｔg;
%低通滤波器设计与实现=＝====＝＝＝====＝＝======＝=＝＝
fｐ=2８0;fs=45０;
wp=2*fｐ/Fs; ws＝2*fs/Fs;rp=０.1;ｒs=6０;
％DF指标（低通滤波器的通、阻带边界频率）
[N,wp］=ellｉpord（wp，ws,ｒp,rs）；
%调用ｅlｌipord计算椭圆DF阶数Ｎ和通带截止频率wp
［B,A]＝ｅlｌｉp（Ｎ，rp,ｒs,wp）;
%调用ellｉp计算椭圆带通DＦ系统函数系数向量B和Ａ
y1t=fiｌter（B,Ａ,st）; ％滤波器软件实现
% 低通滤波器设计与实现绘图部分
fｉgure（２）;subｐlot（3，１,1）；
myｐloｔ（Ｂ,A）;
％调用绘图函数mｙｐlｏt绘制损耗函数曲线
ｙt='y_１（t）';
sｕbpｌｏt（3,1，2）;ｔｐlｏt（y１t,T,yt）;
%调用绘图函数tploｔ绘制滤波器输出波形
%======＝======＝=====＝=====＝＝===＝＝========＝=＝
％带通滤波器设计与实现=＝===＝=========＝==＝＝＝====
ｆｐl=4４0；fｐu=５60;fsl=２75;fｓu=9０0;
wｐ=[2*fpl/Fs,2*fｐu/Fs]；
ｗｓ=[２*fｓl／Fs,２*ｆsu/Fs]; ｒｐ=0.1；rs＝60；
[Ｎ，wp]＝ellipord（ｗp，ws,rp，rs）;
％调用elliporｄ计算椭圆DＦ阶数N和通带截止频率wｐ
[Ｂ,Ａ]=elliｐ（N,rp,rs,wｐ）；
％调用ｅｌlｉｐ计算椭圆带通DF系统函数系数向量B和A
y2t=filtｅｒ（B,A,st）; %滤波器软件实现
% 带通滤波器设计与实现绘图部分（省略）
％补充１
figure（２）； sｕbplot（3，1，3）;
ｍｙpｌot（Ｂ,A）；
%--
ｆｉgｕrｅ（3）;
yt='y＿2（t）＇；
subｐlot（3,１，1）；tplot（y2t，T，yt）;
%========＝=＝==＝========＝=========＝＝=＝=======
%高通滤波器设计与实现===============＝=＝=========
fp=890;ｆｓ=６0０;
ｗp=2＊fｐ/Fｓ; ｗs=２*fs/Fｓ; rp=０.1；rs＝６０；
％DF指标（低通滤波器的通、阻带边界频率）
［Ｎ,wｐ］=ｅｌliｐord（wp,ws,rp,rs）;
%调用elliｐord计算椭圆ＤF阶数N和通带%截止频率wp
［B,A］=ｅlｌip（N,ｒp,rｓ,wp,'ｈｉgh'）;
%调用eｌｌiｐ计算椭圆带通DF系统函数系数向量Ｂ和A
ｙ3t=fｉｌｔｅr（B,Ａ,sｔ）;%滤波器软件实现
%高低通滤波器设计与实现绘图部分（省略）
%补充2
fｉｇｕｒe（3）;subploｔ（３，１,２）;
myｐｌot（Ｂ,A）;
％--
fiｇuｒe（3）;
yt='ｙ＿3（t）';
suｂplｏt（3，１,3）;tploｔ（ｙ3t,T,ｙt）;
%＝
升级会员 