MATLAB实验教程

实验一

1.实现单位采样序列，单位阶跃序列，矩形序列，并用图形显示。写出程序及输出图形。代码：

（1）单位采样序列

x=-10:10;

y=[zeros(1,10),1,zeros(1,10)];

figure(3);

stem(x,y);

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

（2）单位阶跃序列

x=-10:10;

y=[zeros(1,10),1,ones(1,10)];

figure(4);

stem(x,y);

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

（3）矩形序列

x=-10:10;

y=[zeros(1,10),1,ones(1,5),zeros(1,5)];

figure(4);

stem(x,y);

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

图：

Figure 8.5.1. 1 单位采样序列Figure 8.5.1. 2 单位阶跃序列

Figure 8.5.1. 3 矩形序列

实验心得：

实验二

1．模拟信号x(t)=2sin(4*pi*t)+5cos(8*pi*t)，以t=0.01n(n=0:N-1)进行采样，求：（1）N=40点FFT的幅度频谱，从图中能否观察出信号的2个频谱分量？

（2）提高采样点数，如N=128，再求该信号的频谱幅度，此时幅度频谱发生了什么变化？代码：

N=40;n=0:N-1;

t=0.01*n;

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

k=0:N/2;w=2*pi/N*k;

X=fft(x,N);

magX=abs(X(1:N/2+1));

subplot(2,1,1);stem(n,x,'.');title('signal x(n)');

subplot(2,1,2);plot(w/pi,magX);title('FFT N=40');

xlabel('f (unit :pi)');ylabel('|X|');grid

N=128;n=0:N-1;

t=0.01*n;

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

k=0:N/2;w=2*pi/N*k;

X=fft(x,N);

magX=abs(X(1:N/2+1));

subplot(2,1,1);stem(n,x,'.');title('signal x(n)');

subplot(2,1,2);plot(w/pi,magX);title('FFT N=128');

xlabel('f (unit :pi)');ylabel('|X|');grid

图：

5

10

15

2025

30

35

40

signal x(n)

50

100

FFT N=40

Figure 8.2.1. 1 N=40 的FFT 频谱图

20

40

60

80

100

120

140

signal x(n)

0.1

0.2

0.3

0.4

0.50.60.7

0.8

0.9

1

010*******

400FFT N=128

f (unit :pi)

|X |

Figure 8.2.1. 2 N=128 的FFT 频谱图

实验心得： 实验六

双线性变换法 代码

%Filter clc;%clear;

%digital filter specifications; lfp=100; lfs=300;

ws=2*pi*lfs;

wp=2*pi*lfp;

Rp=1;

Rs=25;

fs=1000;

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

[z,p,k] = buttap(n);

[b,a]=zp2tf(z,p,k);

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(b,a,x);

N=56;n=0:N-1;

subplot(2,1,1);plot(n,x);subplot(2,1,2);plot(n,Y);

[bt,at]=lp2lp(b,a,Wn);

[bbs,abs]=bilinear(bt,at,fs);

[dB,mag,pha,grd,w]=freqz_m(bbs,abs);

Ripple=10^(-Rp/20);

Attn=10^(-Rs/20);

%»-ͼ

subplot(2,2,1);

plot(w/pi,mag);

title('·ùƵÏìÓ¦');

xlabel('¹éÒ»»¯ÆµÂÊ');

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

set(gca,'XTickMode','manual','XTick',[0,2*lfp/fs,2*lfs/fs,1]);

set(gca,'XTickMode','manual','XTick',[Attn,Ripple,1]);

grid;

subplot(2,2,2);

plot(w/pi,pha/pi);

title('ÏàƵÏìÓ¦');

xlabel('¹éÒ»»¯ÆµÂÊ');

ylabel('µ¥Î»:pi');

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

set(gca,'XTickMode','manual','XTick',[0,2*lfp/fs,2*lfs/fs,1]);

set(gca,'XTickMode','manual','XTick',[-1,-0.5,0,0.5,1]);

grid;

subplot(2,2,3);

plot(w/pi,dB,'red');

title('·ùƵÏìÓ¦:dB');

xlabel('¹éÒ»»¯ÆµÂÊ');

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

set(gca,'XTickMode','manual','XTick',[0,2*lfp/fs,2*lfs/fs,1]);

set(gca,'XTickMode','manual','XTick',[-60,-40,-20,5]);

grid;

subplot(2,2,4); plot(w/pi,grd,'red'); title('ȺʱÑÓ'); xlabel('¹éÒ»»¯ÆµÂÊ') axis([0,1,0,6]);

set(gca,'XtickMode','manual','Xtick',[0,2*lfp/fs,2*lfs/fs,1]); set(gca,'YtickMode','manual','Ytick',[0,1,2,4,6]); grid;

%²ÉÓÃÊý×ÖÓòƵÂʱ任·¨Éè¼ÆÊý×Ö´ø×èÂ˲¨Æ÷(using bilinear and Chebyshev-I)

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(b,a,x); N=56;n=0:N-1;figure(2)

subplot(2,1,1);plot(n,x);grid on ; subplot(2,1,2);plot(n,Y);grid on ;

图

双线性变换法设计的滤波器特性

0.0562

0.89131

0.20.40.60.8

1幅频响应

归一化频率0

0.5

1

-1-0.500.5

1相频响应

归一化频率单位:p i

-60

-40-20

0幅频响应:dB

归一化频率

0.2

0.6

1

01

24

6群时延

归一化频率

原始心电图信号和经过滤波后的心电图信号

0102030405060

-100

-50

50

10

20

30

40

50

60

-30-20-10010

20

实验心得：

脉冲响应不变法

代码：

%Âö³åÏìÓ¦²»±ä·¨

% digital filter specifications: lfp=100; lfs=300; ws=2*pi*lfs; wp=2*pi*lfp; Rp=1; Rs=25; fs=1000;

%Êý×ÖÖ¸±êµ½Ä£ÄâÖ¸±êµÄ±ä»¯¡ªË«ÏßÐԱ任 %¼ÆËã½×ÊýºÍ½ØֹƵÂÊ

[n,Wn]=buttord(wp,ws,Rp,Rs,'s'); %Éè¼ÆÄ£ÄâµÍͨÔ-ÐÍ [z,p,k]=buttap(n); [b,a]=zp2tf(z,p,k);

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-69,-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(b,a,x);

%ÓÉÄ£ÄâµÍͨÔ-Ð;-ƵÂʱ任»ñµÃÄ£ÄâµÍͨ¡¢¸ßͨ¡¢´øͨ»ò´ø×èÂ˲¨Æ÷

[bt,at]=lp2lp(b,a,Wn);

%ÀûÓÃÂö³åÏìÓ¦²»±ä·¨»òË«ÏßÐÔ²»±ä·¨£¬ÊµÏÖÄ£ÄâÂ˲¨Æ÷µ½Êý×ÖÂ˲¨Æ÷µÄÓ³Éä[bbs,abs]=impinvar(bt,at,fs);

%Â˲¨Æ÷ƵÂÊÏìÓ¦

[dB,mag,pha,grd,w]=freqz_m(bbs,abs);

Ripple=10^(-Rp/20);

Attn=10^(-Rs/20);

%»-ͼ

subplot(2,2,1);

plot(w/pi,mag) %Êý×ÖÂ˲¨Æ÷·ù¶ÈÏìÓ¦

title('·ùƵÏìÓ¦');

xlabel('¹éÒ»»¯ÆµÂÊ');

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

set(gca,'XTickMode','manual','XTick',[0,2*lfp/fs,2*lfs/fs,1]);

set(gca,'YTickMode','manual','YTick',[Attn,Ripple,1]);

grid;

subplot(2,2,2);

plot(w/pi,pha/pi); %Êý×ÖÂ˲¨Æ÷ÏàλÏìÓ¦

title('ÏàƵÏìÓ¦');

xlabel('¹éÒ»»¯ÆµÂÊ');

ylabel('µ¥Î»£ºpi');

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

set(gca,'XTickMode','manual','XTick',[0,2*lfp/fs,2*lfs/fs,1]);

set(gca,'YTickMode','manual','YTick',[-1,-0.5,0,0.5,1]);

grid;

subplot(2,2,3);

plot(w/pi,dB,'red');

%Êý×ÖÂ˲¨Æ÷·ù¶ÈÏìÓ¦£¨dB£©

title('·ùƵÏìÓ¦£ºdB');

xlabel('¹éÒ»»¯ÆµÂÊ');

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

set(gca,'XTickMode','manual','XTick',[0,2*lfp/fs,2*lfs/fs,1]);

set(gca,'YTickMode','manual','YTick',[-60,-40,-20,5]);

grid;

subplot(2,2,4);

plot(w/pi,grd,'red');

%Êý×ÖÂ˲¨Æ÷ÑÓʱ

title('ȺÑÓʱ');

xlabel('¹éÒ»»¯ÆµÂÊ');

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

set(gca,'XTickMode','manual','XTick',[0,2*lfp/fs,2*lfs/fs,1]);

set(gca,'YTickMode','manual','YTick',[0,1,2,4,6]);

grid;

%²ÉÓÃÊý×ÖÓòƵÂʱ任·¨Éè¼ÆÊý×Ö´ø×èÂ˲¨Æ÷£¨using bilinear and Chebyshev-1£©

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,-2,-4,0,0,0,-2,-2,0 ,0,-2,-2,-2,-2,0]

Y =filter(b,a,x);

N=56;n=0:N-1;figure(2)

subplot(2,1,1);plot(n,x);grid on; subplot(2,1,2);plot(n,Y);grid on;

子函数：

function[dB,mag,pha,grd,w]=freqz_m(b,a)

%computation of s-domain frequency response:modified version

%---------------------------------------------------------

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

% dB=Relative magintude in dB over [0 to 2pi]

% mag=Absolute magintude over [0 to 2pi]

% pha=Phase response in radians over [0 to 2pi]

% w=array of 500 frequency samples between [0 to 2pi]

% b=numerator polynomial coefficients of H(z)

% a=denominator polynomial coefficients of H(z)

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

%w ranging from 0 to 2pi radians per sample.

%uses 1000 sample points around the entire unit circle to calculate the %frequency response

H=(H(1:501))';

w=(w(1:501))';

mag=abs(H);

dB=20*log10((mag+eps)/max(mag));

pha=angle(H);

grd=grpdelay(b,a,w);

图：

0.2

0.6

1

0.0562

0.8913

1幅频响应

归一化频率0

0.2

0.6

1

-1-0.500.5

1

相频响应

归一化频率单位：p i

00.2

0.6

1-60

-40-20

5幅频响应：dB

归一化频率

00.2

0.6

1

124

6群延时

归一化频率

0102030405060

-100

-50

50

10

20

30

40

50

60

-30-20-10010

20

实验七

代码：

主程序：

% design a digital FIR lowpass filter with Hamming window

wp=0.3*pi;

ws=0.5*pi;

tr_width=ws-wp;

N=ceil(6.6*pi/tr_width)+1; % N=34

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

wc=(ws-wp)/2;

hd = ideal_lp(wc,N);

w_ham=(hamming(N))';

h=hd.*w_ham;

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

delta_w=2*pi/1000;

Rp=-(min(dB(1:1:wp/delta_w+1)));

As=-round(max(dB(ws/delta_w+1:1:501)));

%plots

subplot(2,2,1);

stem(n,hd);

title('ÀíÏëÂö³åÏìÓ¦');

axis([0,N-1,-0.09,0.5]);

xlabel('n');

set(gcr,'YTickMode','manual','YTick',[0,0.1,0.2,0.3]);

subplot(2,2,2);

stem(n,w_ham);

title('hamming´°');

axis([0,N-1,0,1.1]);

xlabel('n');

subplot(2,2,3);

stem(n,h);

title('ʵ¼ÊÂö³åÏìÓ¦');

axis([0,N-1,-0.09,0.5]);

xlabel('n');

set(gcr,'YTickMode','manual','YTick',[0,0.1,0.2,0.3]);

subplot(2,2,4);

plot(w/pi,dB);

title('·ùƵÏìÓ¦:dB');

grid;

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

xlabel('pi');

ylabel('dB');

set(gcr,'YTickMode','manual','YTick',[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0 .8,0.9,1.0]);

set(gcr,'YTickMode','manual','YTick',[-100,-90,-80,-70,-60,-50,-40,-3 0,-20,-10,0]);

子程序：

function[dB,mag,pha,grd,w]=freqz_m(b,a)

%computation of s-domain frequency response:modified version

%---------------------------------------------------------

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

% dB=Relative magintude in dB over [0 to 2pi]

% mag=Absolute magintude over [0 to 2pi]

% pha=Phase response in radians over [0 to 2pi]

% w=array of 500 frequency samples between [0 to 2pi]

% b=numerator polynomial coefficients of H(z)

% a=denominator polynomial coefficients of H(z)

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

%w ranging from 0 to 2pi radians per sample.

%uses 1000 sample points around the entire unit circle to calculate the %frequency response

H=(H(1:501))';

w=(w(1:501))';

mag=abs(H);

dB=20*log10((mag+eps)/max(mag));

pha=angle(H);

grd=grpdelay(b,a,w);

函数ideal_lp;

function hd=ideal_lp(wc,N);

alpha=(N-1)/2;

n=0:1:N-1;

m=n-alpha+eps;

hd=sin(wc*m)./(pi*m);

汉明窗法设计滤波器的

10

20

30

理想脉冲响应

n

010

2030

hamming 窗

n

10

20

30

实际脉冲响应

n 0

0.51

-100

-50

幅频响应:dB

pi

d B