数字信号处理实验程序2.docx

数字信号处理实验程序2.docx

《数字信号处理实验程序2.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验程序2.docx（17页珍藏版）》请在冰豆网上搜索。

数字信号处理实验程序2

2.1

clc

closeall;

n=0:

15;

p=8;q=2;

x=exp（-（n-p.^2/q;

figure（1;

subplot（3,1,1;

stem（n,x;

title（'exp（-（n-p^2/q,p=8,q=2';

xk1=fft（x,16;

q=4;

x=exp（-（n-p.^2/q;

subplot（3,1,2;

xk2=fft（x,16;

stem（n,x;

title（'exp（-（n-p^2/q,p=8,q=4';

q=8;

x=exp（-（n-p.^2/q;

xk3=fft（x,16;

subplot（3,1,3;

stem（n,x;

title（'exp（-（n-p^2/q,p=8,q=8';%时域特性figure（2;

subplot（3,1,1;

stem（n,abs（xk1;

title（'exp（-（n-p^2/q,p=8,q=2';

subplot（3,1,2;

stem（n,abs（xk2;

title（'exp（-（n-p^2/q,p=8,q=4';

subplot（3,1,3;

stem（n,abs（xk3;

title（'exp（-（n-p^2/q,p=8,q=8';%频域特性%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p=8;q=8;

figure（3;

subplot（3,1,1;

stem（n,x;

title（'exp（-（n-p^2/q,p=8,q=8';

xk1=fft（x,16;

p=13;

x=exp（-（n-p.^2/q;

subplot（3,1,2;xk2=fft（x,16;

stem（n,x;

title（'exp（-（n-p^2/q,p=13,q=8';

p=14;

x=exp（-（n-p.^2/q;

xk3=fft（x,16;

subplot（3,1,3;

stem（n,x;

title（'exp（-（n-p^2/q,p=14,q=8';%时域特性figure（4;

subplot（3,1,1;

stem（n,abs（xk1;

title（'exp（-（n-p^2/q,p=8,q=8';

subplot（3,1,2;

stem（n,abs（xk2;

title（'exp（-（n-p^2/q,p=13,q=8';

subplot（3,1,3;

stem（n,abs（xk3;

title（'exp（-（n-p^2/q,p=14,q=8';%频域特性2.2

clc

closeall

n=0:

15;

f=0.0625;a=0.1;

xb1=exp（-a*n.*sin（2*pi*f*n;

xk1=fft（xb1,16;

f=0.4375;

xb2=exp（-a*n.*sin（2*pi*f*n;

xk2=fft（xb2,16;

f=0.5625;

xb3=exp（-a*n.*sin（2*pi*f*n;

xk3=fft（xb3,16;

figure（1;

subplot（3,1,1;

stem（n,xb1;

title（'f=0.0625';

subplot（3,1,2;

stem（n,xb2;

title（'f=0.4375';

subplot（3,1,3;

stem（n,xb3;

title（'f=0.5625';

figure（2;

subplot（3,1,1;

stem（n,abs（xk1;

title（'f=0.0625';

subplot（3,1,2;

stem（n,abs（xk2;

title（'f=0.4375';

subplot（3,1,3;

stem（n,abs（xk3;

title（'f=0.5625';

2.3

clc

closeall;

N=8;

n=0:

N-1;

xc=0:

7;

xd=0:

7;

form=0:

7;

if（0<=m&m<=3

xc（m+1=m;

xd（m+1=4-m;

elseif（4<=m&m<=7

xc（m+1=8-m;

xd（m+1=m-4;

end

end

figure（1;

subplot（2,1,1;

stem（n,xc;

title（'三角波序列';

subplot（2,1,2;

stem（n,xd;

title（'反三角波序列';

xck=fft（xc,N;

xdk=fft（xd,N;

figure（2;

subplot（2,1,1;

stem（n,abs（xck;

title（'三角波序列8点fft';subplot（2,1,2;

stem（n,abs（xdk;

title（'反三角波序列8点fft';figure（3;

xck1=fft（xc,32;

xdk1=fft（xd,32;

subplot（2,1,1;stem（0:

4*N-1,abs（xck1;

title（'三角波序列32点fft';

subplot（2,1,2;

stem（0:

4*N-1,abs（xdk1;

title（'反三角波序列32点fft';

2.4

clc

closeall

N=16;detaf=1/16;

n=0:

N-1;

x=sin（2*pi*0.125*n+cos（2*pi*（0.125+detaf*n;

xk=fft（x,N;

figure（1;

subplot（2,1,1;

stem（n,abs（xk;

title（'N=16,Δf=1/16';

detaf=1/64;

x=sin（2*pi*0.125*n+cos（2*pi*（0.125+detaf*n;

xk1=fft（x,N;

subplot（2,1,2;

stem（n,abs（xk1;

title（'N=16,Δf=1/64';

N=128;detaf=1/16;

n=0:

N-1;

x=sin（2*pi*0.125*n+cos（2*pi*（0.125+detaf*n;

xk=fft（x,N;

figure（2;

subplot（2,1,1;

stem（n,abs（xk;

title（'N=128,Δf=1/16';

detaf=1/64;

x=sin（2*pi*0.125*n+cos（2*pi*（0.125+detaf*n;

xk1=fft（x,N;

subplot（2,1,2;

stem（n,abs（xk1;

title（'N=128,Δf=1/64';

2.5

clc

closeall

N=16;

p=8;q=2;

a=.1;

f=.0625;

n=0:

N-1;

xa=exp（-（n-p.^2/q;

xb=exp（-a*n.*sin（2*pi*f*n;

xak=fft（xa,N;

xbk=fft（xb,N;

f=ifft（xak.*xbk,N;

figure（1;

subplot（4,1,1;

stem（0:

N-1,xa;

title（'xa（n=exp（-（n-p^2/q';subplot（4,1,2;

stem（0:

N-1,xb;

title（'xb（n=exp（-a*n*sin（2*pi*f*n';subplot（4,1,3;

stem（0:

N-1,f;

title（'xa（n⊙xb（n';

xak1=fft（xa,2*N;

xbk1=fft（xb,2*N;

h=ifft（xak1.*xbk1,2*N;

subplot（4,1,4;

stem（0:

2*N-2,h（1:

2*N-1;

title（'xa*xb';

2.6

clc

closeall

N=512;

xe=rand（1,N-0.5;

form=0:

7;

if（0<=m&m<=3

xc（m+1=m;

elseif（4<=m&m<=7

xc（m+1=8-m;

end

end

xck=fft（xc,71;

f=0:

518;

fm=0:

141;

form=1:

8;

xm=xe（（m-1*64+1:

m*64;

xmk=fft（xm,71;

fm（1:

71=fm（72:

142;

fm（72:

142=ifft（xmk.*xck,71;

ifm==2

f（1:

7=fm（1:

7;

end

ifm>1

f（（m-2*64+8:

（m-1*64=fm（8:

64;

f（（m-1*64+1:

（m-1*64+7=fm（65:

71+fm（72:

78;

end

ifm==8

f（（m-1*64+8:

m*64+7=fm（79:

142;

end

end

stem（0:

512+6,f;%%%%重叠相加法

title（'重叠相加法';

form=1:

8

xm1（8:

71=xe（（m-1*64+1:

m*64;

ifm==1

xm1（1:

7=0;

elsexm1（1:

7=xe（（m-1*64-6:

（m-1*64;

end

xmk1=fft（xm1,71;

fm1=ifft（xmk1.*xck,71;

f1（（m-1*64+1:

m*64=fm1（8:

71;

end

xmk2=fft（xe（506:

512,71;

h=ifft（xmk2.*xck,71;

f1（513:

519=h（8:

14;

figure（2

stem（0:

518,f1;

title（'重叠保留法';

fa=conv（xc,xe;

detaf=f1-f;

figure（3;

stem（0:

518,detaf;

2.7

clc

closeall

N=16;

p=8;

q=2;

a=.1;

n=0:

N-1;

f=0.0625;

xa=exp（-（n-p.^2/q;

xb=exp（-a*n.*sin（2*pi*f*n;

xak=fft（xa,N*2;

xbk=fft（xb,2*N;

rm=real（ifft（conj（xak.*xbk;

rm=[rm（N+2:

2*Nrm（1:

N];

m=（-N+1:

N-1;

stem（m,rm;

title（'线性相关';

figure（2;

xak1=fft（xa,N;

xbk1=fft（xb,N;

rm1=real（ifft（xak1.*xbk1,N;

stem（n,rm1;

3.1

clear

closeall

wc=2*1000*tan（2*pi*300/（2*1000;wt=2*1000*tan（2*pi*200/（2*1000;[N,wn]=cheb1ord（wc,wt,0.8,20,'s';[B,A]=cheby1（N,0.8,wn,'high','s';[num,den]=bilinear（B,A,1000;[h,w]=freqz（num,den;

f=w/pi*500;

plot（f,20*log10（abs（h;

axis（[0,500,-80,10];

grid;

xlabel（'频率/Hz';

ylabel（'幅度/dB';

3.2

clear;

closeall;

wc=2*pi*200;

wr=2*pi*300;

[n1,wn1]=buttord（wc,wr,1,25,'s';[B,A]=butter（n1,wn1,'low','s';[num1,den1]=impinvar（B,A,1000;[h1,w]=freqz（num1,den1;

wc=2*1000*tan（2*pi*200/（2*1000;wr=2*1000*tan（2*pi*300/（2*1000;[n,wn]=buttord（wc,wr,1,25,'s';[B,A]=butter（n,wn,'low','s';

[num2,den2]=bilinear（B,A,1000;[h2,w]=freqz（num2,den2;

f=w/pi*500;

plot（f,20*log10（abs（h1,'-.',f,20*log10（abs（h2,'-';

legend（'脉冲响应不变法','双线性变换法';axis（[0,500,-80,10]

grid;

xlabel（'频率/Hz';

ylabel（'幅度/dB';

3.3

clear

closeall

wc=2*8000*tan（2*pi*1200/（2*8000;

wr=2*8000*tan（2*pi*2000/（2*8000;

[n,wn]=buttord（wc,wr,0.5,40,'s';

[B,A]=butter（n,wn,'low','s';

[num,den]=bilinear（B,A,8000;

[h1,w]=freqz（num,den;

[n,wn]=cheb1ord（wc,wr,0.5,40,'s';

[B,A]=cheby1（n,0.5,wn,'low','s';

[num,den]=bilinear（B,A,8000;

[h2,w]=freqz（num,den;

[n,wn]=ellipord（wc,wr,0.5,40,'s';

[B,A]=ellip（n,0.5,40,wn,'low','s';

[num,den]=bilinear（B,A,8000;

[h3,w]=freqz（num,den;

f=w/pi*4000;

plot（f,20*log10（abs（h1,'-.',f,20*log10（abs（h2,'.-',f,20*log10（abs（h3,'-';

axis（[0,4000,-100,10];

grid;

xlabel（'频率/Hz';

ylabel（'幅度/dB';

legend（'Butterworth','Chebyshev','Ellicptic'holdon

x=1200;

plot（x

3.4

clear

closeall;

[n,wc]=buttord（[20003000],[15006000],3,20,'s';

[B,A]=butter（n,wc,'bandpass','s';

[num,den]=bilinear（B,A,30000;

[h1,w]=freqz（num,den;

[num,den]=impinvar（B,A,30000;

[h2,w]=freqz（num,den;

f=w/pi*15000;

plot（f,20*log10（abs（h1,'-.',f,20*log10（abs（abs（h2;

gridon;

axis（[0,3000,-100,0];

legend（'脉冲响应不变法','双线性变换';xlabel（'频率/Hz';

ylabel（'幅度/dB';

3.5

clear

closeall

w1=2*10000*tan（2*pi*1000/（2*10000;

w2=2*10000*tan（2*pi*2000/（2*10000;

w3=2*10000*tan（2*pi*500/（2*10000;

w4=2*10000*tan（2*pi*3000/（2*10000;[n,wn]=cheb1ord（[w3w4],[w1w2],3,18,'s'[B,A]=cheby1（n,3,wn,'stop','s';

[num,den]=bilinear（B,A,10000;

[h,w]=freqz（num,den;

f=w/pi*5000;

plot（f,20*log10（abs（h;

axis（[0,5000,-100,0];

grid

xlabel（'频率/Hz';

ylabel（'幅度/dB';

4.1

%N=45,计算并画出矩形框、汉明窗、布莱克曼窗的归一化的幅度谱,并比较各自的主要特点

%（1矩形窗（RectangleWindow调用格式:

w=boxcar（n,根据长度n产生一个矩形窗w。

%（2三角窗（TriangularWindow调用格式:

w=triang（n,根据长度n产生一个三角窗w。

%（3汉宁窗（HanningWindow调用格式:

w=hanning（n,根据长度n产生一个汉宁窗w。

%（4海明窗（HammingWindow调用格式:

w=hamming（n,根据长度n产生一个海明窗w。

%（5布拉克曼窗（BlackmanWindow调用格式:

w=blackman（n,根据长度n产生一个布拉克曼窗w。

%（6恺撒窗（KaiserWindow调用格式:

w=kaiser（n,beta,根据长度n和影响窗函数旁瓣的β参数产生一个恺撒窗w

clearall

closeall

N=45;

w1=boxcar（N;

w2=hamming（N;

w3=blackman（N;

[h,w]=freqz（w1,N;

figure（1

plot（w/pi,20*log10（abs（h;

axis（[0,1,-80,10];

gridon

xlabel（'归一化频率/π';

ylabel（'幅度/dB';

title（'矩形窗';

figure（2

[h,w]=freqz（w2,N;

plot（w/pi,20*log10（abs（h;

axis（[0,1,-80,10];

gridon

xlabel（'归一化频率/π';

ylabel（'幅度/dB';

title（'汉明窗';

figure（3

[h,w]=freqz（w3,N;

plot（w/pi,20*log10（abs（h;

axis（[0,1,-150,10];

gridon

xlabel（'归一化频率/π';

ylabel（'幅度/dB';

title（'布莱克曼窗';

4.2

%N=15,带通滤波器的两个通带边界分别是

w1=0.3π,w2=0.5π。

用汉宁窗设计此线性相位滤波器，观察%它的实际3dB和20dB带宽。

N=45，重复这一设计，观察幅频和相位特性的变化，注意N变化的影响。

closeallclearallN=15;w1=0.3;w2=0.5;w=hanning（N;n=0:

N-1;alfa=（N-1/2;h=fir1（N-1,[w1w2],w;[h1,w3]=freqz（h,1;figure（1subplot（2,1,1;plot（w3/pi,20*log10（abs（h1;gridon;axis（[0,1,-80,10];xlabel（'归一化频率/π';ylabel（'幅度/dB';subplot（2,1,2;plot（w3/pi,angle（h1;gridon;axis（[0,1,-4,4];xlabel（'归一化频率/π';ylabel（'角度/rad';N=45;w=hanning（N;n=0:

N-1;alfa=（N-1/2;h=fir1（N-1,[w1w2],w;[h1,w3]=freqz（h,1;figure（2subplot（2,1,1;plot（w3/pi,20*log10（abs（h1;gridon;axis（[0,1,-80,10];xlabel（'归一化频率/π';ylabel（'幅度/dB';subplot（2,1,2;plot（w3/pi,angle（h1;gridon;axis（[0,1,-4,4];xlabel（'归一化频率/π';ylabel（'角度/rad';4.3closeallclearallN=15;w1=0.3;w2=0.5;wn1=boxcar（N;wn2=blackman（N;hn1=fir1（N-1,[w1w2],wn1;hn2=fir1（N-1,[w1w2],wn2;[h1,w3]=freqz（hn1,1;figure（1plot（w3/pi,20*log10（abs（h1;gridon;axis（[0,1,-80,10];xlabel（'归一化频率/π';ylabel（'幅度/dB';title（'矩形窗,N=15';[h1,w3]=freqz（hn2,1;figure（2plot（w3/pi,20*log10（abs（h1;gridon;axis（[0,1,-80,10];xlabel（'归一化频率/π';ylabel（'幅度/dB';title（'布莱克曼窗,N=15';N=45;wn1=boxcar（N;wn2=blackman（N;hn1=fir1（N-1,[w1w2],wn1;hn2=fir1（N-1,[w1w2],wn2;[h1,w3]=freqz（hn1,1;figure（3plot（w3/pi,20*log10（abs（h1;gridon;axis（[0,1,-80,10];xlabel（'归一化频率/π';ylabel（'幅度/dB';title（'矩形窗,N=45';[h1,w3]=freqz（hn2,1;figure（4

plot（w3/pi,20*log10（abs（h1;gridon;axis（[0,1,-110,10];xlabel（'归一化频率/π';ylabel（'幅度/dB';title（'布莱克曼窗,N=45';4.4closeallclearallN=40;%beta=4forn=1:

3ifn==1beta=4;elseifn==2beta=6;elsebeta=10;end;w=kaiser（N,beta;h=fir1（N-1,[0.20.40.60.8],w;[h1,w1]=freqz（h,1;figure（nsubplot（2,1,1;plot（w1/pi,20*log10（abs（h1;gridon;axis（[0,1,-80,10];xlabel（'归一化频率/π';ylabel（'幅度/dB';ifn==1title（'beta=4';elseifn==2title（'beta=6';elsetitle（'beta=10';end;subplot（2,1,2;plot（w1/pi,angle（h1;gridon;axis（[0,1,-4,4];xlabel（'归一化频率/π';ylabel（'角度/rad';end4.5clearallcloseallN=45;k=0:

N-1;fork=0:

N-1w=2*pi/N*k;hk（1,k+1=0;if（（w>=0.2*pi（w<=0.4*pi||（w>=0.6*piw<=0.8*pi||（w>=1.2*piw<=1.4*pi||（w>=1.6*pi&&w<=1.8*pihk（1,k+1=1;endendk=0:

N-1;hk（1,5=0.5;hk（1,11=0.5;hk（1,14=0.5;hk（1,20=0.5;hk（1,27=0.5;hk（1,33=0.5;hk（1,36=0.5;hk（1,42=0.5;thetak=-k*2*pi/N*（（N-1/2;hk1=hk.*exp（j*thetak;hn=ifft（hk1;[h1,w1]=freqz（hn,1;plot（w1/pi,20*log10（abs（h1;gridon;axis（[