数字信号处理第七章作业

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

1,4jHde

0,4现希望用窗函数法设计该滤波器，要求具有线性相位。假定滤波器系数的长度为29

点，即M/2=14。

试计算并打印滤波器的系数、幅频响应及相频响应。滤波器系数的计算先用手算，然后再调用有关的MATLAB文件来计算。

解：线性相位理想低通滤波器的频率特性：

jwe,4jHde

0,4

1hdn2ejejn1d2ejnde2jj1ejn2jn441enj4n4n

sinn4,N114n2（1） 矩形窗：

1,0n28 窗函数:wrn0,其它sinn144,0n28 h1nhdnwrnn140,其它n

系数算法一：

clc clear

for n=0:1:28

x=sin(0.25*pi*(n-14))/(pi*(n-14)) n=n+1 end

系数算法二： clc clear n=0:1:28

b=sin(0.25*pi*(n-14))./(pi*(n-14))

求得滤波器系数为：h(0)——h(28)=

-0.0227 -0.0173 0.0000 0.0205 0.0318 0.0250 -0.0000 -0.0322 -0.0531 -0.0450 0.0000 0.0750 0.1592 0.2251 0.2500 0.2251 0.1592 0.0750 0.0000 -0.0450 -0.0531 -0.0322 -0.0000 0.0250 0.0318 0.0205 0.0000 -0.0173 -0.0227 （2） 汉明（Hamming）窗：

2n0.0.46cos,n0,1,,28

窗函数:whnN0,其它nsinn1440.0.46cos2n,0n28 h2nhdnwhnn14N0,其它n

Matlab程序：

clc clear n=0:1:28

b=sin(0.25*pi*(n-14))./(pi*(n-14)).*(0.-0.46*cos(2*pi*n/29))

求得滤波器系数为：h(0)——h(28)=

-0.0018 -0.0016 0.0000 0.0036 0.0077 0.0081 -0.0000 -0.0166 -0.0326 -0.0320 0.0000 0.0656 0.1487 0.2197 0.2493 0.2245 0.1553 0.0701 0.0000 -0.0359 -0.0377 -0.0198 -0.0000 0.0104 0.0103 0.0050 0.0000 -0.0021 -0.0021 （3） 布莱克曼（Blackman）窗：

2n4n0.420.5cos0.08cos,n0,1,,28

窗函数:wbnNN0,其它nh3nhdnwbnsinn1440.420.5cos2n0.08cos4n,0n28 n14NN0,其它n

Matlab程序：

clc clear n=0:1:28

b=sin(0.25*pi*(n-14))./(pi*(n-14)).*(0.42-0.5*cos(2*pi*n/29)+0.08*cos(4*pi*n/29))

求得滤波器系数为：h(0)——h(28)=

0.0000 -0.0001 0.0000 0.0009 0.0027 0.0035 -0.0000 -0.0101 -0.0226 -0.0246 0.0000 0.0591 0.1410 0.2155 0.2488 0.2240 0.1524 0.0665 0.0000 -0.0302 -0.0290 -0.0137 -0.0000 0.00 0.0045 0.0017 0.0000 -0.0003 -0.0001 （4） Matlab代码及运行结果：

clc; clear; N=28;

br=0.25*sinc((-14:1:14)/4);%计算滤波器系数

bb1=fir1(N,0.25,boxcar(N+1));%用fir1计算滤波器系数 [H1,wr]=freqz(br,1,100,2);%频域变换 [Hb1,wr]=freqz(bb1,1,100,2); subplot(331);

plot(wr,abs(Hb1),'-r',wr,abs(H1),'-b');grid;%画出幅频特性 title('29点矩形窗') subplot(334);

plot(wr,unwrap(angle(Hb1)),'-r',wr,unwrap(angle(H1)),'-b');grid;%画出相频特性 subplot(337);

plot(wr,20*log10(abs(Hb1)/sum(Hb1)),'-r',wr,20*log10(abs(H1)/sum(H1)),'-b');grid;

title('对数幅频特性')

bh=br.*(0.-0.46*cos(2*pi*(0:1:28)/29)); bb2=fir1(N,0.25,hamming(N+1)); [H2,wh]=freqz(bh,1,100,2); [Hb2,wh]=freqz(bb2,1,100,2); subplot(332);

plot(wh,abs(Hb2),'-r',wh,abs(H2),'-b');grid; title('29点汉明窗')

subplot(335);

plot(wh,unwrap(angle(Hb2)),'-r',wh,unwrap(angle(H2)),'-b');grid; subplot(338);

plot(wh,20*log10(abs(Hb2)/sum(Hb2)),'-r',wh,20*log10(abs(H2)/sum(H2)),'-b');grid;

title('对数幅频特性')

bb=br.*(0.42-0.5*cos(2*pi*(0:1:28)/29)+0.08*cos(4*pi*(0:1:28)/29)); [H3,wb]=freqz(bb,1,100,2); subplot(333);

plot(wb,abs(Hb2),'-r',wb,abs(H3),'-b');grid; title('29点布莱克曼窗') subplot(336);

plot(wb,unwrap(angle(Hb2)),'-r',wb,unwrap(angle(H3)),'-b');grid; subplot(339);

plot(wh,20*log10(abs(Hb2)/sum(Hb2)),'-r',wb,20*log10(abs(H3)/sum(H3)),'-b');grid;

title('对数幅频特性')

29点矩形窗1.511.5129点汉明窗1.5129点布莱克曼窗0.50.50.5000.20.40.60.81000.20.40.60.81000.20.40.60.810-50-5-100-10-20-3000.20.40.60.81-4000.20.40.60.81-10-1500.20.40.60.81-20-15对数幅频特性0-20-40-60-8000.20.40.60.81-1000-50对数幅频特性0-50对数幅频特性-100-15000.20.40.60.81-15000.20.40.60.81

—Matlab自动设计 —窗函数设计

图1 长度为29点时滤波器的频率特性

分析：

由图可知，矩形窗设计的滤波器，幅频接近Matlab自动设计的滤波器的幅频特性，相频特性曲线和对数幅频特性曲线几乎与之重合。汉明窗设计的滤波器，幅频特性曲线与Matlab自动设计的滤

波器的幅频特性重合，相频特性曲线和对数幅频特性曲线接近Matlab自动设计的滤波器的特性曲线。对于Hamming窗和Blackman窗来说，当长度为29点时，滤波器的相频特性不理想。

（5） 将滤波器长度增大为41点，重新用这三种窗函数设计滤波器：

1hdn2ejejn1ejn2jn441jnded2 sinn4,N120n21,0n40矩形窗：wrn

0,其它sinn204,0n40 h1nhdnwrnn200,其它n2n0.0.46cos,n0,1,,40

汉明（Hamming）窗：whnN0,其它nsinn2040.0.46cos2n,0n40 h2nhdnwhnn20N0,其它n2n4n0.420.5cos0.08cos,n0,1,,40

布莱克曼（Blackman）窗：wbnNN0,其它nh3nhdnwbnsinn2040.420.5cos2n0.08cos4n,0n40 n20NN0,其它nMatlab代码及运算结果： clc; clear; N=40;

br=0.25*sinc((-20:1:20)/4);%计算滤波器系数

bb1=fir1(N,0.25,boxcar(N+1));%用fir1计算滤波器系数

[H1,wr]=freqz(br,1,100,2);%频域变换 [Hb1,wr]=freqz(bb1,1,100,2); subplot(331);

plot(wr,abs(Hb1),'-r',wr,abs(H1),'-b');grid;%画出幅频特性 title('41点矩形窗') subplot(334);

plot(wr,unwrap(angle(Hb1)),'-r',wr,unwrap(angle(H1)),'-b');grid;%画出相频特性 subplot(337);

plot(wr,20*log10(abs(Hb1)/sum(Hb1)),'-r',wr,20*log10(abs(H1)/sum(H1)),'-b');grid; title('对数幅频特性')

bh=br.*(0.-0.46*cos(2*pi*(0:1:40)/41)); bb2=fir1(N,0.25,hamming(N+1)); [H2,wh]=freqz(bh,1,100,2); [Hb2,wh]=freqz(bb2,1,100,2); subplot(332);

plot(wh,abs(Hb2),'-r',wh,abs(H2),'-b');grid; title('41点汉明窗') subplot(335);

plot(wh,unwrap(angle(Hb2)),'-r',wh,unwrap(angle(H2)),'-b');grid; subplot(338);

plot(wh,20*log10(abs(Hb2)/sum(Hb2)),'-r',wh,20*log10(abs(H2)/sum(H2)),'-b');grid; title('对数幅频特性')

bb=br.*(0.42-0.5*cos(2*pi*(0:1:40)/41)+0.08*cos(4*pi*(0:1:40)/41)); [H3,wb]=freqz(bb,1,100,2); subplot(333);

plot(wb,abs(Hb2),'-r',wb,abs(H3),'-b');grid; title('41点布莱克曼窗') subplot(336);

plot(wb,unwrap(angle(Hb2)),'-r',wb,unwrap(angle(H3)),'-b');grid; subplot(339);

plot(wh,20*log10(abs(Hb2)/sum(Hb2)),'-r',wb,20*log10(abs(H3)/sum(H3)),'-b');grid; title('对数幅频特性')

41点矩形窗1.511.5141点汉明窗141点布莱克曼窗0.50.50.5000.20.40.60.81000.20.40.60.81000.20.40.60.810-5-10-15-2000.20.40.60.810-100-10-20-20-3000.20.40.60.81-3000.20.40.60.81对数幅频特性0-20-40-60-8000.20.40.60.81-10000.2-500对数幅频特性0-50对数幅频特性-1000.40.60.81-15000.20.40.60.81

—Matlab自动设计 —窗函数设计

图表 1长度为41点时滤波器的频率特性

与图1对比，当N=41时，滤波器的幅频特性的过度带变得更窄，相频特性有所改善。还可看出，对于矩形窗来说，N的大小对于其幅频谱的过冲大小没有影响。

7.2 给定一理想带阻滤波器的频率特性

1,6Hdej0,

361,3重复7.1题各项要求。

解：

为了使设计出的滤波器hd(n)的对称中心在M/2，并保证滤波器具有线性相位，我们按照下

面的方法规定滤波器的相位。

Mexp（-j）,，Hdej263

0,其他对Hdej做IDTFT有

1hd(n)2Hd(ej)ejnd

MM12nn1/6j1j22ededed/6/322MMM sinnsinnsinn62232Mn2/3Mjn2对该序列加一个长度M=29的矩形窗，就得到了按矩形窗设计的FIR滤波器

sinhRnMMMnsinnsinn62232Mn2,n0,1,,28

clear all; close all; n=0:1:28

br=sin(pi/6*(n-14))./(pi*(n-14))+sin(pi*(n-14))./(pi*(n-14))-sin(pi/3*(n-14))./(pi*(n-14)) br(15)=0.833 stem(n,br)

求得滤波器系数为：h(0)——h(28)=

-0.0000 -0.0090 -0.0000 0.0106 -0.0000 -0.03 -0.06 -0.0621 -0.0000 0.0870 0.1378 0.1061 -0.0000 -0.1165 0.8330 -0.1165 -0.0000 0.1061 0.1378 0.0870 -0.0000 -0.0621 -0.06 -0.03 -0.0000 0.0106 -0.0000 -0.0090 -0.0000

若施加一个长度M=29的Hamming窗，就得到了按Hamming窗设计的FIR滤波器

sinhHnMn62MMsinnsinn232Mn20.0.46cos2n,n0,1,,28Mclear all; close all; n=0:1:28

br=sin(pi/6*(n-14))./(pi*(n-14))+sin(pi*(n-14))./(pi*(n-14))-sin(pi/3*(n-14))./(pi*(n-14)) br(15)=0.833

bh=br.*(0.-0.46*cos(2*pi*(0:1:28)/29)); stem(n,bh)

求得滤波器系数为：h(0)——h(28)=

-0.0000 -0.0090 -0.0000 0.0106 -0.0000 -0.03 -0.06 -0.0621 -0.0000 0.0870 0.1378 0.1061 -0.0000 -0.1165 0.8330 -0.1165 -0.0000 0.1061 0.1378 0.0870 -0.0000 -0.0621 -0.06 -0.03 -0.0000 0.0106 -0.0000 -0.0090 -0.0000

clear all; close all; n=0:1:28;

br=sin(pi/6*(n-14))./(pi*(n-14))+sin(pi*(n-14))./(pi*(n-14))-sin(pi/3*(n-14))./(pi*(n-14)) br(15)=0.833;

[H1,wr]=freqz(br,1,100,2);

f=[0 0.16 0.17 0.33 0.34 1]; % 给定频率轴分点； m=[1 1 0 0 1 1];% 给定在这些频率分点上理想的幅频响应 N1=28;

b1=fir2(N1,f,m);

[H2,wr]=freqz(b1,1,100,2)

subplot(321); stem(n,br,'-b') title('29点矩形窗系数')

subplot(323);

plot(wr,abs(H2),'-r',wr,abs(H1),'-b');grid; title('幅频特性') subplot(325);

plot(wr,unwrap(angle(H2)),'-r',wr,unwrap(angle(H1)),'-b');grid; title('相频特性')

bh=br.*(0.-0.46*cos(2*pi*(0:1:28)/29)); subplot(322); stem(n,bh)

title('29点汉明窗系数')

[H3,wr]=freqz(bh,1,100,2); subplot(324);

plot(wr,abs(H2),'-r',wr,abs(H3),'-b');grid; title('幅频特性') subplot(326);

plot(wr,unwrap(angle(H2)),'-r',wr,unwrap(angle(H3)),'-b');grid; title('相频特性')

29点矩形窗系数10.510.529点汉明窗系数00-0.5051015幅频特性202530-0.5051015幅频特性2025301.511.510.50.5000.20.40.60.81000.20.40.60.81相频特性0-200-20相频特性-40-40-6000.20.40.60.81-6000.20.40.60.81

—Matlab自动设计 —窗函数设计

图表 2 长度为29点的滤波器的频率特性

7.3 试用频域抽样法完成7.1及7.2题两个滤波器的设计。

解：

j利用频域抽样法设计数字滤波器时，需要将理想的连续频域相应Hde转换为理想的离

散频域响应Hdk，例如，对低通滤波器，

当N为偶数时，有

ejN1k/N,k0,1,,N/21 Hdk0,kN/2jN1k/N,kN/21,,N1ejN1k/N或令Hdke,k0,1,,N/21 jN1k/N则HdNke,k0,1,,N/21

Hdk0,kN/2

jN1k/N当N为奇数时，有Hdke,k0,1,,N1 jN1k/N,k0,1,,N1/2Hdke或令

HdNkHdk,k1,2,,N1/22jnk1N1根据hnnHdkeN,n0,1,,N1

Nk0即可求出所要设计的滤波器。上述赋值方法保证了所设计的滤波器具有线性相位。

j根据线性相位的条件，本题中N=29，为奇数，并且在02区间内，He关于偶

对称，选择第1种线性相位进行设计。

根据指标要求，在02内有29个采样点，所以第K点对应的频率为几个特殊的频率：位于

2k。 2942273和4之间，同理，位于第25个和第26个采样点之间，29294求得Hk：

ad =

1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1

kN122kN28k,0k28 292jnk1N1hnHkeN,n0,1,,N1Nk01HkeNk0

N1jN12ej2nkN 将H(k)代入，得到h(n)：

,n0,1,,N1clc;

clear;

k1=[0:3];k2=[4:25];k3=[26:28]; a1=0*k1+1; a2=0*k2; a3=0*k3+1;

ad(1:29)=[a1,a2,a3] N=29; k=0:(N-1); wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd)) stem(k,h);

title('29点抽样h(n)图像');

h =

-0.0321 -0.0147 0.0114 0.0328 0.0376 0.0215 -0.0097 -0.0415 -0.0562 -0.0405 0.00 0.0823

0.1602 0.2193 0.2414 0.2193 0.1602 0.0823 0.00 -0.0405 -0.0562 -0.0415 -0.0097 0.0215 0.0376 0.0328 0.0114 -0.0147 -0.0321

clc; clear;

k1=[0:3];k2=[4:25];k3=[26:28]; a1=0*k1+1; a2=0*k2; a3=0*k3+1;

ad(1:29)=[a1,a2,a3] N=29; k=0:(N-1);

wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd)) subplot(221); stem(k,h);

title('29点抽样h(n)的图像'); [Hd,p]=freqz(h,1,100,2); figure(1); subplot(222); plot(p,abs(Hd));

title('29点频率抽样幅频特性'); subplot(223); plot(p,angle(Hd));

title('29点频率抽样相频特性'); subplot(224);

plot(p,20*log10(abs(Hd)/sum(Hd)));grid on; title('29点频率抽样对数幅频特性')

29点矩形窗2100-10-200-50-1002100-10-200-50-10029点汉明窗10.500-20-400-100-20029点布莱克曼窗00.5100.5100.5100.5对数幅频特性100.5对数幅频特性100.5对数幅频特性100.5100.5100.51

比较两图可以看出，用频域抽样法设计的滤波器的幅频响应没有用窗函数法设计的幅频响

应的曲线平滑，因此，对频域抽样法需要改进，改进的方法是在Hdk由1变为0的频率出加过度点。

插一个点：

插两个点：

29点抽样h(n)的图像0.30.20.10-0.101020300.51.5129点频率抽样幅频特性000.5129点频率抽样相频特性420-2-400.51-1000-5029点频率抽样对数幅频特性-15000.51

结论：差值要注意，不一定越多越好，规律有待总结。 clc; clear;

k1=[0:2];k2=[4:24];k3=[26:28]; a1=0*k1+1; a2=0.4 a3=0*k2; a4=0.4

a5=0*k3+1;

ad(1:29)=[a1,a2,a3,a4,a5] N=29; k=0:(N-1); wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd)) subplot(221); stem(k,h,'g-');hold

title('29点抽样h(n)的图像'); [Hd,p]=freqz(h,1,100,2); figure(1); subplot(222);

plot(p,abs(Hd),'g-');hold

title('29点频率抽样幅频特性'); subplot(223);

plot(p,angle(Hd),'g-');hold

title('29点频率抽样相频特性'); subplot(224);

plot(p,20*log10(abs(Hd)/sum(Hd)),'g-');grid on;hold title('29点频率抽样对数幅频特性') clc; clear;

k1=[0:2];k2=[5:24];k3=[27:28]; a1=0*k1+1; a2=0.4 a3=0.3 a4=0*k2; a5=0.3 a6=0.4

a7=0*k3+1;

ad(1:29)=[a1,a2,a3,a4,a5,a6,a7] N=29; k=0:(N-1);

wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd)) subplot(221);

stem(k,h,'r-');grid on;hold title('29点抽样h(n)的图像'); [Hd,p]=freqz(h,1,100,2); figure(1); subplot(222);

plot(p,abs(Hd),'r-');grid on;hold title('29点频率抽样幅频特性'); subplot(223);

plot(p,angle(Hd),'r-');grid on;hold title('29点频率抽样相频特性'); subplot(224);

plot(p,20*log10(abs(Hd)/sum(Hd)),'r-');grid on;hold title('29点频率抽样对数幅频特性')

clc; clear;

k1=[0:2];k2=[3:4];k3=[5:9];k4=[10:25];k5=[26:28]; a1=0*k1+1; a2=0*k2; a3=0*k3+1; a4=0*k4; a5=0*k5+1;

ad(1:29)=[a1,a2,a3,a4,a5]; N=29; k=0:(N-1); wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd)) stem(k,h);

title('29点抽样h(n)图像');

clc; clear;

k1=[0:2];k2=[3:4];k3=[5:9];k4=[10:25];k5=[26:28]; a1=0*k1+1; a2=0*k2; a3=0*k3+1; a4=0*k4; a5=0*k5+1;

ad(1:29)=[a1,a2,a3,a4,a5]; N=29; k=0:(N-1); wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd)) subplot(221); stem(k,h);

title('29点抽样h(n)的图像'); [Hd,p]=freqz(h,1,100,2); figure(1); subplot(222); plot(p,abs(Hd));

title('29点频率抽样幅频特性'); subplot(223); plot(p,angle(Hd));

title('29点频率抽样相频特性'); subplot(224);

plot(p,20*log10(abs(Hd)/sum(Hd)));grid on; title('29点频率抽样对数幅频特性')

比较两图可以看出，用频域抽样法设计的滤波器的幅频响应没有用窗函数法设计的幅频响应的曲线平滑，因此，对频域抽样法需要改进，改进的方法是在Hdk由1变为0的频率出加过度点。

7.4一滤波器的理想频率响应如图5所示：

图5

（1） 试用窗函数法设计该滤波器，要求具有线性相位，滤波器长度为33，用汉宁窗。 （2） 用频率抽样法设计，仍要求具有线性相位，滤波器长度为33，过渡点自行设置。 要求先用手算得出hn，然后上机求Hej。

解：线性相位理想滤波器的频率特性：

2j1e,22j1,0e2j Hde

21ej,0221ej,21hdn212212jjn1eed22jjn1ede202012jjn1eed22jjn1ede2

1144cosn8cosn,n222n21,n2N1其中,,当N33时,162（1） 窗函数法

(a) 采用33点汉宁窗设计滤波器：

2n0.50.5cos,n0,1,,32

窗函数:wh1nN0,其它nh1nhdnwh1n144cosn168cosn16211cos2n,n0,1,,15,17,18,,32222332n162 12n0.50.5cos,n162330,其它n

clc clear n=0:1:32

x=(4+4*cos((n-16)*pi)-8*cos(0.5*pi*(n-16)))./(2*pi*pi*(n-16).^2) b=x.*(0.5-0.5*cos(2*pi*n/33))

b(17)=0.5 stem(n,b);grid;

滤波器系数h1(0)——h1(32)：

0 0 0.0003 0 0 0 0.0029 0 0 0 0.0163 0 0 0 0.1958 0 0.5000 0 0.1958 0 0 0 0.0163 0 0 0 0.0029 0 0 0 0.0003 0 0

滤波器的频率特性：

clc clear n=0:1:32

x=(4+4*cos((n-16)*pi)-8*cos(0.5*pi*(n-16)))./(2*pi*pi*(n-16).^2) b=x.*(0.5-0.5*cos(2*pi*n/33)) b(17)=0.5 stem(n,b);grid;

[H1,wr]=freqz(b,1,100,2); subplot(121);

plot(wr,abs(H1));grid; title('幅频特性') subplot(122);

plot(wr,unwrap(angle(H1)));grid; title('相频特性')

图6 33点汉宁窗设计的滤波器频率特性

(b) 采用33点矩形窗设计滤波器：

1,0n32 窗函数:wrn0,其它h2nhdnwrn144cosn168cosn162,n0,1,,15,17,18,,3222n162 12,n160,其它n

clc clear n=0:1:32

x=(2+2*cos((n-16)*pi)-4*cos(0.5*pi*(n-16)))./(pi*pi*(n-16).^2)

滤波器系数：h2(0)——h2(32) =

0 0 0.0041 0 0 0 0.0081 0 0 0 0.0225 0 0 0 0.2026 0 0.5000 0 0.2026 0 0 0 0.0225 0 0 0 0.0081 0 0 0 0.0041 0 0

（c）采用33点三角窗

2n,n0,1,,16 窗函数:wtn32w(32n),n17,18,,32h3nhdnwtn144cosn168cosn1622n,n0,1,,152322n1621 ,n1621n168cosn1644cos2232n,n17,18,,322322n1620,其它n滤波器系数：h3(0)——h3(32)=

0 0 0.0007 0 0 0 0.0033 0 0 0 0.0146 0 0 0 0.1788 0 0.5000 0 0.1788 0 0 0 0.0146 0 0 0 0.0033 0 0 0 0.0007 0 0

33点矩形窗设计滤波器10.80.60.90.80.70.60.50.40.200.40.30.200.20.40.60.810.100.20.40.60.8133点三角窗设计滤波器相频特性(wraped)43210-1-2-3-400.20.40.60.8143210-1-2-3-400.2相频特性(wraped)0.40.60.81

图7 33点矩形窗和三角窗设计的滤波器的频率特性

（d）采用65点矩形窗和三角窗设计滤波器 方法与33点相同，得到滤波器的频率特性：

65点矩形窗设计滤波器10.80.610.80.665点三角窗设计滤波器0.40.200.40.2000.20.40.60.8100.20.40.60.81相频特性(wraped)43210-1-2-3-400.20.40.60.8143210-1-2-3-400.2相频特性(wraped)0.40.60.81

图8 65点矩形窗和三角窗设计的滤波器的频率特性

结论：对于本题来说，采用多点矩形窗设计的FIR数字滤波器，其频率特性接近于理想滤波器。 （2） 用频率抽样法设计滤波器 （a）33点抽样

j根据线性相位的条件，本题中N=33，为奇数，并且在02区间内，He关于偶

对称，选择第1种线性相位进行设计。

根据指标要求，在02内有33个采样点，所以第K点对应的频率为几个特殊的频率：

2k。 33228和9之间，同理，位于第16个和第17个采样点之间，位于333323位于第24个和第26个采样点之间。根据对称性，求得Hk： 2 H0= 1.0000 H1=H32=0.8788 H2=H31=0.7576 H3=H30=0.63

H4=H29=0.5152 H5=H28=0.3939 H6=H27=0.2727 H7=H26=0.1515

H8=H25=0.0303 H9=H24=0.0909 H10=H23=0.2121 H11=H22=0.3333

H12=H21=0.45 H13=H20=0.5758 H14=H19=0.6970 H15=H18=0.8182 H16=H17=0.9394

kN122k32N33k,0k32 hn1N1j2nkNHkeN,n0,1,,N1k0N1 将H(k)代入，得到h(n)：

1jN12j2NHkeeNnk,n0,1,,N1k0

k1=[0:8];k2=[9:16];k3=[17:24];k4=[25:32]; a1=(-4/33)*k1+1; a2=(4/33)*k2-1; a3=(-4/33)*k3+3; a4=(4/33)*k4-3;

ad(1:33)=[a1,a2,a3,a4]; N=33; k=0:(N-1); wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd)) stem(k,h);

title('33点抽样h(n)图像');

h(n)=

0.0005 0.0027 0.0043 0.0006 0.0005 0.0018 0.0083 0.0005 0.0014 0.0227 0.0007 0.0005 0.0011 0.2028 0.5005 0.0009 0.2028 0.0011 0.0005 0.0007 0.0227 0.0005 0.0006 0.0083 0.0018 0.0005 0.0006 0.0043 0.0005 h(n)的图像：

0.0006 0.0009 0.0014 0.0027

图9 33点频率抽样时h(n)的图像

MATLAB代码：

k1=[0:8];k2=[9:16];k3=[17:24];k4=[25:32]; a1=(-4/33)*k1+1; a2=(4/33)*k2-1; a3=(-4/33)*k3+3; a4=(4/33)*k4-3;

ad(1:33)=[a1,a2,a3,a4]; N=33; k=0:(N-1); wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd)) subplot(221); stem(k,h);

title('33点抽样h(n)的图像 '); [Hd,p]=freqz(h,1,100,2); figure(1); subplot(222); plot(p,abs(Hd));

title('33点频率抽样幅频特性'); subplot(223); plot(p,angle(Hd));

title('33点频率抽样相频特性'); subplot(224);

plot(p,20*log10(abs(Hd)/sum(Hd)));grid on; title('33点频率抽样对数幅频特性')

滤波器频率特性：

图10 33点频率抽样时滤波器的频率特性

（b）65点抽样

j 根据线性相位的条件， N=65，为奇数，并且在02区间内，He关于偶对

称，选择第1种线性相位进行设计。

根据指标要求，在02内有65个采样点，所以第K点对应的频率为几个特殊的频率：间，

2k。 652216和17之间，同理，位于第32个和第33个采样点之位于653323位于第48个和第49个采样点之间。根据对称性，求得Hk： 2H0——H=

1.0000 0.9385 0.8769 0.81 0.7538 0.6923 0.6308 0.5692 0.5077 0.4462 0.3846 0.3231 0.2615 0.2000 0.1385 0.0769 0.01 0.0462 0.1077 0.1692 0.2308 0.2923 0.3538 0.41 0.4769 0.5385 0.6000 0.6615 0.7231 0.7846 0.8462 0.9077 0.9692 0.9692 0.9077 0.8462 0.7846 0.7231 0.6615 0.6000 0.5385 0.4769 0.41 0.3538 0.2923 0.2308 0.1692 0.1077 0.0462 0.01 0.0769 0.1385 0.2000 0.2615 0.3231 0.3846 0.4462 0.5077 0.5692 0.6308 0.6923 0.7538 0.81 0.8769 0.9385

kN122kNk,0k 652jnk1N1hnHkeN,n0,1,,N1Nk01HkeNk0N1N1j2e2jnkN 将H(k)代入，得到h(n)：

,n0,1,,N1k1=[0:16];k2=[17:32];k3=[33:48];k4=[49:]; a1=(-4/65)*k1+1; a2=(4/65)*k2-1; a3=(-4/65)*k3+3; a4=(4/65)*k4-3;

ad(1:65)=[a1,a2,a3,a4]; N=65; k=0:(N-1); wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd)) stem(k,h);

title('65点抽样h(n)图像');

h(0)—— h()=

0.0001 0.0007 0.0009 0.0001 0.0005 0.0017 0.0001 0.0004 0.0042 0.0001 0.0003 0.0226 0.5001 0.0002 0.2027 0.0001 0.0002 0.0081 0.0001 0.0002 0.0025 0.0001 0.0002 0.0012 0.0001

0.0001 0.0001 0.0002 0.0001 0.0002 0.0001 0.0002 0.0001 0.0003 0.0001 0.0003 0.0001 0.0004 0.0001 0.0006 0.0001 0.0006 0.0012 0.0004 0.0025 0.0003 0.0081 0.0003 0.2027 0.0002 0.0226 0.0002 0.0042 0.0002 0.0017 0.0001 0.0009 0.0002 0.0002 0.0002 0.0002 0.0003 0.0004 0.0005 0.0007

k1=[0:16];k2=[17:32];k3=[33:48];k4=[49:]; a1=(-4/65)*k1+1; a2=(4/65)*k2-1; a3=(-4/65)*k3+3; a4=(4/65)*k4-3;

ad(1:65)=[a1,a2,a3,a4]; N=65; k=0:(N-1); wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd)) subplot(221); stem(k,h);

title('65点抽样h(n)的图像 '); [Hd,p]=freqz(h,1,100,2); figure(1); subplot(222); plot(p,abs(Hd));

title('65点频率抽样幅频特性'); subplot(223); plot(p,angle(Hd));

title('65点频率抽样相频特性'); subplot(224);

plot(p,20*log10(abs(Hd)/sum(Hd)));grid on; title('65点频率抽样对数幅频特性')

滤波器的频率特性：

对比

结论：用频率抽样法设计数字滤波器时，增加采样点的数量，可以减小逼近误差。

clc; clear;

ka1=[0:8];ka2=[9:16];ka3=[17:24];ka4=[25:32];

aa1=(-4/33)*ka1+1; aa2=(4/33)*ka2-1; aa3=(-4/33)*ka3+3; aa4=(4/33)*ka4-3;

ada(1:33)=[aa1,aa2,aa3,aa4]; Na=33; ka=0:(Na-1); wma=2*pi*ka/Na;

hda=ada.*exp(-j*(Na-1)/2*wma); ha=real(ifft(hda))

k1=[0:16];k2=[17:32];k3=[33:48];k4=[49:]; a1=(-4/65)*k1+1; a2=(4/65)*k2-1; a3=(-4/65)*k3+3; a4=(4/65)*k4-3;

ad(1:65)=[a1,a2,a3,a4]; N=65; k=0:(N-1); wm=2*pi*k/N;

hd=ad.*exp(-j*(N-1)/2*wm); h=real(ifft(hd))

subplot(221);

stem(ka,ha,'b-');hold; stem(k,h,'g-');

title('抽样h(n)的图像')

[Hda,pa]=freqz(ha,1,100,2); [Hd,p]=freqz(h,1,100,2); subplot(222);

plot(pa,abs(Hda),'b-',p,abs(Hd),'g-'); title('频率抽样幅频特性'); subplot(223);

plot(pa,angle(Hda),'b-',p,angle(Hd),'g-'); title('频率抽样相频特性'); subplot(224);

plot(pa,20*log10(abs(Hda)/sum(Hda)),'b-',p,20*log10(abs(Hd)/sum(Hd)),'g-');grid on;

title('频率抽样对数幅频特性')