数字信号处理MATLAB实验50题修改

1-1

clear ;close all;clc;

b=[1,1];a=[1,-0.5];

subplot(3,1,1);zplane(b,a);title('因果系统零极图'); n=0:50;

x=3*cos(pi*n/3);

y=filter(b,a,x);

subplot(3,1,2);stem(n,x,'.');title('输入x的波形'); subplot(3,1,3);stem(n,y,'.');title('输出y的波形');

1-2

clear ;close all;clc;

b=[1,1,1];a=[1,0.5,-0.25];

subplot(3,1,1);zplane(b,a);title('因果系统零极图'); n=0:50;

x=3*cos(pi*n/3);

y=filter(b,a,x);

subplot(3,1,2);stem(n,x,'.');title('输入x的波形'); subplot(3,1,3);stem(n,y,'.');title('输出y的波形');

2

clear ;close all;clc;

b=[0,1];a=[1,-1,-1];

x=impseq(0,-5,50);n=-5:50;

h=filter(b,a,x);

stem(n,h,'.');

title('单位脉冲响应')

sum(abs(h))

3

clear ;close all;clc;

b=[2];

a=[1 -0.8 -0.5];

subplot(4,1,1);

zplane(b,a);

title('系统的零极图');

[H,w]=freqz(b,a,100,'whole');

magH=abs(H);

phaH=angle(H);

subplot(4,1,2);

plot(w/pi,magH);

title('系统的幅频响应');

subplot(4,1,3);

plot(w/pi,phaH/pi);

title('系统的相频响应');

n=0:100;

x=impseq(0,0,100);

h=filter(b,a,x);

subplot(4,1,4);

stem(n,h,'.');

title('系统的冲激响应');

4

clear ;close all;clc;

b=[1 1];

a=[1 -0.9 0.81];

[H,w]=freqz(b,a,400,'whole'); magH=abs(H);

phaH=angle(H);

subplot(4,1,1);

plot(w/pi,magH);

title('系统的幅频响应'); subplot(4,1,2);

plot(w/pi,phaH/pi);

title('系统的相频响应');

n=0:200;

x=sin(pi*n/3)+5*cos(pi*n); y=filter(b,a,x);

subplot(4,1,3);

plot(n,x);

title('输入信号X');

subplot(4,1,4);

plot(n,y);

title('输出信号Y');

grid;

5

clear ;close all;clc;

x11=[1 1 1 1];

n=0:5;

x12=cos(pi*n/4);

y11=circonvt(x11,x12,8)

y12=conv(x11,x12)

y13=[y11(1:1:8),zeros(1,1)] e1=y13-y12

x21=[1 -1 1 -1];

x22=[1 0 -1 0];

y21=circonvt(x21,x22,5)

y22=conv(x21,x22)

y23=[y21(1:1:5),zeros(1,2)]

e2=y23-y22

n=0:15;

x31=cos(2*pi*n/32);

x32=sin(2*pi*n/32);

y31=circonvt(x31,x32,32)

y32=conv(x31,x32)

y33=[y32(1:1:31),zeros(1,1)]

e3=y31-y33

n=0:9;

x41=(0.8).^n;

x42=(-0.8).^n;

y41=circonvt(x41,x42,15)

y42=conv(x41,x42)

y43=[y41(1:1:15),zeros(1,4)]

e4=y43-y42

6

clear ;close all;clc;

x1=[2 1 1 2];

x2=[1 -1 -1 1];

n=[0:8-1];

y11=circonvt(x1,x2,4)

y12=circonvt(x1,x2,7)

y13=circonvt(x1,x2,8)

y2=conv(x1,x2)

%N最小值7

7--1

clear ;close all;clc;

x=[2,2,2,2,2,2,2,2];w=[0:1:500]*2*pi/500;

[H]=freqz(x,1,w);

magH=abs(H);phaH=angle(H);

subplot(2,2,1);plot(w/pi,magH);grid

xlabel('');ylabel('|x|');title('DTFT的幅度')

subplot(2,2,2);plot(w/pi,phaH/pi*180);grid

xlabel('以pi为单位的频率');ylabel('度');title('DTFT的相角') N=8;w1=2*pi/N;k=0:N-1;

X=dft(x,N);

magX=abs(X),phaX=angle(X)*180/pi

subplot(2,2,3);plot(w*N/(2*pi),magH,'--');

axis([-0.1,8.1,0,20]);

hold on

stem(k,magX);

ylabel('|x(k)|');title('DFT的幅度:N=8');text(4.3,-1,'k')

hold off

subplot(2,2,4);plot(w*N/(2*pi),phaH*180/pi,'--');

axis([-0.1,8.1,-200,200]);

hold on

stem(k,phaX);

ylabel('度');title('DFT的相角:N=8');text(4.3,-200,'k')

7--2

clear ;close all;clc;

x=[2,2,2,2,2,2,2,2];w=[0:1:500]*2*pi/500;

[H]=freqz(x,1,w);

magH=abs(H);phaH=angle(H);

subplot(2,2,1);plot(w/pi,magH);grid

xlabel('');ylabel('|x|');title('DTFT的幅度')

subplot(2,2,2);plot(w/pi,phaH/pi*180);grid

xlabel('以pi为单位的频率');ylabel('度');title('DTFT的相角') N=16;w1=2*pi/N;k=0:N-1;

X=fft(x,N);

magX=abs(X),phaX=angle(X)*180/pi

subplot(2,2,3);plot(w*N/(2*pi),magH,'--');

axis([-0.1,16.1,0,20]);

hold on

stem(k,magX);

ylabel('|x(k)|');title('DFT的幅度:N=16');text(4.3,-1,'k')

hold off

subplot(2,2,4);plot(w*N/(2*pi),phaH*180/pi,'--');

axis([-0.1,16.1,-200,200]);

hold on

stem(k,phaX,'.');

ylabel('度');title('DFT的相角:N=16');text(4.3,-250,'k')

8--1

clear ;close all;clc;

N=12;w1=2*pi/N;k=0:N-1;

x=[1,2,3,4,5,6,6,5,4,3,2,1];

X=dft(x,N);

magX=abs(X),phaX=angle(X)*180/pi

subplot(2,1,1);

axis([-0.1,12.1,0,50]);

hold on

stem(k,magX);

ylabel('|x(k)|');title('DFT的幅度:N=12');

hold off

subplot(2,1,2);

axis([-0.1,12.1,-400,400]);

hold on

stem(k,phaX);

ylabel('度');title('DFT的相角:N=12');

8--2

clear ;close all;clc;

x=[1,2,3,4,5,6,6,5,4,3,2,1];w=[0:1:500]*2*pi/500;

[H]=freqz(x,1,w);

magH=abs(H);phaH=angle(H);

subplot(2,2,1);plot(w/pi,magH);grid

axis([0,2,0,50]);

xlabel('');ylabel('|x|');title('DTFT的幅度')

subplot(2,2,2);plot(w/pi,phaH/pi*180);grid

axis([0,2,-400,400]);

xlabel('以pi为单位的频率');ylabel('度');title('DTFT的相角')

N=12;w1=2*pi/N;k=0:N-1;

x=[1,2,3,4,5,6,6,5,4,3,2,1];

X=dft(x,N);

magX=abs(X),phaX=angle(X)*180/pi

subplot(2,2,3);plot(w*N/(2*pi),magH,'--');

axis([-0.1,12.1,0,50]);

hold on

stem(k,magX);

ylabel('|x(k)|');title('DFT的幅度:N=12');

hold off

subplot(2,2,4);plot(w*N/(2*pi),phaH*180/pi,'--');

axis([-0.1,12.1,-400,400]);

hold on

stem(k,phaX);

ylabel('度');title('DFT的相角:N=12');

9

clear ;close all;clc;

N1=40;

n=0:1:N1-1;

t=0.01*n;

x=2*sin(4*pi*t)+5*cos(16*pi*t);

x1=fft(x);

magx1=abs(x1);

w=2*pi/N1*n;

subplot(3,1,1);

plot((w*100)/(2*pi),magx1);title('DFT幅度N=40 ');

axis([0,25,0,200]);

N2=60;

n=0:1:N2-1;

t=0.01*n;

x=2*sin(4*pi*t)+5*cos(16*pi*t);

x2=fft(x);

magx2=abs(x2);

w=2*pi/N2*n;

subplot(3,1,2);

plot((w*100)/(2*pi),magx2);title('DFT幅度N=60');

axis([0,25,0,200]);

N3=128;

n=0:1:N3-1;

t=0.01*n;

x=2*sin(4*pi*t)+5*cos(16*pi*t);

x3=fft(x);

magx3=abs(x3);

w=2*pi/N3*n;

subplot(3,1,3);

plot((w*100)/(2*pi),magx3);title('DFT幅度N=128');

axis([0,25,0,400]);

10

clear ;close all;clc;

N=128;

n=0:1:N-1;

t=0.01*n;

x=2*sin(4*pi*t)+5*cos(16*pi*t);

y=x+0.8*randn(1,length(t));

x1=fft(x);

magx1=abs(x1);

w=2*pi/N*n;

subplot(2,1,1); plot((w*100)/(2*pi),magx1);title('DFT幅度');

axis([0,40,0,400]);

y1=fft(y);

magy1=abs(y1);

w=2*pi/N*n;

subplot(2,1,2);plot((w*100)/(2*pi),magy1);title('被噪声污染后DFT幅度') axis([0,100,0,400]);

11

clear ;close all;clc;

N=512;n=0:N-1;t=0.01*n;

x=sin(2*pi*5*t)+sin(2*pi*15*t)+sin(2*pi*30*t); X=fft(x,N);

magx=abs(X);

k=[0:1:N-1];w=2*pi/N*k;

plot(k/N*100,magx);

title('FFT N=512')

xlabel('频率（单位：Hz)');

ylabel('|X|');grid

axis([0,100,0,300])

12

clear ;close all;clc;

N1=128;

n1=0:N1-1;

t1=0.01*n1;

x1=0.5*sin(2*pi*15*t1)+2*sin(2*pi*40*t1);

k1=0:1:127;

w1=2*pi/N1*k1;

X1=fft(x1);

magX1=abs(X1);

subplot(2,1,1);

plot((w1*100)/(2*pi),magX1);

axis([0,50,0,150]);

title('DFT N=128');

xlabel('频率（单位：pi）');

ylabel('X(k)');

grid;

N2=1024;

n2=0:N2-1;

t2=0.01*n2;

x2=0.5*sin(2*pi*15*t2)+2*sin(2*pi*40*t2);

k2=0:1:1023;

w2=2*pi/N2*k2;

X2=fft(x2);

magX2=abs(X2);

subplot(2,1,2);

plot((w2*100)/(2*pi),magX2);

axis([0,50,0,900]);

title('DFT N=1024');

xlabel('频率（单位：pi）');

ylabel('X(k)');

grid;

13

clear ;close all;clc;

t=0:0.001:1;

x=sin(2*pi*60*t)+sin(2*pi*200*t); subplot(2,1,1);plot(t,x);

title('signial x(n)');grid;

y=x+1.5*randn(1,length(t));

Y=fft(y,1024);

p=Y.*conj(Y)/1024;

N=1:1024;n=N/1000*1024; subplot(2,1,2);plot(n,p);

axis([0,600,0,280]);

title('signial y(n)');grid;

xlabel('频率（单位：Hz)'); ylabel('p');grid;

14

clear ;close all;clc;

n=[0:1:9];

x=cos(0.48*pi*n)+cos(0.52*pi*n); X=fft(x);

magx=abs(X(1:1:10));

k=0:1:9;

w=2*pi/10*k;

subplot(3,1,1);

stem(w/pi,magx);

title('N=10点DFT幅度');

xlabel('频率(单位:pi)');

axis([0,1,0,10]);

n=[0:1:9];

y=cos(0.48*pi*n)+cos(0.52*pi*n); n1=[0:1:99];

x=[y(1:1:10) zeros(1,90)];

x1=fft(x);

magx1=abs(x1(1:1:50));

k1=0:1:49;

w1=2*pi/100*k1;

subplot(3,1,2);

stem(w1/pi,magx1);

title('补零到一百点DFT幅度'); xlabel('频率(单位:pi)');

axis([0,1,0,10]);

n=[0:1:99];

x=cos(0.48*pi*n)+cos(0.52*pi*n); X=fft(x);

magx=abs(X(1:1:50));

k=0:1:49;

w=2*pi/100*k;

subplot(3,1,3);

stem(w/pi,magx);

title('N=100点DFT幅度'); xlabel('频率(单位:pi)');

axis([0,1,0,60]);

15

clear ;close all;clc;

n=0:10;

x=10*(0.8.^n);

x1=fft(x);

k=0:10;

y1=x1.*(exp(8*j*pi*k/11));

y=ifft(y1);

subplot(2,2,1);

stem(n,x);

title('原序列x(n)');

xlabel('n');

axis([0,10,0,12]);

subplot(2,2,2);

stem(n,y);

title('移位序列y(n)');

axis([0,10,0,12]);

n=0:10;

y=10*(0.8.^n);

x=[y(1:1:11) zeros(1,4)];

n1=0:14;

subplot(2,2,3);

stem(n1,x);

title('15点序列x(n)');

xlabel('n');

axis([0,14,0,12]);

x1=fft(x);

k=0:14;

y1=x1.*(exp(8*j*pi*k/15));

y=ifft(y1);

subplot(2,2,4);

stem(n1,y);

title('15点移位序列y(n)');

axis([0,14,0,12]);

16

clear ;close all;clc;

N=31;

n=[0:N];

x=n.*(stepseq(0,0,N)-stepseq(16,0,N));

y=stepseq(0,0,N)-stepseq(8,0,N);

X=fft(x);

Y=fft(y);

Z=X.*Y;

z=ifft(Z);

subplot(3,1,1);

stem(n,z);

title('线性卷积');

axis([0,25,0,100]);

N1=15;

n1=[0:31];

x1=n1.*(stepseq(0,0,31)-stepseq(16,0,31)); y1=stepseq(0,0,N1)-stepseq(8,0,N1);

X1=fft(x1,16);

Y1=fft(y1);

Z1=X1.*Y1;

z1=ifft(Z1);

subplot(3,1,2);

n1=[0:15];

stem(n1,z1);

title('16点圆周卷积');

axis([0,20,0,100]);

N=31;

n=[0:N];

x=n.*(stepseq(0,0,N)-stepseq(16,0,N));

y=stepseq(0,0,31)-stepseq(8,0,31);

X=fft(x);

Y=fft(y);

Z=X.*Y;

z=ifft(Z);

subplot(3,1,3);

stem(n,z);

title('32点圆周卷积');

axis([0,25,0,100]);

17

clear ;close all;clc;

Rp=0.5;

T=0.001;

ws=200*2*pi*T;

ws1=(2/T)*tan(ws/2);

[b,a]=cheby1(9,Rp,ws1,'high','s');

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az); subplot(2,1,1);

plot(w/pi,db);

grid;

axis([0,1,-400,100]);

title('系统的幅频响应');

subplot(2,1,2);

plot(w/pi,pha);

title('系统的相频响应');

18

clear ;close all;clc;

Wn=2*pi*100;

fs=1000;

[b,a]=butter(6,Wn,'s');

[bz,az]=impinvar(b,a,fs);

[db,mag,pha,grd,w]=freqz_m(bz,az); subplot(2,2,1);

plot(w/pi,db);

title('系统的幅频响应');

axis([0,1,-50,5]);

subplot(2,2,2);

plot(w/pi,pha);

title('系统的相频响应');

%Filter

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,...

-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,...

-4,2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,...

-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];

y=filter(bz,az,x);

N=56;

n=0:N-1;

subplot(2,2,3);

plot(n,x);

title('输入波形');

subplot(2,2,4);

plot(n,y);

title('输出波形');

19

%最高f=30Hz，可取fs=100Hz，即t=0.01n

%s(n)=sin(0.1*pi*n)+sin(0.3*pi*n)+sin(0.6*pi*n);

%s(n)的样本取301点

%注意：这不是双线性变换法，是完全设计法，不过，效果一样。参考P clear ;close all;clc;

Rp=0.1;

Rs=40;

wp1=0.2*pi;

wp2=0.4*pi;

wn=[wp1,wp2]/pi;

n=4;

[b,a]=ellip(n,Rp,Rs,wn);%默认时表带通

[db,mag,pha,grd,w]=freqz_m(b,a);

subplot(3,1,1);

plot(w/pi,db);

axis([0,1,-100,5]);

n=0:300;

s=sin(0.1*pi*n)+sin(0.3*pi*n)+sin(0.6*pi*n);

subplot(312);

plot(n,s)

hold on

y=filter(b,a,s);

subplot(313);

plot(n,y)

20

%可抽象成一低通或带阻滤波器。抽象成低通来设计

%抽样频率取fs=1000Hz；

%验证看指标

clear ;close all;clc;

fp=100;

fs=130;

Rp=2;

Rs=50;

T=0.001;

wp=2*pi*fp*T;

ws=2*pi*fs*T;

wp1=(2/T)*tan(wp/2);

ws1=(2/T)*tan(ws/2);

[n,wn]=cheb1ord(wp1,ws1,Rp,Rs,'s');

[b,a]=cheby1(n,Rp,wn,'low','s');

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);

grid on

plot(w/pi,db);

axis([0,1,-80,5]);

21

clear ;close all;clc;

b=[1,1];a=[1,5,6];T=1;

[bz,az]=impinvar(b,a,1/T)

[bz1,az1]=bilinear(b,a,1/T)

22

clear ;close all;clc;

Rp=2;Rs=30;T=0.001;

wp1=2*pi*100*T;wp2=2*pi*250*T;ws1=2*pi*50*T;ws2=2*pi*300*T;

wp3=(2/T)*tan(wp1/2);wp4=(2/T)*tan(wp2/2);

ws3=(2/T)*tan(ws1/2);ws4=(2/T)*tan(ws2/2);

wp=[wp3,wp4];ws=[ws3,ws4];

[n,wn]=cheb1ord(wp,ws,Rp,Rs,'s');[z,p,k]=cheb1ap(n,Rp);[b,a]=zp2tf(z,p,k); w0=sqrt(wp3*wp4);Bw=wp4-wp3;

[b1,a1]=lp2bp(b,a,w0,Bw);

[bz,az]=bilinear(b1,a1,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi/T/2,db);axis([0,400,-50,2]); 23

clear ;close all;clc;

Rp=3;

Rs=18;

fs=2000;

ws1=0.2*pi;

wp1=0.3*pi;

wp2=0.4*pi;

ws2=0.5*pi;

wp3=(2*fs)*tan(wp1/2);

wp4=(2*fs)*tan(wp2/2);

ws3=(2*fs)*tan(ws1/2);

ws4=(2*fs)*tan(ws2/2);

wp=[wp3,wp4];

ws=[ws3,ws4];

[n,wn]=buttord(wp,ws,Rp,Rs,'s');

[b,a]=butter(n,wn,'bandpass','s');

[bz,az]=bilinear(b,a,fs);

[db,mag,pha,grd,w]=freqz_m(bz,az);

plot((w*2000)/(2*pi),db);

grid;

24

clear ;close all;clc;

Rp=1;Rs=15;fs=2000;

wp1=0.2*pi;ws1=0.3*pi;

wp2=(2*fs)*tan(wp1/2);

ws2=(2*fs)*tan(ws1/2);

[n,wn]=cheb1ord(wp2,ws2,Rp,Rs,'s');

[b,a]=cheby1(n,Rp,wn,'s');

[bz,az]=bilinear(b,a,fs);

[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);

25

clear ;close all;clc;

Rp=1;Rs=25;wp1=0.2*pi;ws1=0.4*pi;T=0.001;

wp=(2/T)*tan(wp1/2);ws=(2/T)*tan(ws1/2);

[n,wn]=cheb2ord(wp,ws,Rp,Rs,'s') %从此处可以计算阶数n [b,a]=cheby2(n,Rs,wn,'low','s') %由b,a的值可以得到系统函数[bz,az]=bilinear(b,a,1/T);

[b0,B,A]=dir2cas(bz,az)

[db,mag,pha,grd,w]=freqz_m(bz,az);

subplot(2,1,1);plot(w/pi,db);

axis([0,1,-65,10]);

xlabel('');

ylabel('');

title('幅频相应')

subplot(2,1,2);

plot(w/pi,pha);

xlabel('频率（单位：pi）');

ylabel('相位');title('相频相应');

26

clear ;close all;clc;

Rp=1;

Rs=15;

wp1=0.2*pi;

ws1=0.3*pi;

T=0.001;

wp=(2/T)*tan(wp1/2);

ws=(2/T)*tan(ws1/2);

[n,wn]=ellipord(wp,ws,Rp,Rs,'s') %从此处可以计算阶数n

[b,a]=ellip(n,Rp,Rs,wn,'low','s') %由b,a的值可以得到系统函数[bz,az]=bilinear(b,a,1/T);

[b0,B,A]=dir2cas(bz,az)

[db,mag,pha,grd,w]=freqz_m(bz,az);

subplot(2,1,1);

plot(w/pi,db);

axis([0,1,-65,10]);

xlabel('');

ylabel('幅频');

title('幅频相应(相对幅度)')

subplot(2,1,2);

plot(w/pi,pha);

xlabel('频率（单位：pi）');

ylabel('相位');

title('相频相应')

27

clear ;close all;clc;

Rp=1.5;

Rs=20;

T=0.001;

wp=0.4*pi;

ws=0.6*pi;

wp1=(2/T)*tan(wp/2);

ws1=(2/T)*tan(ws/2);

[n,wn]=buttord(wp1,ws1,Rp,Rs,'s')

[b,a]=butter(n,wn,'s')

[bz,az]=bilinear(b,a,1/T);

[b0,B,A]=dir2par(bz,az)

[db,mag,pha,grd,w]=freqz_m(bz,az);

subplot(2,1,1);

plot(w/pi,db);

axis([0,1,-200,5]);

grid

subplot(2,1,2);

plot(w/pi,pha/pi);

axis([0,1,-1.2,1.2]);

grid

28

clear ;close all;clc;

Rp=1;Rs=15;T=0.001;

wp=0.6*pi;ws=0.4*pi;

wp1=(2/T)*tan(wp/2);ws1=(2/T)*tan(ws/2);

[n,wn]=cheb1ord(wp1,ws1,Rp,Rs,'s')

[b,a]=cheby1(n,Rp,wn,'high','s')

[bz,az]=bilinear(b,a,1/T);

[b0,B,A]=dir2par(b,a)

[db,mag,pha,grd,w]=freqz_m(bz,az);

subplot(2,1,1);plot(w/pi,db);axis([0,1,-300,5]);grid

subplot(2,1,2);plot(w/pi,pha/pi);axis([0,1,-1.2,1.2]);grid

29

clear ;close all;clc;

Rp=1;Rs=40;fs=1000;

ws1=0.25*pi;ws2=0.8*pi;

wp1=0.4*pi;wp2=0.7*pi;

wp3=(2*fs)*tan(wp1/2);

wp4=(2*fs)*tan(wp2/2);

ws3=(2*fs)*tan(ws1/2);

ws4=(2*fs)*tan(ws2/2);

wp=[wp3,wp4];

ws=[ws3,ws4];

[n,wn]=cheb2ord(wp,ws,Rp,Rs,'s')

[b,a]=cheby2(n,Rs,wn,'bandpass','s');

[bz,az]=bilinear(b,a,fs);

[db,mag,pha,grd,w]=freqz_m(bz,az);

subplot(2,1,1);plot(w/pi,db);axis([0,1,-100,10]);

subplot(2,1,2);plot(w/pi,pha);

30

clear ;close all;clc;

Rp=1;Rs=60;wp1=0.4*pi;ws1=0.5*pi;T=0.001;

wp=(2/T)*tan(wp1/2);ws=(2/T)*tan(ws1/2);

[n,wn]=ellipord(wp,ws,Rp,Rs,'s') ;

[b,a]=ellip(n,Rp,Rs,wn,'low','s') ;

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);

subplot(3,1,1);plot(w/pi,db);

axis([0,1,-100,5]);

xlabel('');ylabel('相对幅度');

title('幅频相应(相对幅度)')

subplot(3,1,3);

plot(w/pi,pha);xlabel('频率（单位：pi）');ylabel('相位');title('相频相应') 31

clear ;close all;clc;

Rp=1;Rs=50;ws1=0.4*pi;ws2=0.48*pi;fs=1000;

ws3=(2*fs)*tan(ws1/2);ws4=(2*fs)*tan(ws2/2);

n=10;

[z,p,k]=ellipap(n,Rp,Rs);[b,a]=zp2tf(z,p,k);

w0=sqrt(ws3*ws4);bw=ws4-ws3;

[b1,a1]=lp2bs(b,a,w0,bw);

[bz,az]=bilinear(b1,a1,fs);

[db,mag,pha,grd,w]=freqz_m(bz,az);

subplot(3,1,1);

plot(w/pi,db);

axis([0,1,-100,2]);

n=0:200;

x=sin(0.44*pi*n);

subplot(3,1,2);

plot(n,x);hold on

subplot(3,1,3);

y=filter(bz,az,x);

plot(n,y);%输出信号逐渐衰减到0，说明w=0.44pi的信号被滤除

32

clear ;close all;clc;

Rp=0.5;Rs=60;T=1/200;

wp1=60*2*pi*T;wp2=80*2*pi*T;ws1=55*2*pi*T;ws2=85*2*pi*T;

wp3=(2/T)*tan(wp1/2);wp4=(2/T)*tan(wp2/2);

ws3=(2/T)*tan(ws1/2);ws4=(2/T)*tan(ws2/2);

wp=[wp3,wp4];ws=[ws3,ws4];

[n,wn]=cheb2ord(wp,ws,Rp,Rs,'s');

[b,a]=cheby2(n,Rs,wn,'bandpass','s');

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);plot((w*200)/(2*pi),db);axis([0,120,-100,2]);

33

clear ;close all;clc;

Rp=0.8;Rs=25;T=0.001;fp=300;fs=200;

wp=2*pi*fp*T;ws=2*pi*fs*T;

wp1=(2/T)*tan(wp/2);ws1=(2/T)*tan(ws/2);

[n,wn]=buttord(wp1,ws1,Rp,Rs,'s');

[b,a]=butter(n,wn,'high','s');

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);plot(0.5*w/(T*pi),db);axis([0,500,-300,5]);grid 34

clear ;close all;clc;

Rp=0.3;Rs=15;T=1/20000;

wp=2*pi*1400*T;ws=2*pi*1000*T;

wp1=(2/T)*tan(wp/2);ws1=(2/T)*tan(ws/2);

[n,wn]=cheb2ord(wp1,ws1,Rp,Rs,'s');

[b,a]=cheby2(n,Rs,wn,'high','s');

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);plot((w*20000)/(2*pi),db);axis([0,1800,-65,5]);grid 35

clear ;close all;clc;

Rp=0.5;Rs=50;fp1=10;fp2=20;fs1=5;fs2=25;T=0.01;

wp1=2*pi*fp1*T;wp2=2*pi*fp2*T;

wp3=(2/T)*tan(wp1/2);wp4=(2/T)*tan(wp2/2);

wp=[wp3,wp4];

ws1=2*pi*fs1*T;ws2=2*pi*fs2*T;

ws3=(2/T)*tan(ws1/2);ws4=(2/T)*tan(ws2/2);

ws=[ws3,ws4];

[n,wn]=ellipord(wp,ws,Rp,Rs,'s');

[b,a]=ellip(n,Rp,Rs,wn,'s');

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);

plot(w/(pi*2*T),db);axis([0 40 -100 5]);

36

clear ;close all;clc;

Rp=3;Rs=30;fp1=100;fp2=250;fs1=50;fs2=300;T=0.001;

wp1=2*pi*fp1*T;wp2=2*pi*fp2*T;

wp3=(2/T)*tan(wp1/2);wp4=(2/T)*tan(wp2/2);

wp=[wp3,wp4];

ws1=2*pi*fs1*T;ws2=2*pi*fs2*T;

ws3=(2/T)*tan(ws1/2);ws4=(2/T)*tan(ws2/2);

ws=[ws3,ws4];

[n,wn]=ellipord(wp,ws,Rp,Rs,'s');

[b,a]=ellip(n,Rp,Rs,wn,'s');

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);

plot(w/(pi*2*T),db);axis([0 400 -80 5]);

37

clear ;close all;clc;

Rp=3;Rs=30;fp1=50;fp2=300;fs1=100;fs2=250;T=0.001;

wp1=2*pi*fp1*T;wp2=2*pi*fp2*T;

wp3=(2/T)*tan(wp1/2);wp4=(2/T)*tan(wp2/2);

wp=[wp3,wp4];

ws1=2*pi*fs1*T;ws2=2*pi*fs2*T;

ws3=(2/T)*tan(ws1/2);ws4=(2/T)*tan(ws2/2);

ws=[ws3,ws4];

[n,wn]=ellipord(wp,ws,Rp,Rs,'s');

[b,a]=ellip(n,Rp,Rs,wn,'stop','s');

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);

plot(w/(pi*2*T),db);axis([0 400 -90 5]);

38

clear ;close all;clc;

Rp=1;Rs=40;T=0.001;

wp1=0.4*pi;wp2=0.6*pi;

wp3=(2/T)*tan(wp1/2);wp4=(2/T)*tan(wp2/2);

wp=[wp3,wp4];

ws1=0.3*pi;ws2=0.75*pi;

ws3=(2/T)*tan(ws1/2);ws4=(2/T)*tan(ws2/2);

ws=[ws3,ws4];

[n,wn]=ellipord(wp,ws,Rp,Rs,'s');

[b,a]=ellip(n,Rp,Rs,wn,'s');

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);

plot(w/pi,db);axis([0 1 -90 5]);

39

clear ;close all;clc;

h=[-4,1,-1,-2,5,6,6,5,-2,-1,1,-4];

M=length(h);

L=M/2;

b=2*h(L:-1:1);

n=[1:L];n=n-0.5;

w=[0:500]'*2*pi/500

Hr=cos(w*n)*b';

subplot(1,2,1);plot(w/pi,Hr);xlabel('w/pi');title('幅度特性'); subplot(1,2,2);plot(w/pi,20*log(abs(Hr)));axis([0 2 -60 70]); xlabel('w/pi');title('频域特性');

40

clear ;close all;clc;

wc=0.25*pi;N=21;

n=[0:1:N-1];

hd=ideal_lp(wc,N);

w_han=(boxcar(N))';%矩形窗

h=hd.*w_han;

[db,mag,pha,grd,w]=freqz_m(h,1);

subplot(2,1,1);stem(n,h,'.');title('h(n)');grid;hold on;

subplot(212);plot(w/pi,db,'y');xlabel('w/pi');ylabel('dB');hold on; w_han=(hamming(N))';%海明窗

h=hd.*w_han;

[db,mag,pha,grd,w]=freqz_m(h,1);

subplot(2,1,1);stem(n,h,'*');title('h(n)');grid;hold on;

subplot(212);plot(w/pi,db,'b');xlabel('w/pi');ylabel('dB');hold on; w_han=(hanning(N))';%汉宁窗

h=hd.*w_han;

[db,mag,pha,grd,w]=freqz_m(h,1);

subplot(2,1,1);stem(n,h,'.');title('h(n)');grid;hold on;

subplot(212);plot(w/pi,db,'r');xlabel('w/pi');ylabel('dB');hold on; w_han=(blackman(N))';%布莱克曼窗

h=hd.*w_han;

[db,mag,pha,grd,w]=freqz_m(h,1);

subplot(2,1,1);stem(n,h,'.');title('h(n)');grid;hold on;

subplot(212);plot(w/pi,db,'g');xlabel('w/pi');ylabel('dB');hold on;

41

clear ;close all;clc;

Rp=3;Rs=20;

N=15;n=[0:1:N-1];%只需要改变N的值

wc=pi/4;

hd=ideal_lp(wc,N);

w_han=(hanning(N))';

h=hd.*w_han;

[db,mag,pha,grd,w]=freqz_m(h,1);

subplot(211);plot(w/pi,db);xlabel('w/pi');ylabel('db');

hold on;

subplot(212);plot(w/pi,pha);xlabel('w/pi');ylabel('pha');

hold on;

N=33;n=[0:1:N-1];%只需要改变N的值

wc=pi/4;

hd=ideal_lp(wc,N);

w_han=(hanning(N))';

h=hd.*w_han;

[db,mag,pha,grd,w]=freqz_m(h,1);

subplot(211);plot(w/pi,db,'r');xlabel('w/pi');ylabel('db'); subplot(212);plot(w/pi,pha,'r');xlabel('w/pi');ylabel('pha');

42

clear ;close all;clc;

Rp=3;Rs=20;

N=33;n=[0:1:N-1];

wc=pi/4;

数字信号处理Matlab实现实例(推荐给学生)

数字信号处理Matlab 实现实例 第１章离散时间信号与系统 例1-1 用MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。 解 MATLAB程序如下： a=[-2 0 1 -1 3]; b=[1 2 0 -1]; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel('n'); ylabel('幅度'); 图1.1给出了卷积结果的图形，求得的结果存放在数组c中为：{-2 -4 1 3 1 5 1 -3}。 例1-2 用MATLAB计算差分方程 当输入序列为时的输出结果。 解 MATLAB程序如下： N=41; a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)];

k=0:1:N-1; y=filter(a,b,x); stem(k,y) xlabel('n');ylabel('幅度') 图 1.2 给出了该差分方程的前41个样点的输出，即该系统的单位脉冲响应。 例1-3 用MATLAB 计算例1-2差分方程 所对应的系统函数的DTFT 。 解 例1-2差分方程所对应的系统函数为： 123 123 0.80.440.360.02()10.70.450.6z z z H z z z z -------++= +-- 其DTFT 为 23230.80.440.360.02()10.70.450.6j j j j j j j e e e H e e e e ωωωω ωωω--------++= +-- 用MATLAB 计算的程序如下： k=256; num=[0.8 -0.44 0.36 0.02]; den=[1 0.7 -0.45 -0.6]; w=0:pi/k:pi; h=freqz(num,den,w); subplot(2,2,1); plot(w/pi,real(h));grid title('实部') xlabel('\omega/\pi');ylabel('幅度')

数字信号处理MATLAB中FFT实现

MATLAB中FFT的使用方法 说明：以下资源来源于《数字信号处理的MATLAB实现》万永革主编 一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MATLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[43267890]; Xk=fft(xn) → Xk= 39.0000-10.7782+6.2929i0-5.0000i 4.7782-7.7071i 5.0000 4.7782+7.7071i0+5.0000i-10.7782-6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。

clf; fs=100;N=128;%采样频率和数据点数 n=0:N-1;t=n/fs;%时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%信号 y=fft(x,N);%对信号进行快速Fourier变换 mag=abs(y);%求得Fourier变换后的振幅 f=n*fs/N;%频率序列 subplot(2,2,1),plot(f,mag);%绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2));%绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%信号 y=fft(x,N);%对信号进行快速Fourier变换 mag=abs(y);%求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag);%绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2));%绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果：

MATLAB数字信号处理

MATLAB 下的数字信号处理实现示例 附录一信号、系统和系统响应 1、理想采样信号序列 （1）首先产生信号x(n),0<=n<=50 n=0:50; A=444.128; a=50*sqrt(2.0)*pi; T=0.001; w0=50*sqrt(2.0)*pi; x=A*exp(-a*n*T).*sin(w0*n*T); close all subplot(3,1,1);stem(x); %定义序列的长度是50 %设置信号有关的参数 %采样率 %pi 是MATLAB 定义的π，信号乘可采用“.*”%清除已经绘制的x(n)图形 %绘制x(n)的图形 title(…理想采样信号序列?); （2）绘制信号x(n)的幅度谱和相位谱 k=-25:25; W=(pi/12.5)*k; X=x*(exp(-j*pi/12.5)).^(n?*k); magX=abs(X); %绘制x(n)的幅度谱subplot(3,1,2);stem(magX);title(…理想采样信号序列的幅度谱?); angX=angle(X); %绘制x(n)的相位谱subplot(3,1,3);stem(angX) ; title (…理想采样信号序列的相位谱?) （3）改变参数为：A = 1,? = 0.4, & 0 = 2.0734, T = 1 n=0:50; A=1; a=0.4; T=1; w0=2.0734; x=A*exp(-a*n*T).*sin(w0*n*T); close all subplot(3,1,1);stem(x); title(…理想采样信号序列?); k=-25:25; W=(pi/12.5)*k; X=x*(exp(-j*pi/12.5)).^(n?*k); magX=abs(X); %定义序列的长度是50 %设置信号有关的参数 %采样率 %pi 是MATLAB 定义的π，信号乘可采用“.*” %清除已经绘制的x(n)图形 %绘制x(n)的图形 %绘制x(n)的幅度谱 subplot(3,1,2);stem(magX);title(…理想采样信号序列的幅度谱?); angX=angle(X); %绘制x(n)的相位谱

实验一 基于Matlab的数字信号处理基本

实验一 基于Matlab 的数字信号处理基本操作 一、 实验目的：学会运用MA TLAB 表示的常用离散时间信号；学会运用MA TLAB 实现离 散时间信号的基本运算。 二、 实验仪器：电脑一台，MATLAB6.5或更高级版本软件一套。 三、 实验内容： （一） 离散时间信号在MATLAB 中的表示 离散时间信号是指在离散时刻才有定义的信号，简称离散信号，或者序列。离散序列通常用)(n x 来表示，自变量必须是整数。 离散时间信号的波形绘制在MATLAB 中一般用stem 函数。stem 函数的基本用法和plot 函数一样，它绘制的波形图的每个样本点上有一个小圆圈，默认是空心的。如果要实心，需使用参数“fill ”、“filled ”，或者参数“.”。由于MATLAB 中矩阵元素的个数有限，所以MA TLAB 只能表示一定时间范围内有限长度的序列；而对于无限序列，也只能在一定时间范围内表示出来。类似于连续时间信号，离散时间信号也有一些典型的离散时间信号。 1. 单位取样序列 单位取样序列)(n δ，也称为单位冲激序列，定义为 ) 0() 0(0 1)(≠=?? ?=n n n δ 要注意，单位冲激序列不是单位冲激函数的简单离散抽样，它在n =0处是取确定的值1。在MATLAB 中，冲激序列可以通过编写以下的impDT .m 文件来实现，即 function y=impDT(n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 【实例1-1】 利用MATLAB 的impDT 函数绘出单位冲激序列的波形图。 解：MATLAB 源程序为 >>n=-3:3; >>x=impDT(n); >>stem(n,x,'fill'),xlabel('n'),grid on >>title('单位冲激序列') >>axis([-3 3 -0.1 1.1]) 程序运行结果如图1-1所示。 图1-1 单位冲激序列

青岛理工大学临沂年数字信号处理及MATLAB试卷

A卷

一、[15分] 1、10 2、f>=2fh

3、()()()y n x n h n =* 4、1 -az -11a 或者-z z ，a 1 -z 或1-1-az -1z 5、对称性 、 可约性 、 周期性 6、191点，256 7、典范型、级联型、并联型 8、T ω = Ω，)2 tan(2ω T = Ω或)2arctan(2T Ω=ω。 二、[20分] 1、C 2、 A 3、 C 4、C 5、B 6、D 7、B 8、A 9、D 10、A (CACCB DBADA) 三、[15分] 1、（5分） 混叠失真：不满足抽样定理的要求。 改善方法：增加记录长度 频谱泄漏：对时域截短，使频谱变宽拖尾，称为泄漏 改善方法：1）增加w (n )长度 2）缓慢截短 栅栏效应：DFT 只计算离散点（基频F0的整数倍处）的频谱，而不是连续函数。 改善方法：增加频域抽样点数N （时域补零），使谱线更密 2、（5分） 3、 （5分） IIR 滤波器: 1）系统的单位抽样相应h (n )无限长 2）系统函数H (z )在有限z 平面（ ）上有极点存在 3）存在输出到输入的反馈，递归型结构 Fir 滤波器: ? 1）系统的单位冲激响应h (n )在有限个n 处不为零； ? 2）系统函数 在||0 z >处收敛，在 处只有零点，即有限z 平面只有零点，而全部极点都在z =0处； ? 3）机构上主要是非递归结构，没有输入到输出的反馈，但有些结构中也包含有反馈的递归部分。 四、计算题（40分） 1、（12分）解： 解： 对上式两边取Z 变换，得： ()H z ||0z >

数字信号处理指导书matlab版

实验1 时域离散信号的产生 一、实验目的 学会运用MATLAB 产生常用离散时间信号。 二、实验涉及的matlab 子函数 1、square 功能：产生矩形波 调用格式： x=square(t)；类似于sin （t ），产生周期为2*pi ，幅值为+—1的方波。 x=square(t ，duty)；产生制定周期的矩形波，其中duty 用于指定脉冲宽度与整个周期的比例。 2、rand 功能：产生rand 随机信号。 调用格式： x=rand （n ，m ）；用于产生一组具有n 行m 列的随机信号。 三、实验原理 在时间轴的离散点上取值的信号，称为离散时间信号。通常，离散时间信号用x （n ）表示，其幅度可以在某一范围内连续取值。 由于信号处理所用的设备主要是计算机或专用的信号处理芯片，均以有限的位数来表示信号的幅度，因此，信号的幅度也必须“量化”，即取离散值。我们把时间和幅度上均取离散值的信号称为时域离散信号或数字信号。 在MATLAB 中，时域离散信号可以通过编写程序直接生成，也可以通过对连续信号的等间隔抽样获得。 下面介绍常用的时域离散信号及其程序。 1、单位抽样序列 ? ? ?≠==000 1)(k k k δ MATLAB 源程序为

1) function [x,n] = impuls (n0,n1,n2) % Generates x(n) = delta(n-n0); n=n0 处建立一个单位抽样序列% [x,n] = impuls (n0,n1,n2) if ((n0 < n1) | (n0 > n2) | (n1 > n2)) error('arguments must satisfy n1 <= n0 <= n2') end n = [n1:n2]; x = [zeros(1,(n0-n1)), 1, zeros(1,(n2-n0))]; 将上述文件存为：impuls.m，在命令窗口输入 n0=0,n1=-10,n2=11; [x,n]=impuls (n0,n1,n2); stem(n,x,’filled’) 2)n1=-5;n2=5;n0=0; n=n1:n2; x=[n==n0]; stem(n,x,'filled','k'); axis([n1,n2,1.1*min(x),1.1*max(x)]); title('单位脉冲序列'); xlabel('时间(n)'); ylabel('幅度x(n)'); 3)n1=-5;n2=5;k=0; n=n1:n2; nt=length(n); %求n点的个数 nk=abs(k-n1)+1; %确定k在n序列中的位置 x=zeros(1,nt); %对所有样点置0 x(nk)=1; %对抽样点置1 stem(n,x,'filled','k'); axis([n1,n2,0,1.1*max(x)]); title('单位脉冲序列'); xlabel('时间(n)'); Ylabel('幅度x(n)');

数字信号处理的MATLAB实现

昆明理工大学信息工程与自动化学院学生实验报告 （2011—2012 学年第二学期） 课程名称：数字信号处理开课实验室：信自楼111 2012 年 5 月 31 日年级、专业、班生医学号姓 名 成绩 实验项目名称数字信号处理的matlab 实现指导教师 教 师 评语教师签名： 年月日 一．实验目的 熟练掌握matlab的基本操作。 了解数字信号处理的MATLAB实现。 二．实验设备 安装有matlab的PC机一台。 三．实验内容 .1.求信号x(n)=cos(6.3Пn/3)+cos(9.7Пn/30)+cos(15.3Пn/30),0≤n≤29的幅度频谱. 2. 用冲击响应不变法设计一个Butterworth低通数字滤波器,要求参数为: Wp=0.2Пαp=1dB Ws=0.3Пαs=15dB 3.用双线性变换法设计一个Chebyshev高通IIR滤波器,要求参数为: Wp=0.6Пαp=1dB Ws=0.4586Пαs=15dB 4.用窗函数法设计一个低通FIR滤波器,要求参数为: Wp=0.2Пαp=0.3dB Ws=0.25Пαs=50dB 5.用频率抽样法设计一个带通FIR滤波器,要求参数为: W1s=0.2П W1p=0.35П W2p=0.65П W2s=0.8П αs=60dB αp=1dB 6.根据 4 点矩形序列，( n ) = [1 1 1 1] 。做 DTFT 变换，再做 4 点 DFT 变换。然后分别补零做 8 点 DFT 及 16 点 DFT。 7.调用filter解差分方程，由系统对u(n)的响应判断稳定性 8编制程序求解下列系统的单位冲激响应和阶跃响应。 y[n]+ 0.75y[n -1]+ 0.125y[n -2] = x[n]- x[n -1] 四．实验源程序 1. n=[0:1:29]; x=cos(6.3*pi*n/30)+cos(9.7*pi*n/30)+cos(15.3*pi*n/30);

数字信号处理MATLAB实验1

实验一熟悉MATLAB环境 一、实验目的 (1)熟悉MATLAB的主要操作命令。 (2)学会简单的矩阵输入和数据读写。 (3)掌握简单的绘图命令。 (4)用MATLAB编程并学会创建函数。 (5)观察离散系统的频率响应。 二、实验内容 认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。在熟悉了MATLAB基本命令的基础上，完成以下实验。 上机实验内容： (1)数组的加、减、乘、除和乘方运算。输入A=[1234]，B=[345 6]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出 A、B、C、D、E、F、G。 (2)用MATLAB实现以下序列。 a）x(n)=0.8n0≤n≤15 b）x(n)=e(0.2+3j)n0≤n≤15 c）x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π)0≤n≤15 (n)=x(n+16),绘出四个d）将c)中的x(n)扩展为以16为周期的函数x 16 周期。 (n)=x(n+10),绘出四个e）将c)中的x(n)扩展为以10为周期的函数x 10 周期。

(3)x(n)=[1,-1,3,5]，产生并绘出下列序列的样本。 a）x 1(n)=2x(n+2)-x(n-1)-2x(n) b）∑=-=5 1k 2) k n (nx (n) x (4)绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注。 a)x(t)=sin(2πt)0≤t≤10s b)x(t)=cos(100πt)sin(πt) 0≤t≤4s (5)编写函数stepshift(n0,n1,n2)实现u(n-n0),n1

MATLAB实现数字信号处理

《数字信号处理》课程设计实例： 声音信号的处理 一．摘要： 这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB或者ＤＳＰ开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。通过对声音的采样，将声音采样后的频谱与滤波。 MATLAB全称是Matrix Laboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。。经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。 在本次课程设计中，主要通过MATLAB来编程对语音信号处理与滤波，设计滤波器来处理数字信号并对其进行分析。 二．课程设计目的： 综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。 三．设计容： 容：录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性

变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音容，分析两段容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除。 四．设计原理： 4.1．语音信号的采集 熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数，在MATLAB环境中，有关声音的函数有： a：y=wavrecord(N,fs,Dtype);利用系统音频输入设备录音，以fs为采样频率，默认值为11025，即以11025HZ进行采样。Dtype为采样数据的存储格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度进行采样，其它三种都是16位采样结果转换为指定的MATLAB数据； b：wavplay(y,fs)；利用系统音频输出设备播放，以fs为播放频率，播放语音信号y； c：wavwrite((y,fs,wavfile);创建音频文件； d：y=wavread(file);读取音频文件； 关于声音的函数还有sound()；soundsc()；等。 4.2滤波器： 4．21.IIR滤波器原理 冲激响应不变法是使数字滤波器在时域上模拟滤波器，但是它们的缺点是产生频率响应的混叠失真，这是由于从s平面到z平面是多值的映射关系所造成的。 双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，再通过变换关系将此横带变换到整个z平面上去，这样就使得s平面与z平面是一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。

数字信号处理基本知识点Matlab实现

数字信号处理（第二版） 绪论 1.4 MATLAB 在信号处理中的应用简介 MATLAB 是美国Mathworks 公司于1984年推出的一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算、信号处理、系统仿真和图形显示于一体，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。 MATLAB 软件包括五大通用功能：数值计算功能（Numeric ），符号运算功能（Symbolic ）；数据可视化功能（Graphic ），数据图形文字统一处理功能（Notebook ）和建模仿真可视化功能（Simulink ）。该软件有三大特点：一是功能强大；二是界面友善、语言自然；三是开放性强。目前，Mathworks 公司已推出30多个应用工具箱。MA TLAB 在线性代数、矩阵分析、数值及优化、数理统计和随机信号分析、电路与系统、系统动力学、信号和图像处理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统、以及财政金融等众多领域的理论研究和工程设计中得到了广泛应用。 2.10 离散时间信号与系统的Matlab 表示 2.10.1 离散时间信号的表示和运算 1、基本序列的Matlab 表示 单位采样序列 在MA TLAB 中，单位采样序列可以通过编写以下的DTimpulse .m 文件来实现，即 function y=DTimpulse (n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 单位阶跃序列 在MA TLAB 中，单位阶跃序列可以通过编写DTu .m 文件来实现，即 function y=DTu (n) y=n>=0; %当参数为非负时输出1 调用该函数时n 必须为整数或整数向量。 矩形序列 用MA TLAB 表示矩形序列可根据公式()()()N R n u n u n N =--并利用DTu 函数生成，即 function y=DTR(n,N) y=DTu(n)-DTu(n-N); 调用该函数时n 必须为整数或整数向量，N 必须为整数。 实指数序列 用MA TLAB 表示实指数序列()(),n x n a u n n N a R =∈∈，即

南京理工大学数字信号处理matlab上机完美版

1.已知3阶椭圆IIR数字低通滤波器的性能指标为：通带截止频率0.4π，通带波纹为0.6dB，最小阻带衰减为32dB。设计一个6阶全通滤波器对其通带的群延时进行均衡。绘制低通滤波器和级联滤波器的群延时。 %Q1_solution %ellip(N,Ap,Ast,Wp) %N--->The order of the filter %Ap-->ripple in the passband %Ast->a stopband Rs dB down from the peak value in the passband %Wp-->the passband width [be,ae]=ellip(3,0.6,32,0.4); hellip=dfilt.df2(be,ae); f=0:0.001:0.4; g=grpdelay(hellip,f,2); g1=max(g)-g; [b,a,tau]=iirgrpdelay(6,f,[0 0.4],g1); hallpass=dfilt.df2(b,a); hoverall=cascade(hallpass,hellip); hFVT=fvtool([hellip,hoverall]); set(hFVT,'Filter',[hellip,hoverall]); legend(hFVT,'Lowpass Elliptic filter','Compensated filter'); clear; [num1,den1]=ellip(3,0.6,32,0.4); [GdH,w]=grpdelay(num1,den1,512); plot(w/pi,GdH); grid xlabel('\omega/\pi'); ylabel('Group delay, samples'); F=0:0.001:0.4; g=grpdelay(num1,den1,F,2); % Equalize the passband Gd=max(g)-g; % Design the allpass delay equalizer [num2,den2]=iirgrpdelay(6,F,[0,0.4],Gd); [GdA,w] = grpdelay(num2,den2,512); hold on; plot(w/pi,GdH+GdA,'r');

基于MATLAB的数字信号处理

数字信号处理课程设计报告题目：语音数字信号处理与分析及 Matlab实现 系别通信工程 专业班级 学生姓名 学号 指导教师 提交日期

摘要 本次课程设计综合利用数字信号处理的理论知识进行语音信号的频谱分析，通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。本次课程设计要求利用MATLAB对语音信号进行分析和处理，要求学生采集语音信号后，在MATLAB软件平台进行频谱分析；并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。待处理语音信号是一个在20Hz～20kHz 频段的低频信号。采用了高效快捷的开发工具——MATLAB，实现了语音信号的采集，对语音信号加噪声及设计滤波器滤除噪声的一系列工作。利用采样原理设计了高通滤波器、低通滤波器、带通滤波器、带阻滤波器。同学通过查阅资料自己获得程序进行滤波器的设计，能过得到很好的锻炼。 关键词：MATLAB滤波器数字信号处理

目录 第一章绪论 (1) 1.1设计的目的及意义 (1) 1.2设计要求 (1) 1.3设计内容 (1) 第二章系统方案论证 (3) 2.1设计方案分析 (3) 2.2实验原理 (3) 第三章信号频谱分析 (6) 3.1原始信号及频谱分析 (6) 3.2加入干扰噪声后的信号及频谱分析 (7) 第四章数字滤波器的设计与实现 (11) 4.1高通滤波器的设计 (11) 4.2低通滤波器的设计 (12) 4.3带通滤波器的设计 (15) 4.4带阻滤波器的设计 (16) 第五章课程设计总结 (19) 参考文献 (20) 附录Ⅰ..................................................................................I 附录Ⅱ................................................................................II

数字信号处理实验 matlab版 线性相位FIR数字滤波器

实验23 线性相位FIR数字滤波器 (完美格式版，本人自己完成，所有语句正确，不排除极个别错误，特别适用于山大，勿用冰点等工具下载，否则下载之后的word格式会让很多部分格式错误，谢谢) XXXX学号姓名处XXXX 一、实验目的 1 加深对线性相位FIR数字滤波器特性的理解。 2 掌握线性相位滤波器符幅特性和零极点分布的研究方法。 3 了解用MATLAB研究线性相位滤波器特性时程序编写的思路和方法。 二、实验内容 1 线性相位FIR滤波器的特性 2 第一类线性相位滤波器(类型Ⅰ) 3 第二类线性相位滤波器(类型Ⅱ) 4 第三类线性相位滤波器(类型Ⅲ) 5 第四类线性相位滤波器(类型Ⅳ) 6 线性相位FIR数字滤波器零点分布特点 三、实验环境 MATLAB7.0 四、实验原理 1.线性相位FIR滤波器的特性 与IIR滤波器相比，FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到有严格的线性相位特性。设FIR滤波器单位脉冲响应h(n)长度为N，其系统函数为 ∑-=- = 1 N n n z h(n) H(z) 当滤波器的系数N满足一定的对称条件时，就可以获得线性相位。线性相位FIR滤波器共分为四种类型，分别为： (1)类型Ⅰ，系数对称，即h(n)＝h(N－1－n)，N为奇数。 (2)类型Ⅱ，系数对称，即h(n)＝h(N－1－n)，N为偶数。 (3)类型Ⅲ，系数反对称，即h(n)＝－h(N－1－n)，N为奇数。 (4)类型Ⅳ，系数反对称，即h(n)＝－h(N－1－n)，N为偶数。 对于上述四类线性相位FIR滤波器，参考文献［1］中提供了一段通用程序，对考虑正负号的幅度频率特性(简称符幅特性)进行求解，程序名为amplres.m，程序如下：function[A,w,type,tao]=amplres(h) N=length(h);tao=(N-1)/2; L=floor((N-1)/2); n=1:L+1; w=[0:500]*2*pi/500; if all(abs(h(n)-h(N-n+1))<1e-10)

数字信号处理MATLAB仿真(DOC)

实验一 数字信号处理的Matlab 仿真 一、实验目的 1、掌握连续信号及其MATLAB 实现方法； 2、掌握离散信号及其MATLAB 实现方法 3、掌握离散信号的基本运算方法，以及MATLAB 实现 4、了解离散傅里叶变换的MATLAB 实现 5、了解IIR 数字滤波器设计 6、了解FIR 数字滤波器设计1 二、实验设备 计算机，Matlab 软件 三、实验内容 （一）、 连续信号及其MATLAB 实现 1、 单位冲击信号 ()0,0()1,0t t t dt εεδδε-?=≠??=?>??? 例1.1：t=1/A=50时，单位脉冲序列的MATLAB 实现程序如下： clear all; t1=-0.5:0.001:0; A=50; A1=1/A; n1=length(t1); u1=zeros(1,n1); t2=0:0.001:A1;

t0=0; u2=A*stepfun(t2,t0); t3=A1:0.001:1; n3=length(t3); u3=zeros(1,n3); t=[t1 t2 t3]; u=[u1 u2 u3]; plot(t,u) axis([-0.5 1 0 A+2]) 2、 任意函数 ()()()f t f t d τδττ+∞ -∞=-? 例1.2：用MATLAB 画出如下表达式的脉冲序列 ()0.4(2)0.8(1) 1.2() 1.5(1) 1.0(2)0.7(3)f n n n n n n n δδδδδδ=-+-+++++++ clear all; t=-2:1:3; N=length(t); x=zeros(1,N); x(1)=0.4; x(2)=0.8 x(3)=1.2; x(4)=1.5; x(5)=1.0;

数字信号处理及matlab实现部分作业详解

1、把序列?????== ,01 ,20 ,1)(其他＝n n n x 表示为单位阶跃之和的形式。解：) 2(2)1()( )]2()1([2)1()()1(2)()(---+=---+--=-+=n u n u n u n u n u n u n u n n n x δδ2、判断下列系统线性性、因果性、稳定性。(a))()(n nx n y =； (b)b n ax n y +=)()(,其中a ，b 为常数； 解：(a)线性性：对于两个输入序列)(1n x 和)(2n x ，相应的输出分别为)(1)(1n nx n y =) (2)(2n nx n y =这两个输出的线性组合为 ) (2)(1)](2)(1)(3n bnx n anx n by n ay n y +=+=这两个输入信号的线性组合产生的输出为 ) (2)(1)](2)(1[)](2)(1[)(4n bnx n anx n bx n ax n n bx n ax T n y +=+=+=现在)(4)(3n y n y =，所以系统为线性系统； 因果性：因为系统只与当前输入有关，所以系统是因果的； 稳定性：若)(n x 有界，即∞<≤M n x |)(|，则nM n x n n nx n y ≤≤=|)(||)(||)(|，当∞→n 时，∞→)(n y ，所以不稳定。 (b)线性性：对于两个输入序列)(1n x 和)(2n x ，相应的输出分别为 b n ax n y +=)(1)(1b n ax n y +=)(2)(2这两个输出的线性组合为 db n dax cb n cax n dy n cy n y +++=+=)(2)(1)](2)(1)(3这两个输入信号的线性组合产生的输出为 b n dax n cax b n dx n cx a n dx n cx T n y ++=++=+=)(2)(1)](2)(1[)](2)(1[)(4现在)(4)(3n y n y ≠，所以系统为非线性系统； 因果性：因为系统只与当前输入有关，所以系统是因果的； 稳定性：若)(n x 有界，即∞<≤M n x |)(|， 则|||||||)(||)(||)(|b M a b n ax b n ax n y +≤+≤+=，即|)(|n y 有界，所以稳定。 3、一个线性移不变系统的冲击响应)(n h 和输入信号)(n x 分别为 ?????=== 2 ,11 ,20 ,1)(n n n n h ＝，?? ???==-=1 ,30 ,21 ,1)( n n n n x ＝，求系统对输入的响应。 解：根据卷积公式∑-=-= 2 1)()()(k k n h k x n y ，得到?????????====-=3 ,32 ,81 ,80 ,41 ,1)(n n n n n n y ＝

数字信号处理原理及应用的matlab程序

复指数的matlab程序： n=[0:40]; x=exp((0.02+0.3j)*n); subplot(2,2,1); stem(n,real(x)); xlabel('n');ylabel('Re(x)'); subplot(2,2,2); stem(n,imag(x)); xlabel('n');ylabel('Im(x)'); subplot(2,2,3); stem(n,abs(x)); xlabel('n');ylabel('Magnitude'); subplot(2,2,4); stem(n,(180/pi)*angle(x)); xlabel('n');ylabel('Phase');

移位的matlab程序： m=[-1:5]; x=(1/3)*(1/3).^m; n=m-1; x1=x; n1=m+1; x2=x; subplot(3,1,1);stem(m,x);ylabel('x(n)');xlabel('n'); subplot(3,1,2);stem(n,x1);ylabel('x(n+1)');xlabel('n'); subplot(3,1,3);stem(n1,x2);ylabel('x(n-1)');xlabel('n');

反褶的matlab程序： n=[-1:5]; x=(1/3)*(1/3).^n; y=fliplr(x); m=-fliplr(n); subplot(2,1,1);stem(n,x);ylabel('x(n)');xlabel('n'); subplot(2,1,2);stem(m,y);ylabel('x(-n)');xlabel('n'); 加减乘积的matlab程序： m1=[-2:-1];y1=2..^m1; m2=[0:5];y2=m2+1;

数字信号处理的matlab实现

数字信号处理的matlab实现 简谐振动的特性完全取决于振幅、频率、初相位角。 x1 = (0.5).^t; x1 = 0.5* sin(2*pi*f*t+pi/4); x1 = [(n - n0) >= 0]; %阶跃信号 x1 = [(n-n0) == 0]; %脉冲信号 y = sin(pi * x +eps)./(pi * x +eps);%sinc function，eps是matlab系统的精度，这里防止被零除 n = [-10:10]; alpha = -0.1 + 0.3 * j; x = exp(alpha * n);%复指数信号 rx = real(x); ix = imag(x); mx = abs(x); %振幅 px = (180/pi) * angle(x); % 相位，并转换为度 x = rand(1,10); x = randn(1,10); %guass series numbers xp = [x x x x x]; x2 = fliplr(x);%左右折叠 第二章信号 %test sampling rule

dt = 0.01; %samping frequence for draw 100Hz t = 0:dt:1; f = 10; x = sin(2 * pi * f * t + 0.3); dt1 = 0.1; t1 = 0:dt1:1;%10Hz x1 = sin(2 * pi * f * t1 + 0.3); dt2 = 0.05; t2 = 0:dt2:1;%10Hz x2 = sin(2 * pi * f * t2 + 0.3); subplot(3,1,1), plot(t, x); title('f = 10Hz, fs = 100Hz'); subplot(3,1,2), plot(t1, x1), hold on, stem(t1, x1), plot(t, x); title('f = 10Hz, fs = 10Hz'); subplot(3,1,3), plot(t2, x2), hold on, stem(t2, x2), plot(t, x); title('f = 10Hz, fs = 20Hz'); 基本信号： x1 = (0.5).^t; %指数信号 x1 = 0.5* sin(2*pi*f*t+pi/4); %余弦信号 x1 = [(n - n0) >= 0]; %阶跃信号 x1 = [(n-n0) == 0]; %脉冲信号 y = sin(pi * x +eps)./(pi * x +eps);%sinc function，eps是matlab系统的精度，这里防止被零除，sinc函数信号 n = [-10:10];

数字信号处理DFTMATLAB程序

页脚内容1 实验三 频域信号处理 1. 实验目的 （1） 学习信号DFT 变换的matlab 实现； （2） 学习fft 的matlab 实现； （3） 验证DFT 的相关性质。 2. 思考题 （1） 若()()()sin sin 4x n n n ππ=+是一个128点的有限长序列，求其128点DFT 结果； 程序如下： 求DFT 变换矩阵A ： clc; clear; N=128; A=dftmtx(N) Ai=conj(dftmtx(N)); n=0:(N-1);

k=0:(N-1); nk=n'*k; Wn=(sin(pi/8)+sin(pi/4)).^nk Wk=conj(Wn)/N; 求128点的DFT（分别用FFT函数和dftmtx函数） clc; clear; N=128; n=0:N-1; x=sin(pi/8*n)+sin(pi/4*n); subplot(3,1,1) plot(n,x); grid on title('原图') y1=fft(x,N); A=dftmtx(N); 页脚内容2

y2=x(1:N)*A; subplot(3,1,2) plot(n,y1) grid on title('FFT') subplot(3,1,3) plot(n,y2) grid on title('dftmtx') 程序运行结果如图1所示： 原图 -13FFT -13dftmtx 页脚内容3

页脚内容4 图 1 （2） 对模拟信号()()()2sin 45sin 8x t t t ππ=+，以0.01t n =，()0:1n N =-进行采样，求 a ） N =40点的FFT 幅度谱，从图中能否观察出两个频谱分量； b ） 提高采样点数值N=128,再求该信号的幅度频谱，此时幅度频谱发生了什么变化？信号的两个模拟频率和数字频率分别为多少？FFT 频谱分析结果和理论上是否一致？ 程序如下： clc; clear; N=40; n=0:N-1; t=0.01*n x=2*sin(4*pi*t)+5*sin(8*pi*t); subplot(2,1,1) plot(x(1:N)) grid on title('原图') y1=fft(x,N);