数字滤波器设计
题:
对模拟信号进行低通滤波处理,要求通带0≤f≤5kHz,通带衰减小于0.5dB,阻带5.5kHz≤f<∞,阻带衰减大于50dB,设采样频率Fs=20kHz。
(1)设计巴特沃斯模拟低通滤波器,求出Ha(s)的分子、分母多项式系数B和A,并画出幅频响应损耗函数曲线。
分别用脉冲响应不变法和双线性变换法设计IIR低通数字滤波器,求出Ha(z) 的分子、分母多项式系数Bz和Az,并画出幅频响应损耗函数曲线
采用窗函数法(分别用汉宁窗、哈明窗、布莱克曼窗函数)设计满足要求的FIR低通滤波器,求出h(n),并画出幅频响应损耗函数曲线.
用频率采样法设计满足要求的FIR低通滤波器,求出h(n),并画出幅频响应损耗函数曲线。
具体内容如下:
(1)设计巴特沃斯模拟低通滤波器,求出Ha(s)的分子、分母多项式系数B和A,并画出幅频响应损耗函数曲线。
程序:
wp=2*pi*5000;
ws=2*pi*5800;
Rp=0.5;
As=50;
[N,wc]=buttord(wp,ws,Rp,As,'s');
[B,A]=butter(N,wc,'s');
k=0:511;
fk=0:20000/512:20000;
wk=2*pi*fk;
Hk=freqs(B,A,wk);
plot(fk/1000,20*log10(abs(Hk)));
grid on
xlabel('频率/kHz');
ylabel('幅度/dB');
axis([0,6,-65,5]);
波形图:
A = 1.0e+207 *
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0020 2.1576
B = 1.0e+207 *
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.1576
N = 46
(2)分别用脉冲响应不变法和双线性变换法设计IIR低通数字滤波器,求出Ha(z) 的分子、分母多项式系数Bz和Az,并画出幅频响应损耗函数曲线
脉冲响应不变法
程序:
Fs=20000;
wp=10000*pi;
ws=11600*pi;
Rp=0.5;
As=50;
[N,wc]=buttord(wp,ws,Rp,As,'s');
[B,A]=butter(N,wc,'s');
[Bz,Az]=impinvar(B,A);
k=0:511;
fk=0:20000/512:20000;
wk=2*pi*fk;
Hk=freqs(B,A,wk);
plot(fk/1000,20*log10(abs(Hk)));
grid on;
xlabel('频率/kHz');
ylabel('幅值/dB');
axis([0,6,-65,5]);
波形图:
Bz = 1.0e-007 *
0 -0.0000 0.0000 -0.0000 0.0001 -0.0006 0.0035 -0.0116 0.0279 -0.0745 0.1490 -0.2235 0.3353 -0.3725 0.4470 -0.4098 0.3353 -0.2235 0.1304 -0.0698 0.0291 -0.0093 0.0026 -0.0006 0.0001 -0.0000 0.0000 -0.0000
Az = 1.0e+007 *
0.0000 -0.0000 0.0000 -0.0003 0.0018 -0.0081 0.0296 -0.0888 0.2219 -0.4684 0.8431 -1.3028 1.7370 -2.0041 2.0040 -1.7367 1.3024 -0.8427 0.4681 -0.2217 0.0887 -0.0296 0.0081 -0.0018 0.0003 -0.0000 0.0000 -0.0000
N = 46
双线性变换法
程序:
Fs=20000;wpz=10000 /Fs;
wsz=11600/Fs;
Rp=0.5;
As=50;
wp=2*tan(wpz*pi /2);
ws=2*tan(wsz*pi /2);[N,wc]=buttord(wp,ws,Rp,As,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs);
[Nd,wdc]=buttord(wpz,wsz,Rp,As);
[Bdz,Adz]=butter(Nd,wdc);
k=0:511;
fk=0:20000/512:20000;
wk=2*pi*fk;Hk=freqs(B,A,wk);plot(fk/1000,20*log10(abs(Hk)));
grid on;xlabel('频率/kHz');
ylabel('幅值/dB');axis([0,16,-2800,5])
波形图:
Bz = 1.0e-007 *
0 -0.0000 0.0000 -0.0000 0.0001 -0.0006 0.0035 -0.0116 0.0279 -0.0745 0.1490 -0.2235 0.3353 -0.3725 0.4470 -0.4098 0.3353 -0.2235 0.1304 -0.0698 0.0291 -0.0093 0.0026 -0.0006 0.0001 -0.0000 0.0000 -0.0000
Az = 1.0e+007 *
0.0000 -0.0000 0.0000 -0.0003 0.0018 -0.0081 0.0296 -0.0888 0.2219 -0.4684 0.8431 -1.3028 1.7370 -2.0041 2.0040 -1.7367 1.3024 -0.8427 0.4681 -0.2217 0.0887 -0.0296 0.0081 -0.0018 0.0003 -0.0000 0.0000 -0.0000
N = 27
(3)采用窗函数法(分别用汉宁窗、哈明窗、布莱克曼窗函数)设计满足要求的FIR低通滤波器,求出h(n),并画出幅频响应损耗函数曲线.
a )汉宁窗
程序:
Fs=20000;
fp=5000;
fs=5800;
m=[1 1 0 0];
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Rp=0.5;
As=50;
Bt=ws-wp;
N0 = ceil (6.6* pi /Bt);
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;
hn=fir1(N-1,wc,hanning(N));
freqz(hn,1,512);
plot (w,20*log(abs(hn)));
grid;
axis ([0 ,1, -1000 , 100]);
xlabel ('频率/kHz');
ylabel ('幅值/dB' );
波形图:
hn =
0.0000 -0.0000 -0.0000 0.0002 -0.0000 -0.0004 0.0002 0.0007 -0.0006 -0.0010 0.0012 0.0012 -0.0021 -0.0010 0.0032 0.0005 -0.0045 0.0006 0.0057 -0.0025 -0.0066 0.0050 0.0070 -0.0083 -0.0065 0.0123 0.0047 -0.0168 -0.0012 0.0215 -0.0046 -0.0262 0.0134 0.0307 -0.0270 -0.0345 0.0497 0.0375 -0.0974 -0.0394 0.3154 0.5400 0.3154 -0.0394 -0.0974 0.0375 0.0497 -0.0345 -0.0270 0.0307 0.0134 -0.0262 -0.0046 0.0215 -0.0012 -0.0168 0.0047 0.0123 -0.0065 -0.0083 0.0070 0.0050 -0.0066 -0.0025 0.0057 0.0006 -0.0045 0.0005 0.0032 -0.0010 -0.0021 0.0012 0.0012 -0.0010 -0.0006 0.0007 0.0002 -0.0004 -0.0000 0.0002 -0.0000 -0.0000 0.0000
b )哈明窗
程序:
Fs=20000;
fp=5000;
fs=5800;
m=[1 1 0 0];
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Rp=0.5;
Rs=50;
Bt=ws-wp;
N0 = ceil (6.6* pi /Bt);
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;
hn=fir1(N-1,wc,hamming(N));
freqz(hn,1,512);
plot (w,20*log(abs(hn)));
grid on;
axis ([0 ,1, -1000 , 100]);
xlabel ('频率/kHz');
ylabel ('幅值/dB' );
波形图:
hn =
0.0003 -0.0006 -0.0001 0.0008 -0.0001 -0.0010 0.0004 0.0012 -0.0008 -0.0014 0.0016 0.0015 -0.0025 -0.0012 0.0037 0.0005 -0.0049 0.0007 0.0061 -0.0026 -0.0069 0.0052 0.0072 -0.0086 -0.0066 0.0125 0.0048 -0.0170 -0.0012 0.0217 -0.0046 -0.0264 0.0135 0.0308 -0.0271 -0.0346 0.0498 0.0375 -0.0975 -0.0394 0.3154 0.5401 0.3154 -0.0394 -0.0975 0.0375 0.0498 -0.0346 -0.0271 0.0308 0.0135 -0.0264 -0.0046 0.0217 -0.0012 -0.0170 0.0048 0.0125 -0.0066 -0.0086 0.0072 0.0052 -0.0069 -0.0026 0.0061 0.0007 -0.0049 0.0005 0.0037 -0.0012 -0.0025 0.0015 0.0016 -0.0014 -0.0008 0.0012 0.0004 -0.0010 -0.0001 0.0008 -0.0001 -0.0006 0.0003
c )布莱克曼窗
程序:
Fs=20000;
fp=5000;
fs=5800;
m=[1 1 0 0];
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Rp=0.5;
Rs=50;
Bt=ws-wp;
N0 = ceil (6.6* pi /Bt);
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;
hn=fir1(N-1,wc,blackman(N));
freqz(hn,1,512);
plot (w,20*log(abs(hn)));
grid;
axis ([0 ,1, -1000 , 100]);
xlabel ('频率/kHz');
ylabel ('幅值/dB' );
波形图:
hn =
-0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0001 0.0001 0.0002 -0.0002 -0.0004 0.0005 0.0005 -0.0009 -0.0005 0.0016 0.0002 -0.0024 0.0004 0.0034 -0.0015 -0.0043 0.0034 0.0049 -0.0060 -0.0049 0.0095 0.0037 -0.0137 -0.0010 0.0186 -0.0040 -0.0237 0.0124 0.0288 -0.0257 -0.0333 0.0485 0.0369 -0.0965 -0.0392 0.3150 0.5400 0.3150 -0.0392 -0.0965 0.0369 0.0485 -0.0333 -0.0257 0.0288 0.0124 -0.0237 -0.0040 0.0186 -0.0010 -0.0137 0.0037 0.0095 -0.0049 -0.0060 0.0049 0.0034 -0.0043 -0.0015 0.0034 0.0004 -0.0024 0.0002 0.0016 -0.0005 -0.0009 0.0005 0.0005 -0.0004 -0.0002 0.0002 0.0001 -0.0001 -0.0000 0.0000 -0.0000 -0.0000 -0.0000
(4)用频率采样法设计满足要求的FIR低通滤波器,求出h(n),并画出幅频响应损耗函数曲线。
程序:
T=input('T=')
fp=5000;
fs=5800;
Fs=20000;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Rp=0.5;
As=50;
Bt=ws-wp;
m=1;
N=ceil((m+1)*2*pi/Bt);
N=N+mod(N+1,2);
Np=fix(wp/(2*pi/N));
Ns=N-2*Np-1;
Hk=[ones(1,Np+1),zeros(1,Ns),ones(1,Np)];
Hk(Np+2)=T;
Ak(N-Np)=T;
thetak=-pi*(N-1)*(0:N-1)/N;
hdk=Hk.*exp(j*thetak);
hn=real(ifft(hdk));
hw=fft(hn,1024);
wk=2*pi*[0:1023]/1024;
hgw=hw.*exp(j*wk*(N-1)/2);
Rp=max(20*log10(abs(hgw)));
hgmin=min(real(hgw));
As=20*log10(abs(hgmin));
[N,wc]=buttord(wp,ws,Rp,As,'s');
[B,A]=butter(N,wc,'s');
k=0:511;fk=0:20000/512:20000;wk=2*pi*fk;
Hk=freqs(B,A,wk);
plot(fk/1000,20*log10(abs(Hk)));
grid on;
xlabel('频率/kHz');
ylabel('幅度/dB');
T=1
T =
1
Rp =
0.2801
As =
-29.6745
波形图:
hn =
0.0004 0.0012 -0.0023 -0.0027 0.0045 0.0040 -0.0069 -0.0051 0.0098 0.0061 -0.0133 -0.0070 0.0175 0.0078 -0.0230 -0.0084 0.0305 0.0089 -0.0416 -0.0093 0.0609 0.0096 -0.1044 -0.0097 0.3178 0.5098 0.3178 -0.0097 -0.1044 0.0096 0.0609 -0.0093 -0.0416 0.0089 0.0305 -0.0084 -0.0230 0.0078 0.0175 -0.0070 -0.0133 0.0061 0.0098 -0.0051 -0.0069 0.0040 0.0045 -0.0027 -0.0023 0.0012 0.0004