最新数字信号处理实验报告x

前 言 错 误 !未定义书签。

TOC \o "1-5" \h \z \o "Current Document" 实验一 MATLAB 简介 2

\o "Current Document" 实验二 用 FFT 实现信号的谱分析 5

\o "Current Document" 实验三 IIR 数字巴特沃思滤波器的设计 16

\o "Current Document" 实验四 FIR 数字滤波器的设计 19

实验一 MATLAB简介

实验目的

.熟悉MATLAB软件的使用方法；

. MATLAB 的绘图功能；

.用MATLAB语句实现信号的描述及变换。

实验原理

.在MATLAB下编辑和运行程序

在MATLAB 中，对于简单问题可以在命令窗( comma nd win dows )直接输入命令， 得到结果；对于比较复杂的问题则可以将多个命令放在一个脚本文件中， 这个脚本文件是以

m为扩展名的，所以称之为 M文件。用M文件进行程序的编辑和运行步骤如下：

(1 )打开MATLAB，进入其基本界面；

在菜单栏的File项中选择新建一个 M文件；

在M文件编辑窗口编写程序；

(4 )完成之后，可以在编辑窗口利用 Debug工具调试运行程序，在命令窗口查看输

出结果；也可以将此文件保存在某个目录中，在 MATLAB 的基本窗口中的 File项中选择

Run The Script ，然后选择你所要运行的脚本文件及其路径，即可得出结果；也可以将此

文件保存在当前目录中，在 MATLAB命令窗口， “ >> ”提示符后直接输入文件名。

. MATLAB 的绘图功能

plot(x,y)基本绘图函数，绘制x

plot(x,y)

基本绘图函数，绘制

x和y之间的坐标图。

figure(n ) 开设一个图形窗口 n

subplot(m,n,N) 分割图形窗口的 MATLAB 函数，用于在一个窗口中显示多个图形, 将图形窗口分为 m行n列，在第N个窗口内绘制图形。

axis([aO,bO,a1,b1]) title(

axis([aO,bO,a1,b1]) title('') xlabel ( ‘ ‘)

ylabel ( ‘ ‘)

给图形加题注

给x轴加标注

给y轴加标注

grid 给图形加网格线

?信号描述及变换

信号描述及变换包括连续时间信号和离散时间信号内容，详细内容请见课本第 第2章。

实验内容

1 .试用

MATLAB绘制出下列信号的波形：

(Sig nal 1.6)

(1)

X1(t)

1.5t

e；

(2)

X2(t)

3sin(0.5 t)

(3)

X3(t)

0.5 0.5sgn(t)；

X4(t) u(t) u(t 1) 2u(t 2)；

X5(t) :[u(t) u(t 4)]

2

【程序代码】

clear all ;close all ;clc; syms t;

x1=exp(-1.5*t)

x2=3*s in (0.5*pi*t) x3=0.5+0.5*sym(( 'sign(t)'))

x4=sym( 'heaviside(t)' )+sym( 'heaviside(t-1)' )-sym( '2*heaviside(t-2)' ) x5=.5*t*(sym( 'heaviside(t)' )-sym( 'heaviside(t-4)'))

subplot(2,3,1);

ezplot(x1);

axis([-6 3 -500 7000]);

title( 'x1(t)=exp(-1.5t)');

grid on

subplot(2,3,2);

ezplot(x2);

title( 'x2(t)=3sin(0.5|Dt)' );

grid on

subplot(2,3,3);

fplot( 'sign(t)/2 + 1/2' ，[-10 10],1e-8);

ezplot(x3,[-10 10]);

axis([-10 10 -.2 1.2]);

TOC \o "1-5" \h \z xlabel( 't' );title( 'x3(t)=0.5+0.5sgn(t)' );

grid on

subplot(2,3,4);

ezplot(x4,[-1 3]);

title( 'x4(t)=u(t)+u(t-1)-2u(t-2)' );

grid on

subplot(2,3,5);

ezplot(x5,[-2 6]);

title( 'x5(t)=0.5t[u(t)-u(t-4)]' );

grid on

subplot(2,3,6);

axis off

2 .已知连续时间信号

2 .

已知连续时间信号(Signal 1.7)

X1 (t)

(4 t)[u(t) u(t 4)] , X2(t)

试用

MATLAB

绘制出下列信号的波形：

(1)

X4(t)

xNt/2)；

(2)

X5(t)

X4(t 2)；

(3)

xe(t)

X2( t)；

(4)

X7(t)

X2(t) X6(t)；

2t

e u(t), X3(t)

sin(2 t)

(5) x8(t) x7(t) x3(t) 。

【程序代码】

clear all ;close all ;clc;

figure(2)

syms t;

x1=(4-t)*(sym( 'heaviside(t)' )-sym( 'heaviside(t-4)' )) x2=exp(-2*t)*sym( 'heaviside(t)' )

x3=sin(2*pi*t)

x4=subs(x1,t,t/2)

x5=subs(x4,t,t-2)

x6=subs(x2,t,-t)

x7=x2+x6

x8=x7*x3

subplot(2,3,1)

text(0,0.9, 'x1(t)=(4-t)[u(t)-u(t-4)] ' );

text(0,0.7, 'x2(t)=exp(-2t)u(t)' );

text(0,0.5, 'x3(t)=sin(2|Dt)' );

axis off ;box off;

subplot(2,3,2);

ezplot(x5,[-3,10]);grid

title( 'x4(t)=x1(t/2)' );

subplot(2,3,3); ezplot(x5,[-1,12]);grid title( 'x5(t)=x4(t-2)' );

subplot(2,3,4);

ezplot(x6,[-10,5]);

axis([-5 5 -0.1 1.1]);

grid

title( ' x6(t)=x2(-t)' );

subplot(2,3,5);

ezplot(x7,[-5 5]);

axis([-5 5 -0.1 1.1]);

grid;

title( 'x7(t)=x2(t)+x6(t)' );

subplot(2,3,6);

tv8=-2.5:0.05:2.5;

xv8=subs(x8,tv8);

plot(tv8,xv8)

axis([-2.5 2.5 -1 1]);

grid

xlabel( 't' );title( 'x8(t)=x7(t)x3(t)' );

clear tv8 xv8

3 ．列出单位冲激信号、单位阶跃信号、正弦信号的 MATLAB 表达式，并绘出信号波 形。

【程序代码】

clear all ;close all ;clc

syms t;

x1=sym( 'dirac(t)');

x2=sym( 'Heaviside(t)');

x3=si n(t);

tn=[-6.3:0.163];

xn 1=subs(x1,t,t n);

xn 2=subs(x2,t,t n);

xn 3=subs(x3,t,t n);

plot(tn,xn1, 'k' ,tn,xn2, 'r' ,tn,xn3, 'm');

grid

xlabel( 't' );ylabel( 'x(t)');

legend( 'dirac' , 'Heaviside' ,'sin')

hold on

plot(tn,xn1, 'k.' ,tn,xn2, 'r.' ,tn,xn3, 'm.')

实验二 用FFT实现信号的谱分析

实验目的

1?了解FFT在信号谱分析中的作用；

2 ?了解谱分析的一般步骤和方法。

实验原理

关于信号谱分析的步骤和方法参见教材第 3章相关内容。为了解信号的特点，了解信

号频谱分布情况，应该对信号进行谱分析，计算出信号的幅度谱、 相位谱和功率谱。信号的

谱分析可以用FFT实现，讨论如下：

?谱分析中的参数选择；

TOC \o "1-5" \h \z A 若已知信号的最高频率 fc，为防止混叠，选定采样频率 fs :

fs 2 fc (1)

B根据实际需要，选定频率分辨 f，一但选定后，即可确定 FFT所需的点数N

N fs/ f (2)

我们希望 f越小越好，但 f越小，N越大，计算量、存储量也随之增大。一般取

N为2的整次幕，以便用 FFT计算，若已给定 N,可用补零方法便 N为2的整次幕。

C fs和N确定后，即可确定所需相应模拟信号 x(t)的长度

T N/fs NTs (3)

分辨率f反比于T,而不是N,在给定的T的情况下，靠减小 Ts来增加N是不能提 高分辨率的，因为T NTs为常数

?谱分析步骤；

A数据准备

x(n) xa(t)t nT Xa(nT) (4)

B使用FFT计算信号的频谱

N 1

N 1

kn

X(k) x(n)WN

n 0

TOC \o "1-5" \h \z X(k) Xr(k) jXi(k) ⑹

C由频谱计算幅度谱 X(k)、相位谱k和功率谱G(k)

X(k) Jx；(k) X：(k) ⑺

k arctan Xi(k) (8)

Xr(k)

G(k) X(k)2 Xr2(k) Xi2(k) (9)

3 ?实验中用到的一些基本函数简介

y=fft(x ,n)

; 计算n点的FFT。

abs(x)

; 取绝对值。

an gle(z)

; 取相角。

[Pxx, f]= periodogram (xn, nfft, fs, window) ; %周期图谱估计

[Pxx, f]=pwelch (xn, nfft, fs, window, noverlap) ; %平均周期图法

Pxx=psd (xn)

;功率谱密度

实验内容

1.已知序列 x(n)=2sin(0.48 n n)+cos(0.52 0< n<100,试绘制 x(n)及它的频谱

图。

若 x(n)=sin(0.56 tT )+2cos(0.25 nn ),结果又如何？

【程序代码】

clear all;close all;clc

N=100;

n=0:N-1;

xn=2*si n(0.48*pi* n)+cos(0.52*pi* n);

XK=fft(x n, N);

magXK=abs(XK);

phaXK=a ngle(XK);

subplot(1,2,1)

plot( n,xn)

xlabel( 'n');ylabel( 'x(n)');

title( 'x(n)N=100');

subplot(1,2,2)

k=O:le ngth(magXK)-1;

stem(k,magXK, '.');

xlabel( 'k' );ylabel( '|X(K)|');

title( 'X(K)N=100');

2.对下面信号进行频谱分析，求幅度谱 X(k)和相位谱(k)。

(1)花⑴⑵ X2(t)at , a 0.8, 0 t 4ms, f

(1)花⑴

⑵ X2(t)

sin t/t , T 0.125s , N 16

【程序代码】

clear all ;close all ;clc fs1=5000;

dt1=1/fs1;

N1=0.004/dt1;

n1=0:N1-1;

xn1= 0.8.A( n1*dt1);

Xk1=fft(x n1,N1);

mag1=2*abs(Xk1)/N1;

pha1=a ngle(Xk1);

f1=n 1*fs1/N1;

subplot(221)

stem(f1,mag1, 'fill');

title( 'Magnitude 1') xlabel( 'f');ylabel( '|Xk(1)|');

subplot(222)

stem(f1,pha1, 'fill');

title(卩hase 1')

xlabel( 'f');ylabel( '|Phk(1)|');

fs2=128;

dt2=1/fs2;

N2=18;

n2=0:N2-1;

xn 2=si n(n 2*dt2)./( n2*dt2+eps);

Xk2=fft(x n2,N2);

mag2=2*abs(Xk2)/N2;

pha2=a ngle(Xk2);

f2=n 2*fs2/N2;

subplot(223)

stem(f2,mag2, 'fill');

title( 'Magnitude 2')

xlabel( 'f');ylabel( '|Xk(2)|');

subplot(224)

stem(f2,pha2, 'fill');

title(卩hase 2')

xlabel( 'f');ylabel( '|Phk(2)|')

3.给定信号 x(t) sin(2 f1t) 2sin(2 f2t) , h 15Hz , h 18Hz，现在对 x(t)采

样，采样点数N 16，采样频率fs=50Hz，设采样序列为x(n)，编写程序计算x(n)的频

谱，并绘图；改变采样频率，得到序列 x,n)，计算X,n)的频谱，并绘图；增大采样点数，

得到序列x2(n)，计算x2(n)的频谱，并绘图；采样点数 N=64，采样频率fs =300Hz，在

采样点后补零得到新序列 X3(n)，计算X3(n)的频谱，并绘图。(Signal 3.18)

【程序代码】

clear all ;close all ; clc

N=16;

n=0:N-1;

f1=15;f2=18;

fs=50;

dt=1/fs;

xn=si n(2*pi*f1* n*dt)+2*si n(2*pi*f2* n*dt);

Xk=fft(x n, N);

mag=2*abs(Xk)/N;

pha=a ngle(Xk);

f=fs/N*n;

figure(1); subplot(121) plot(f,mag);grid title( 'Magnitude' );xlabel( 'f' );ylabel( '|Xk|' );

subplot(122) plot(f,pha);grid title( 'Phase' );xlabel( 'f' );ylabel( '? § Xk);' N=16;

n=0:N-1; f1=15;f2=18;

fs=36; dt=1/fs;

xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt); Xk=fft(xn,N);

mag=2*abs(Xk)/N; pha=angle(Xk); f=fs/N*n;

figure(2); subplot(221) plot(f,mag, 'r' );grid title( 'Magnitude(36Hz Sampling requency)' );

xlabel( 'f' );ylabel( '|Xk|' );

subplot(222)

plot(f,pha, 'r' );grid

title( 'Phase(36Hz Sampling requency)' );

xlabel( 'f' );ylabel( '? § Xk);' fs=100;

dt=1/fs; xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);

Xk=fft(xn,N); mag=2*abs(Xk)/N; pha=angle(Xk); f=fs/N*n;

subplot(223) plot(f,mag, 'r' );grid title( 'Magnitude(100Hz Sampling Frequency)' ); xlabel( 'f' );ylabel( '|Xk|' );

subplot(224)

plot(f,pha, 'r' );grid

title( 'Phase(100Hz Sampling Frequency)' ); xlabel( 'f' );ylabel( '? § Xk);'

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

f1=15;f2=18;

fs=50;

dt=1/fs;

xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);

Xk=fft(xn,N); mag=2*abs(Xk)/N; pha=angle(Xk); f=fs/N*n; figure(3);

subplot(221)

plot(f,mag, 'm' );grid

title( 'Magnitude(64 Sampling Points)' );xlabel( 'f' );ylabel( '|Xk|' ); subplot(222) plot(f,pha, 'm' );grid

title( 'Phase(64 Sampling Points)' );xlabel( 'f' );ylabel( '? § Xk);' N=512;

n=0:N-1;

dt=1/fs;

xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);

Xk=fft(xn,N); mag=2*abs(Xk)/N; pha=angle(Xk); f=fs/N*n;

subplot(223) plot(f,mag, 'm' );grid title( 'Magnitude(512 Sampling Points)' );xlabel( 'f' );ylabel( '|Xk|' );

subplot(224) plot(f,pha, 'm' );grid

title( 'Phase(512 Sampling Points)' );xlabel( 'f' );ylabel( '? § Xk);'

N=64;

n=0:N-1;

f1=15;f2=18;

fs=300;

dt=1/fs;

xn0=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt); xn=[xn0,zeros(1,64)];

Xk=fft(xn,N); mag=2*abs(Xk)/N; pha=angle(Xk); f=fs/N*n; figure(4);

subplot(221) plot(f,mag, 'g' );grid title( 'Magnitude(filled with 64 zeros)' );xlabel( 'f' );ylabel( '|Xk|' );

subplot(222)

plot(f,pha, 'g' );grid

title( 'Phase (filled with 64 zeros)' );xlabel( 'f');ylabel( '? § Xk);' xn=[xn0,zeros(1,448)];

Xk=fft(xn,N); mag=2*abs(Xk)/N;

pha=angle(Xk);

f=fs/N*n;

subplot(223)

plot(f,mag, 'g' );grid

title( 'Magnitude(filled with 448 zeros)' );xlabel( 'f' );ylabel( '|Xk|' );

4. 试求下列差分方程所描述的输出序列

4. 试求下列差分方程所描述的输出序列

x( n) 的功率谱并作图。

x(n) 0.81x(n

x(n) 0.81x(n 2) w(n) w(n 1)

x(n) w(n) w(n 2)

x(n) 0.81x(n 2) w(n)

式中， w( n )

式中， w( n )是方差为 w2 (例如，

w2 =1/12)

的白噪声。

【程序代码】

clear all ;close all ;clc nfft=512;

window=hanning(256); noverlap=128;

w=randn(1,1000);

for i=1:1000

if i<=2 x(i)=randn(1,1); else x(i)=-0.81*x(i-2)+w(i)-w(i-1); end end

'mean' );[Pxx f]=psd(x,nfft,2,window,noverlap, figure(1)

'mean' );

xlabel( 'f' );ylabel( 'Power Spectrum,dB' ); title( '(a): x(n)=-0.81x(n-2)+w(n)-w(n-1)' ) w=randn(1,1000);

for i=1:1000 if i<=2 x(i)=randn(1,1);

else

x(i)=w(i)-w(i-2);

end

end

[Pxx f]=psd(x,nfft,2,window,noverlap, 'mean' );

figure(2)

plot(f,10*log10(Pxx), 'm' );grid

xlabel( 'f' );ylabel( 'Power Spectrum,dB' );

title( '(b): x(n)=w(n)-w(n-2)' )

w=randn(1,1000);

for i=1:1000

if i<=2

x(i)=randn(1,1);

else

x(i)=-0.81*x(i-2)+w(2);

end

end

[Pxx f]=psd(x,nfft,2,window,noverlap);

figure(3)

plot(f,10*log10(Pxx), 'r' );grid

xlabel( 'f' );ylabel( 'Power Spectrum,dB' );

title( '(c): x(n)=-0.81x(n-2)+w(n)' )

5. 一序列 x(n) 是由两个频率相距为 f 的模拟信号采样得来的，即

x( n ) si n2 (0.135) n cos 2 (0.135 f )n n=0,1,…,15

已知序列长度 N=16 ，试采用周期图法，应用 DFT 分别计算当 f =0.06 及 f =0.01 时的

功率谱估计，并通过作图说明从功率谱估计的分布是否能分辨出这两个正弦信号的真实频

谱？若 N=64 又有什么变化？

【程序代码】

clear all ;close all ;clc

nfft=512;window=hanning(256);noverlap=128;df=[0.06 0.01];

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

x1=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(1))*n);

x2=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(2))*n);

[Pxx1 f1]=psd(x1,nfft,1,window,noverlap, 'mean' );

[Pxx2 f2]=psd(x2,nfft,1,window,noverlap, 'mean' );

figure(1)

plot(f1,Pxx1);grid

xlabel( 'f');ylabel( 'Power Spectrum' );

title( 'sin2(0.135)n+cos2(0.135+0.06)n N=16' figure(2) plot(f2,Pxx2);grid

xlabel( 'f');ylabel( 'Power Spectrum ' );

title( '(sin2(0.135)n+cos2(0.135+0.01)n N=16'

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

x1=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(1))*n);

x2=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(2))*n);

'mean' );'mean' );

'mean' );

'mean' );

[Pxx2 f2]=psd(x2,nfft,1,window,noverlap, figure(3)

plot(f1,Pxx1, 'm' );grid

xlabel( 'f');ylabel( 'Power Spectrum' );

title( 'sin2(0.135)n+cos2(0.135+0.06)n N=64'

figure(4);plot(f2,Pxx2, 'm' );grid xlabel( 'f');ylabel( 'Power Spectrum' );

title( '(sin2(0.135)n+cos2(0.135+0.01)n N=64'

6. 用

6. 用 MATLAB 产生 256 点白噪声序列，

应用 Welch 法估计其功率谱， 每段长 64 点，

重叠 32

重叠 32 点，输出平均后的功率谱图以及对

256 点一次求周期图的功率谱图。

程序代码】

N=256; x=rand(1,N);

window=hanning(64); noverlap=32;nfft=128;

[Pxx1 f1]=periodogram(x); figure(1) plot(f1/pi,10*log10(Pxx1));grid xlabel( 'f');ylabel( 'Power Spectrum,dB'title( '256 卩? o ?' o ? u ? u);[1?1|?

[Pxx1 f1]=periodogram(x); figure(1) plot(f1/pi,10*log10(Pxx1));grid xlabel( 'f');ylabel( 'Power Spectrum,dB'

title( '256 卩? o ?' o ? u ? u

);

[1

?1|?

e )? X] ?'

[Pxx2 f2]=psd(x,nfft,2,window,noverlap, figure(2)

plot(f2,10*log10(Pxx2), 'm' );grid xlabel( 'f');ylabel( 'Power Spectrum,dB' title( 'Welch ? ？？？？ o ? u ? u

'line'

);

);

[1

?1|?

7. 离散信号序列 x(n) 2sin(2

fn/ fs) ，

长度 N=512 ， f s 32Hz ，令 f 取值分别

为 16Hz 和 15Hz ，计算序列的功率谱，

比较谱图的差别。试采用不同的窗函数，再比较谱

图的变化。

程序代码】

clear all ;close all ;clc

N=512;

n=0:N-1;

nfft=512;

noverlap=128; window=hanning(256);

fs=32;

f=16;

x=2*sin(2*pi*f*n/fs);

[Pxx f]=psd(x,nfft,fs,window,noverlap, subplot(321)

'mean' );

plot(f,Pxx);grid

xlabel( 'f' );ylabel( 'Power Spectrum' ) title( 'hanning window f=16' ) f=15;

x=2*sin(2*pi*f*n/fs);

[Pxx f]=psd(x,nfft,fs,window,noverlap, subplot(322)

;

'mean' );

plot(f,Pxx);grid

xlabel( 'f' );ylabel( 'Power Spectrum' ) title( 'hanning window f=15' ) window=hanning(256);

f=16;

x=2*sin(2*pi*f*n/fs);

[Pxx f]=psd(x,nfft,fs,window,noverlap, subplot(323)

;

'mean' );

plot(f,Pxx);grid

xlabel( 'f' );ylabel( 'Power Spectrum' ) title( ' boxcar window f=16' ) f=15;

x=2*sin(2*pi*f*n/fs);

[Pxx f]=psd(x,nfft,fs,window,noverlap, subplot(324)

;

'mean' );

plot(f,Pxx);grid

xlabel( 'f' );ylabel( 'Power Spectrum' ) title( ' boxcar window f=15' ) window=kaiser(256);

f=16;

x=2*sin(2*pi*f*n/fs);

[Pxx f]=psd(x,nfft,fs,window,noverlap, subplot(325)

;

'mean' );

plot(f,Pxx);grid

xlabel( 'f' );ylabel( 'Power Spectrum' title( ' kaiser window f=16' )

f=15;

x=2*sin(2*pi*f*n/fs);

[Pxx f]=psd(x,nfft,fs,window,noverlap, 'mean' ); subplot(326)

plot(f,Pxx);grid

xlabel( 'f' );ylabel( 'Power Spectrum' );

title( ' kaiser window f=15' )

8.已知一个被白噪声r(t)污染的信号x(t),

f1 =25Hz ，x(t) 2sin(2 f1t) 0.5sin(2 f2t) 0.5sin(2 f3t) r (t ) ，其中， f2 =75Hz ， f3 =150Hz

f1 =25Hz ，

clear;close;clc N=512;n=0:N-1;

f1=25;f2=75;f3=150; Fs=500;Ts=1/Fs;

noverlap=80;nfft=256;rn=randn(1,N); window=hanning(128);

xn=2*sin(2*pi*f1*n*Ts)+0.5*sin(2*pi*f2*n*Ts)+0.5*sin(2*pi*f3*n*Ts)+rn; [Pxx f]=psd(xn,nfft,Fs,window,noverlap, 'mean' );

plot(f,10*log(Pxx));grid

xlabel( 'f/Hz' );ylabel( 'Power Spectrum,dB' );

title( 'Welch ?…1|? e ? X 1 a f1=25 Hz,f2=75Hz,f3=150Hz' )

figure

plot(f,Pxx, 'm' );grid

xlabel( 'f/Hz' );ylabel( 'Power Spectrum' );

title( 'Welch ?…1|? e ? X 1 a f1=25Hz,f2= 75Hz,f3=150Hz')

实验三 IIR 数字巴特沃思滤波器的设计

实验目的

1． 掌握用模拟滤波器原型设计 IIR 滤波器的基本方法；

2． 掌握数字巴特沃思滤波器的设计方法与步骤；

3． 理解系统频率响应的概念，学习编写计算系统频率响应的方法。

实验原理

1 ．数字巴特沃思滤波器设计的详细内容见教材第 4 章，现将设计步骤归纳如下： A 根据给定的频带指标，由双线性变换的频率关系，确定相应的模拟滤波器原型频带 指标；

B 利用原型低通滤波器，选择合适的参数，设计出符合指标的模拟低通滤波器；

C 利用双线性变换， 将所获得的模拟滤波器的 s 域表示转换为相应数字滤波器的 z 域 表示，即它的系统函数，再利用 IIR 滤波器设计方案具体实现该滤波器。

2. 实验中用到的一些基本函数见教材第 4.5.1 节(与 IIR 数字滤波器设计相关 MATLAB 函数)。

实验内容

1 ． 希望设计一个巴特沃斯低通数字滤波器，其 3dB 带宽为 0.2 ，阻带边缘频率为

0.5 ，阻带衰减大于 30dB 。给定采样间隔 Ts 10 s 。

用双线性变换法设计该低通数字滤波器。给出它的 H (z) 及对数幅频响应。

【程序代码】

clear all ;close all ;clc

Ts=10e-6;

fs=1/Ts;

wp=0.2*pi;ws=0.5*pi;

rp=3;rs=30;

[n wn]=buttord(wp/pi,ws/pi,rp,rs);

[numz denz]=butter(n,wp/pi);

Gz=tf(numz,denz,Ts)

[h w]=freqz(numz,denz,128,fs);

plot(w,abs(h), 'r' );

xlabel( 'f' );ylabel( 'H' );

title( 'LowPass Filter' );

grid on

2 ． 给定待设计的数字高通和带通滤波器的技术指标如下：

(1) HP ： f p 400Hz ， fs 300Hz ， Fs 1000Hz ， p 3dB ， s 35dB 。

(2) BP ： fsl 200Hz ， f1 300Hz ， f2 400Hz ， fsh 500Hz ，Fs 2000Hz ， p 3dB ， s 40dB 。

试用双线性变换分别设计满足上述要求的巴特沃斯滤波器，给出其系统函数、对数幅 频及相频曲线。

【程序代码】

clear all ;close all ;clc

Fs=1000;

Ts=1/Fs; fp=400;fs=300;

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

rp=3;rs=35;

[n wn]=buttord(wp/pi,ws/pi,rp,rs);

[numz denz]=butter(n,wp/pi, 'high' );

Gz_HP=tf(numz,denz,Ts)

[h w]=freqz(numz,denz,128,Fs);

figure(1)

subplot(2,2,[1 2]) plot(w,20*log10(abs(h)), 'r' );grid xlabel( 'f' );ylabel( '20log(H)' );box off title( 'High Pass Filter' );

subplot(2,2,3) plot(w,abs(h), 'r' );

xlabel( 'f' );ylabel( 'H' );box off subplot(224)

plot(w,angle(h), 'm' );grid

xlabel( 'f' );ylabel( 'phase' );

grid on ;box off

clear all;

Fs=2000;

Ts=1/Fs;

fs1=200;f1=300;

f2=400;fsh=500;

wp=2*pi*fs1*Ts;ws=2*pi*fsh*Ts; w1=2*pi*f1*Ts;w2=2*pi*f2*Ts;

rp=3;rs=40;

[n wn]=buttord(wp/pi,ws/pi,rp,rs);

[numz denz]=butter(n,[w1/pi w2/pi]);

Gz_BP=tf(numz,denz,Ts)

[h w]=freqz(numz,denz,128,Fs); figure(2)

subplot(2,2,[1 2]); plot(w,20*log10(abs(h)), 'r' );

xlabel( 'f' );ylabel( '20log(H)' ); grid on ;box off title( 'Band Pass Filter' ); subplot(2,2,3);

plot(w,abs(h), 'r' );grid xlabel( 'f' );ylabel( 'H' ); grid on ;box off subplot(2,2,4);

plot(w,angle(h), 'm' ); xlabel( 'f' );ylabel( 'phase' ); grid on ;box off clear all

实验四 FIR数字滤波器的设计

实验目的

?掌握FIR数字滤波器的设计方法与步骤；

?理解系统频率响应的概念，学习编写计算系统频率响应的方法。

实验原理

FIR数字滤波器设计的详细内容见教材第 4章。

窗口法：窗口法设计 FIR数字滤波器的步骤：

A给出希望的滤波器频率响应函数 H d (ej )；

B根据允许的过渡带宽度及阻带衰减确定所采用的窗函数和 N值；

C 做Hd(ej )的逆傅里叶变换得 hd(n)；

D对hd(n)加窗处理得到有限长序列 h(n) hd(n) w(n)

E对h(n)做傅里叶变换得到频率响应 H(ej )，用H(ej )作为H d (ej )的逼近，并

用给定的技术指标来检验。

切比雪夫一致逼近计算机辅助设计方法：用 Remez算法实现FIR滤波器的等波纹

逼近。详细内容教材第 443节。

实验中用到的一些基本函数见教材第第 4.5.2 节(与FIR数字滤波器设计相关

MATLAB 函数)。

实验内容

1?请选择合适的窗函数及 N来设计一个线性相位低通滤波器

Hd

Hd(ej )

0,

要求其最小阻带衰减为—45dB，过渡带宽为8/51

(1)已知c 0.5 ，求出h(n)并画出20lgH(ej )曲线。

(2)保留原有轨迹，画出用满足所给条件的其他几种窗函数设计出的 20lg H(ej )曲线。

【程序代码】

clear all ;close all ;clc rs=45; w0=8*pi/51; w0_ham=3.3;

N=2*pi*w0_ham/w0;

M=ceil(N-1)-1; wc=0.5*pi; hn=fir1(M,wc/pi,hamming(M+1)) subplot(211)

stem(0:M,hn, 'filled' ); xlabel( 'n' );ylabel( 'h(n)' ); title( 'Hamming Window' ); [h w]=freqz(hn,1,256);

figure(1) subplot(212) plot(w,20*log10(abs(h)));

xlabel( 'w');ylabel( '20lg|H(eAjw)|' );

title( 'Freqency Magnitude(dB)' );

grid on figure(2);

plot(w,20*log10(abs(h)), '-.' ); clear all ;

w0=8*pi/51;

w0_kai=5.0;

N=2*pi*w0_kai/w0;

M=ceil(N-1)-1; wc=0.5*pi;

hn=fir1(M,wc/pi,kaiser(M+1));

[h w]=freqz(hn,1,256);

hold on plot(w,20*log10(abs(h)), 'g'); w0_han=5.0;

N=2*pi*w0_han/w0;

M=ceil(N-1)-1; hn=fir1(M,wc/pi,hanning(M+1));

[h w]=freqz(hn,1,256); plot(w,20*log10(abs(h)), 'm:' ); w0_box=0.9;

N=2*pi*w0_box/w0;

M=ceil(N-1)-1; hn=fir1(M,wc/pi,boxcar(M+1));

[h w]=freqz(hn,1,256); plot(w,20*log10(abs(h)), 'r--' );

xlabel( 'w');ylabel( '20lg|H(eAjw)|' );

title( 'Freqency Magnitude(dB)' );

legend( 'hamming' ，‘kaiser' ,'hanning' ，‘boxcar')

grid on

clear all

2 给定一理想低通FIR滤波器的频率特性

Ha(ej )

/4

希望用窗函数（矩形窗和哈明窗）法设计该滤波器，要求具有线性相位。假定滤波器系数的 长度为29，即M/2 = 14。

【程序代码】

clear all ;close all ;clc

N=29;

M=N-1;

wc=0.25*pi;

hn_1=fir1(M,wc/pi,hammi ng(M+1))

figure(1)

subplot(211) stem(0:M,hn_1, 'filled'); xlabel( 'n' );ylabel( 'h(n)');

title( 'Hamming Window' );

[h_1 w_1]=freqz(hn_1,1,256);

subplot(212) plot(w_1,abs(h_1)); xlabel( 'w' );ylabel( 'H');

grid on

hn_2=fir1(M,wc/pi,boxcar(M+1))

figure(2)

subplot(211)

stem(0:M,hn_2, 'filled' ,'m'); xlabel( 'n' );ylabel( 'h(n)'); title( 'Boxcar Window' );

[h_2 w_2]=freqz(hn_2,1,256);

subplot(212) plot(w_2,abs(h_2), 'm'); xlabel( 'w' );ylabel( 'H');

grid on

figure(3)

plot(w_1,abs(h_1),w_2,abs(h_2), 'r-.');

xlabel( 'w' );ylabel( 'H' );grid on legend( 'hamming' , 'boxcar');

title( 'contrast'

title( 'contrast'