数字信号处理课内实验记录
实验一 利用DFT分析信号频谱
一、实验目的
- 加深对DFT 原理的理解。
- 应用DFT 分析信号的频谱。
- 深刻理解利用DFT 分析信号频谱的原理，分析实现过程中出现的现象及解决方法。
二、实验原理
1.DFT与DTFT的关系
有限长序列
x
(
n
)
(
0
⩽
n
⩽
N - 1
)
x(n)(0 \leqslant {\text{n}} \leqslant {\text{N - 1}})
x(n)(0⩽n⩽N - 1)的离散时间傅里叶变换
X
(
e
j
ω
)
X(e^{j\omega})
X(ejω)在频率区间
(
0
⩽
ω
⩽
2
π
)
(0 \leqslant \omega \leqslant {\text{2}}\pi )
(0⩽ω⩽2π)的N个等间隔分布的点
k
ω
=
2
π
k
/
N
(
0
⩽
k
⩽
N - 1
)
k\omega=2\pi k/N(0 \leqslant {\text{k}} \leqslant {\text{N - 1}})
kω=2πk/N(0⩽k⩽N - 1)上的N个取样值可以由下式表示：
X
(
e
j
ω
)
∣
ω
=
2
π
k
/
N
=
∑
k
=
0
N
−
1
x
(
n
)
e
−
j
2
π
N
k
n
=
X
(
k
)
,
0
⩽
k
⩽
N
−
1
X({e^{j\omega }}){|_{\omega = 2\pi k/N}} = \sum\limits_{k = 0}^{N - 1} {x(n){e^{ - j\frac{{2\pi }}{N}kn}} = X(k),0 \leqslant k \leqslant N - 1}
X(ejω)∣ω=2πk/N=k=0∑N−1x(n)e−jN2πkn=X(k),0⩽k⩽N−1
由上市可知，序列
x
(
n
)
x(n)
x(n)的N点DFT
X
(
k
)
X(k)
X(k)，实际上就是
x
(
n
)
x(n)
x(n)序列的DTFT在N个等间隔频率点
k
ω
=
2
π
k
/
N
(
0
⩽
k
⩽
N - 1
)
k\omega=2\pi k/N(0 \leqslant {\text{k}} \leqslant {\text{N - 1}})
kω=2πk/N(0⩽k⩽N - 1)上样本
X
(
k
)
X(k)
X(k)
2.利用DFT求DTFT
方法一：
方法二：
然而在实际MATLAB计算中，上述插值运算不见得是最好的办法。由于DFT是DTFT的取样值，其相邻两个频率样本点的间距为 2 π / N 2\pi/N 2π/N ，所以如果我们增加数据的长度N，使得到的DFT谱线就更加精细，其包络就越接近DTFT的结果，这样就可以利用DFT来近似计算DTFT。如果没有更多的数据，可以通过补零来增加数据长度。
利用DFT分析连续时间信号的频谱
采用计算机分析连续时间信号的频谱，第一步就是把连续时间信号离散化，这里需要进行两个操作：一是采样，二是截断。因此，可以将利用DFT 分析连续非周期信号频谱的步骤归纳如下：
- 确定时域采样间隔T，得到离散序列x(n)；
- 确定截取长度M，得到M点离散序列 x M ( n ) = x ( n ) w ( n ) x_M(n)=x(n)w(n) xM(n)=x(n)w(n)，这里 w ( n ) w(n) w(n)为窗函数。
- 确定频域采样点数N，要求 N ⩾ M N \geqslant M N⩾M。
- 利用FFT计算离散序列的D点DFT，得到 X M ( k ) X_M(k) XM(k)。
- 由 X M ( k ) X_M(k) XM(k)计算 X a ( j Ω ) X_a(j\Omega) Xa(jΩ)采样点的近似值。
采用上述方法计算 x a ( t ) x_a(t) xa(t)的频谱，需要注意如下三个问题：
- 频谱混叠。如果不满足采样定理的条件，频谱会出现混叠误差。对于频谱无限宽的信号，应考虑覆盖大部分主要频率分量的范围。
- 栅栏效应和频谱分辨率。使用DFT 计算频谱，得到的结果只是N 个频谱样本值，样本值之间的频谱是未知的，像通过一个栅栏观察频谱，称为“栅栏效应”。频谱分辨率与记录长度成反比，要提高频谱分表率，就要增加记录时间。
- 频谱泄露。对信号截断会把窗函数的频谱引入信号频谱，照成频谱泄露。解决这个问题的主要办法是采用旁瓣小的窗函数，频谱泄露和窗函数均会引起误差。因此，要合理选取采样间隔和截取长度，必要时还需考虑加适当的窗。
对于连续时间周期信号，我们在采用计算机进行计算时，也总是要进行截断，序列总是有限长的，仍然可以采用上述方法近似计算。
三、实验内容
1.已知 x ( n ) = { 2 → , − 1 , 1 , 1 } x(n) = \{\underrightarrow 2, - 1,1,1\} x(n)={2,−1,1,1}，完成如下要求：
（1）计算其DFT，并画出 [ − π , π ] [-\pi,\pi] [−π,π]区间的波形
clc;
clear;
close all
%%DTFT波形
n = 0:3;
x=[2,-1,1,1];
w = -pi:0.01*pi:pi;
X = x*exp(-j*n'*w)%计算频谱
subplot(211);
plot(w,abs(X));%绘制DTFT波形
xlabel('\Omega/\pi');
title('幅频');
axis tight
subplot(212);
plot(w,angle(X)/pi);
xlabel('\Omega/\pi');
title('相频');
axis tight
（2）计算4点DFT，并把结果显示在（1）所画的图中
clc;
clear;
close all
%%4点DFT及波形
n = 0:3;
x=[2,-1,1,1];
w = -pi:0.01*pi:pi;
X = x*exp(-j*n'*w)%计算频谱
y = fft(x,4);
subplot(211);
hold on
plot(w,abs(X));%绘制DTFT波形
xlabel('\Omega/\pi');
title('幅频');
axis tight
stem(0:3,abs(y),'fill')
subplot(212);
hold on
plot(w,angle(X)/pi);
xlabel('\Omega/\pi');
title('相频');
axis tight
stem(0:3,angle(y)/pi,'fill')
（3）对x(n)补零，计算64点DFT，并显示结果
clc;
clear;
close all
n = 0:3;
x=[2,-1,1,1];
x= [x,zeros(1,60)];
y = fft(x,64);
w = -pi:0.01*pi:pi;
subplot(211);
stem(0:63,abs(y),'fill')
xlabel('\Omega/\pi');
title('幅频');
axis tight
subplot(212);
stem(0:63,angle(y)/pi,'fill')
xlabel('\Omega/\pi');
title('相频');
axis tight
（4）根据实验结果，分析是否可以由DFT计算DTFT，如果可以，如何实现。
可以由DFT计算DTFT。通过补零加长序列，提高采样密度，可以由DFT近似计算DTFT。
2.考察序列
x
(
n
)
=
c
o
s
(
0.48
π
n
)
+
c
o
s
(
0.52
π
n
)
x(n)=cos(0.48\pi n)+cos(0.52\pi n)
x(n)=cos(0.48πn)+cos(0.52πn)
（1）
0
⩽
n
⩽
10
0 \leqslant n \leqslant 10
0⩽n⩽10时，用DFT估计x(n)的频谱；将x(n)补零加长到长度为100点序列用DFT估计x(n)的频谱，要求画出相应波形。
clc;
clear;
close all
%%0<=n<=10
n = 0:10;
x = cos(0.48*pi*n)+cos(0.52*pi*n);
y = fft(x);
subplot(221);
stem(0:10,abs(y),'filled');
title('幅频');
subplot(222);
stem(0:10,angle(y)/pi,'filled');
title('相频')
%%补零
x = [x,zeros(1,89)];
y = fft(x);
subplot(223);
stem(0:99,abs(y),'filled');
title('幅频（补零）');
subplot(224);
stem(0:99,angle(y)/pi,'filled');
title('相频（补零）');
（2） 0 ⩽ n ⩽ 100 0 \leqslant n \leqslant 100 0⩽n⩽100时，用DFT估计x(n)的频谱，并画出波形
clc;
clear;
close all
%%0<=n<=10
n = 0:100;
x = cos(0.48*pi*n)+cos(0.52*pi*n);
y = fft(x);
subplot(211);
stem(0:100,abs(y),'filled');
title('幅频');
subplot(212);
stem(0:100,angle(y)/pi,'filled');
title('相频')
（3）根据实验结果，分析怎样提高频谱分辨率
频谱分辨率与信号采集时间间隔成反比，因此可以通过增加时域内信号采样点数和增大截取长度来增加分辨率。
3.已知信号 x ( t ) = 0.15 s i n ( 2 π f 1 t ) + s i n ( 2 π f 2 t ) − 0.1 s i n ( 2 π f 3 t ) x(t)=0.15sin(2\pi f_1 t)+sin(2\pi f_2 t)-0.1sin(2\pi f_3 t) x(t)=0.15sin(2πf1t)+sin(2πf2t)−0.1sin(2πf3t)，其中， f 1 = 1 H z f_1 = 1Hz f1=1Hz， f 2 = 2 H z f_2 = 2Hz f2=2Hz， f 3 = 3 H z f_3 = 3Hz f3=3Hz。从x(t)的表达式可以看出，它是包含三个频率的正弦波，但是，从其时域比兴来看，似乎是一个正弦信号，利用DFT做频谱分析，确定适合的参数，使得得到的频谱的频率分辨率符合需要。
clc;
clear;
close all
f1 = 1;
f2 = 2;
f3 = 3;
fs = 100;%采样频率
Tp = 1;%截取长度
t = 0:1/fs:Tp;
x = 0.15*sin(2*pi*f1*t)+sin(2*pi*f2*t)-0.1*sin(2*pi*f3*t);
N = Tp*fs+1;%采样点数
X = fft(x,N);
w = 0:fs/(N-1):fs;
subplot(211);
stem(w,abs(X),"filled");
xlabel('f/Hz');
title('幅频');
subplot(212);
stem(w,angle(X),"filled");
xlabel('\Omega/\pi');
title('相频');
由奈奎斯特采样定律，采样频率必须大于等于信号最高频率的两倍，因此采样频率必须大于等于6Hz，得到的频谱的频率分辨率才符合需要。
通过选取合适的采样周期，可以完整地恢复出原信号的频谱波形。
4.利用DFT近似分析连续事件信号 x ( t ) = e − 0.1 t u ( t ) x(t)=e^{-0.1t}u(t) x(t)=e−0.1tu(t)的频谱（幅度谱）。分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定适合的参数。
clc;
clear;
close all
T = 1;%采样间隔
Tp = 100;%截取长度
n = 0:T:Tp;
x = exp(-0.1*n);
X = fft(x);
subplot(211);
stem(n,abs(X),"filled");
xlabel('f/Hz');
title('幅频');
subplot(212);
stem(n,angle(X),"filled");
xlabel('\Omega/\pi');
title('相频');
采样间隔为1s，截取长度为100s时：
采样间隔为5s，截取长度为100s时：
采样间隔为25s，截取长度为100s时：
采样间隔为1s，截取长度为50s时：
采样间隔为1s，截取长度为25s时：
对比上图，最终确定合适的参数为：采样间隔5s，截取长度100s。
四、实验心得体会
通过本次实验，我们掌握并加深对DFT原理的理解并且学会应用DFT分析信号频谱在此基础上利用DFT分析信号频谱的原理,掌握了利用matlab分析现实问题的步骤及办法。并且，通过这次的实验，对信号序列有了更加深刻的认识，单纯的一个信号序列是没有意义的，只有配合他本身的时间序列才是一个完整的信号序列，才可以对其进行分析。
实验二 利用FFT计算线性卷积
一、实验目的
- 掌握利用FFT计算线性卷积的原理及具体实现方法。
- 加深理解重叠相加法和重叠保留法。
- 考察利用FFT计算线性卷积各种方法的适用范围。
二、实验原理
1. 线性卷积与圆周卷积
设
x
(
n
)
x(n)
x(n)为
L
L
L点序列，
h
(
n
)
h(n)
h(n)为
M
M
M点序列，
x
(
n
)
x(n)
x(n)和
h
(
n
)
h(n)
h(n)的线性卷积为：
y
l
(
n
)
=
x
(
n
)
∗
h
(
n
)
=
∑
m
=
−
∞
∞
x
(
m
)
h
(
n
−
m
)
{y_l}(n) = x(n)*h(n) = \sum\limits_{m = - \infty }^\infty {x(m)h(n - m)}
yl(n)=x(n)∗h(n)=m=−∞∑∞x(m)h(n−m)
y
l
(
n
)
y_l(n)
yl(n)的长度为
L
+
M
−
1
L+M-1
L+M−1。
x ( n ) x(n) x(n)和 h ( n ) h(n) h(n)的 N N N点圆周卷积为：
圆周卷积与线性菊娜姬相等而不产生交叠的必要条件为：
N
⩾
L
+
M
−
1
N \geqslant L + M - 1
N⩾L+M−1
圆周卷积定理：根据DFT的性质，
x
(
n
)
x(n)
x(n)和
h
(
n
)
h(n)
h(n)的N点圆周卷积的DFT等于它们DFT的乘积：
2. 快速卷积
快速卷积算法用圆周卷积实现线性卷积，根据圆周卷积定理利用FFT算法实现圆周卷积。可以将快速卷积的步骤归纳如下：
-
为了使线性卷积可以用圆周卷积来计算，必须选择 N ⩾ L + M − 1 N \geqslant L + M - 1 N⩾L+M−1；同时为了能使用基-2FFT完成卷积运算，要求$N=2^\nu $。采用补零的办法使x(n)和h(n)的长度均为N。
-
计算x(n)和h(n)的N点FFT
-
组成乘积
-
利用IFFT计算Y(k)的IDFT，得到线性卷积y(n)
3.分段卷积
考察单位取样响应为h(n)的线性系统，输入为x(n)，输出为y(n)，则
y
l
(
n
)
=
x
(
n
)
∗
h
(
n
)
{y_l}(n) = x(n)*h(n)
yl(n)=x(n)∗h(n)
当输入序列x(n)极长时，如果要等x(n)全部集齐时再开始进行卷积，会使输出相对输入有较大的延时，再者如果序列太长，需要大量的存储单元。为此，我们把x(n)分段，分别求出每段的卷积，合在一起得到最后总的输出。这种方法称为分段卷积。分段卷积可细分为重叠相加法和重叠保留法。
重叠保留法：设x(n)的长度为 N x N_x Nx，h(n)的长度为M。我们把序列x(n)分成多段N点序列 x i ( n ) x_i(n) xi(n)，每段与前一段重叠 M − 1 M-1 M−1个样本。由于第一段没有前一段保留信号,为了修正，我们在第一个输入段前面填充 M − 1 M-1 M−1个零。计算每一段与h(n)的圆周卷积，则其每段卷积结果的前 M − 1 M-1 M−1个样本不等于线性卷积值，不是正确的样本值。所以我们将每段卷积结果的前 M − 1 M-1 M−1个样本舍去，只保留后面的 N − M + 1 N-M+1 N−M+1个正确输出样本，把这些输出样本合起来，得到总得输出。
重叠相加法：设h(n)长度为M ，将信号x(n)分解成长为L的子段，建议L选择与的M数量级相同。以
x
i
(
n
)
x_i(n)
xi(n)表示每段信号，则:
每一段卷积 y i ( n ) y_i(n) yi(n)的长度为 L + M − 1 L+M-1 L+M−1，所以在做求和时，相邻两段序列有 M − 1 M-1 M−1个样本重叠，即前一段的最后 M − 1 M-1 M−1个样本和下一段的前 M − 1 M-1 M−1个序列重叠，这个重叠部分相加，再与不重叠部分共同组成输出y(n)。
三、实验内容
假设要计算序列 x ( n ) = u ( n ) − u ( n − L ) , 0 ⩽ n ⩽ L x(n)=u(n)-u(n-L),0 \leqslant n \leqslant L x(n)=u(n)−u(n−L),0⩽n⩽L和 h ( n ) = c o s ( 0.2 π n ) , 0 ⩽ n ⩽ M h(n)=cos(0.2\pi n),0 \leqslant n \leqslant M h(n)=cos(0.2πn),0⩽n⩽M的线性卷积，完成以下实验内容
1.设 L = M L=M L=M。根据线性卷积的表达式和快速卷积的原理，分别编成实现计算两个序列线性卷积的方法，比较当序列长度分别为8，16，32，64，256，1024时，两种方法计算线性卷积所需时间。
clc;
clear;
close all
for M = [8 16 32 64 256 512 1024]
L = M;
n = 0:1:L-1;
x = ones(1,L);
h = cos(0.2.*pi.*n);
disp('L=M=');disp(M)
disp('线性卷积用时：')
tic
y1 = conv(x,h);
toc
disp('快速卷积用时：')
tic
X = fft(x);
H = fft(h);
Y2 = X.*H;
y2 = ifft(Y2);
toc
end
实验结果如下：
L=M=
8
线性卷积用时：
历时 0.000011 秒。
快速卷积用时：
历时 0.000041 秒。
L=M=
16
线性卷积用时：
历时 0.000004 秒。
快速卷积用时：
历时 0.000024 秒。
L=M=
32
线性卷积用时：
历时 0.000029 秒。
快速卷积用时：
历时 0.000016 秒。
L=M=
64
线性卷积用时：
历时 0.000028 秒。
快速卷积用时：
历时 0.000025 秒。
L=M=
256
线性卷积用时：
历时 0.000047 秒。
快速卷积用时：
历时 0.000093 秒。
L=M=
512
线性卷积用时：
历时 0.000069 秒。
快速卷积用时：
历时 0.000033 秒。
L=M=
1024
线性卷积用时：
历时 0.000218 秒。
快速卷积用时：
历时 0.000038 秒。
>>
对比可知，长度相同的情况下，快速卷积算法所需时间更短。
在本实验中，因为程序启动时设备的运算速率与稳定运行时运算速率有一定差距，因此起初的几组数据运算时间会更长。
2.当L=2048且M=256时，比较直接计算线性卷积和快速卷积所需的时间，进一步考察当L=4096且M=256时两种算法所需的时间。
clc;
clear;
close all
L = 2048;
M = 256;
n1 = 0:1:L;
n2 = 0:1:M;
x = ones(1,L);
h = cos(0.2.*pi.*n2);
disp('L=2048,M=256')
disp('线性卷积用时：')
tic
y1 = conv(x,h);
toc
disp('快速卷积用时：')
N = M+L-1;
tic
X = fft(x,N);
H = fft(h,N);
Y2 = X.*H;
y2 = ifft(Y2,N);
toc
实验结果如下：
当L=2048,M=256
L=2048,M=256
线性卷积用时：
历时 0.000131 秒。
快速卷积用时：
历时 0.000203 秒。
>>
当L=4096,M=256
L=4096,M=256
线性卷积用时：
历时 0.000360 秒。
快速卷积用时：
历时 0.003145 秒。
>>
3.编程实现利用重叠相加法计算两个序列的线性卷积，考察L=2048且M=256时计算线性卷积的时间，与第2题的结果进行比较。
clc;
clear;
close all
%重叠相加法
L = 2048;
M = 256;
n = 0:1:M-1;
x = ones(1,L);
h = cos(0.2.*pi.*n);
disp('L=');
disp(L);
disp('M=');
disp(M);
Li = M;%分段的段长
Num = L/Li;%段数
Ni = Li+M-1;%每段卷积长度
N = L+M-1;%卷积长度
y = zeros(1,N);
disp('重叠相加法卷积用时：');
tic
H = fft(h,Ni);
for i=1:Num
i_st = (i-1)*M+1;
xi = x(i_st:i_st+M-1);
Xi = fft(xi,Ni);
Yi = Xi.*H;
yi = ifft(Yi,Ni);
y(i_st:i_st+Ni-1) = y(i_st:i_st+Ni-1)+yi;%重叠相加
end
toc
结果如下：
L=
2048
M=
256
重叠相加法卷积用时：
历时 0.000255 秒。
>>
与第2题中相比，这种方法更快。
4.编程实现利用重叠保留法计算两个序列的线性卷积，考察L=2048且M=256时计算线性卷积的时间，与第2 题的结果进行比较。
clc;
clear;
close all
%重叠保留法
L = 2048;
M = 256;
n = 0:1:M-1;
x = ones(1,L);
x_add = [zeros(1,M-1),x];
h = cos(0.2.*pi.*n);
disp('L=');
disp(L);
disp('M=');
disp(M);
Li = M;%分段的段长
Num = L/Li;%段数
Ni = Li+M-1;%每段卷积长度
N = L+M-1;%卷积长度
y = zeros(1,N);
disp('重叠保留法卷积用时：');
tic
H = fft(h,Ni);
for i=1:Num
i_st = (i-1)*M+1;
xi = x_add(i_st:i_st+Li-1);
Xi = fft(xi,Ni);
Yi = Xi.*H;
yi = ifft(Yi,Ni);
y(i_st:i_st+Li-1) = yi(M:Ni);
end
toc
结果如下：
L=
2048
M=
256
重叠保留法卷积用时：
历时 0.000271 秒。
>>
与第2题中相比，这种方法更快。
四、实验心得体会
本次实验要求我们掌握利用FFT计算线性卷积的原理及具体实现方法，以及通过对快速卷积和线性卷积运算速度的比较，更加直观的看到利用FFT快速卷积的优点。本次实验让我切实看到了FFT算法的高效性及对于庞大的数据的处理实时性，同时对重叠保留法、重叠相加法也有了更深的认识，对以后的实验有了很好的理论基础，受益颇多。
实验三 IIR数字滤波器设计
一、实验目的
- 掌握利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方法。
- 加深理解数字滤波器和模拟滤波器之间的技术指标转化。
- 掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及适用范围。
二、实验原理
IIR滤波器具有无限长持续时间脉冲响应，而模拟滤波器一般都具有无限长的脉冲响应，因此它与模拟滤波器相匹配。IIR滤波器设计的基本方法就是先设计一个合适的模拟滤波器，然后利用复值映射把模拟滤波器变换成数字滤波器。
1. 数字滤波器和模拟滤波器的一些指标
下图中画出了数字滤波器的幅频特性和指标。其中 [ 0 , ω p ] [0,\omega_p] [0,ωp]称为通带， [ ω s t , π ] [\omega_{st},\pi] [ωst,π]称为阻带， [ ω p , ω s t ] [\omega_p,\omega_{st}] [ωp,ωst]为过滤带， δ 1 {\delta _1} δ1为通带相应中的容限， δ 2 \delta_2 δ2为阻带容限， R p R_p Rp为通带波动， A s A_s As为阻带衰减，且：
下图中画出了模拟滤波器的技术指标。其中 ε \varepsilon ε为通带波动系数， Ω p \Omega_p Ωp是单位弧度/秒的通带截止频率，A为以dB为单位的阻带衰减参数， Ω s t \Omega_{st} Ωst是单位为弧度/秒的阻带截止频率。参数 ε \varepsilon ε和A与滤波器的通带波动 R p R_p Rp及阻带衰减 A s A_s As之间有如下关系：
2. 模拟原型滤波器
IIR滤波器设计方法由已有的模拟滤波器得到数字滤波器，我们将这些模拟滤波器称为原型滤波器。常用的模拟原型滤波器有巴特沃斯（Butterworth）滤波器、切比雪夫（Chebyshev，Ⅰ型和Ⅱ型）滤波器和椭圆（Ellipse）滤波器等。
2.1巴特沃斯滤波器
N阶巴特沃斯低通滤波器的幅度平方响应为：
其中，N为滤波器阶数， Ω c \Omega_c Ωc是截止频率（单位：rad/s）。
系统函数
H
a
(
s
)
H_a(s)
Ha(s)为：
分母多项式的极点为：
2.2切比雪夫低通滤波器
切比雪夫滤波器有两种，切比雪夫Ⅰ型滤波器在通带中具有等波动响应，而切比雪夫Ⅱ型滤波器在阻带中具有等波动响应。相比于单调特性的滤波器，选择等波动特性的滤波器，可以得到较低阶数的滤波器，例如巴特沃斯滤波器在通带和阻带上的响应均是单调的，因此对于相同的指标，切比雪夫滤波器的阶数比巴特沃斯滤波器的阶数低。
切比雪夫Ⅰ型滤波器
切比雪夫Ⅰ型滤波器的幅度平方响应为：
其中
N
N
N为滤波器阶数，
Ω
c
\Omega_c
Ωc为截止频率，$\varepsilon
为通带波动系数，它与
为通带波动系数，它与
为通带波动系数，它与R_p$的关系如下：
c
N
(
x
)
c_N(x)
cN(x)是N阶切比雪夫多项式。
确定系统函数 H a ( s ) H_a(s) Ha(s)的过程与前面巴特沃斯滤波器里介绍的方法类似，先求出 H a ( s ) H_a(s) Ha(s) H a ( − s ) H_a(-s) Ha(−s)的极点，然后把左半平面的极点分配给 H a ( s ) H_a(s) Ha(s)。通过求下式的根得到 H a ( s ) H_a(s) Ha(s) H a ( − s ) H_a(-s) Ha(−s)的极点
如果 p k = σ k + j Ω k , k = 0 , 1 , . . . , N − 1 p_k=\sigma_k+j\Omega_k,k=0,1,...,N-1 pk=σk+jΩk,k=0,1,...,N−1是上式的根（左半平面），则
其中
这些根落在长轴为 b Ω c b\Omega_c bΩc短轴为 a Ω c a\Omega_c aΩc的椭圆上，最后可以得到系统函数为
其中K是归一化因子，可以由下式决定
在MATLAB中，设计阶数为N，截止频率为 Ω c \Omega_c Ωc，通带波动为 R p R_p Rp的切比雪夫Ⅰ型莫您滤波器可以通过两种方法实现：
**方法1：**MATLAB 提供了函数[z,p,k] = cheb1ap(N,Rp)来设计一个阶数为N，通带波动为 R p R_p Rp的归一化（截止频率为1）切比雪夫Ⅰ型模拟原型滤波器，它返回数阻z和p（零点和积点）以及增益k。切比雪夫Ⅰ型滤波器系统函数不存在零点，将归一化切比雪夫Ⅰ型滤波器系统函数的极点倍乘以 Ω c \Omega_c Ωc增益倍乘以s=0时非归一化滤波器系统函数分母多项式和归一化滤波器系统函数多项式的比值，即可实现由归一化切比雪夫Ⅰ型滤波器到非归一化切比雪夫Ⅰ型滤波器的变换。
**方法2：**在MATLAB 中还可以直接利用cheby1函数设计切比雪夫Ⅰ型模拟原型滤波器。cheby1函数既可用于设计模拟切比雪夫Ⅰ型滤波器也可用于设计数字切比雪夫Ⅰ型滤波器，用于设计模拟切比雪夫Ⅰ型波器时，其调用格式为:
[b,a] = cheby1(N,Rp,Wn,‘s’) N表示滤波器的阶数，Rp表示通带波动，Wn表示截止频率，'s’表示设计模拟滤波器。
切比雪夫Ⅱ型滤波器
切比雪夫Ⅱ型滤波器有一个单调的通带和等波动的阻带，这就意味着这种滤波器在s平面上既有零点又有极点，因此通带中的群延迟特性比切比雪夫Ⅰ型滤波器要好。切比雪夫Ⅱ型滤波器的幅度平方响应为
在MATLAB中，设计阶数为N，截止频率为 Ω c \Omega_c Ωc，通带波动为 R p R_p Rp的切比雪夫Ⅱ型莫您滤波器可以通过两种方法实现：
**方法1：**MATLAB 提供了函数[z,p,k] = cheb2ap(N,As)来设计一个阶数为N ，阻带衰减为 A s A_s As的归一化（截止频率为1）切比雪夫Ⅱ型模拟原型滤波器，它返回数阻z 和p（零点和积点）以及增益k。切比雪夫Ⅱ型滤波器系统函数存在零点，可以将归一化切比雪夫Ⅱ型滤波器系统函数的极点和零点倍乘以 Ω c \Omega_c Ωc，增益倍乘以s=0时非归一化的有理函数和归一化的有理函数的比值，即可实现由归一化切比雪夫Ⅱ型滤波器到非归一化切比雪夫Ⅱ型滤波器的变换。
**方法2：**在MATLAB 中还可以直接利用cheby2函数设计切比雪夫Ⅱ型模拟原型滤波器。cheby2函数既可用于设计模拟切比雪夫Ⅰ型滤波器也可用于设计数字切比雪夫Ⅱ型滤波器，用于设计模拟切比雪夫Ⅱ型波器时，其调用格式为:
[b,a] = cheby2(N,Rp,Wn,‘s’) N表示滤波器的阶数，Rp表示通带波动，Wn表示截止频率，'s’表示设计模拟滤波器。
2.3椭圆滤波器
椭圆滤波器的通带和阻带均具有等波动响应，因此相对于巴特沃斯滤波器和切比雪夫滤波器，对于给定相同的指标，它可使阶数最小，换言之对于给定的阶数N，它能使过渡带最陡。
椭圆滤波器的幅度平方响应为
其中N为滤波器阶数， Ω c \Omega_c Ωc为截止频率， ε \varepsilon ε为通带波动系数， U N ( x ) U_N(x) UN(x)是N阶雅可比椭圆函数。
在MATLAB中，设计阶数为N，截止频率为 Ω c \Omega_c Ωc，通带波动为 R p R_p Rp，阻带衰减为 A s A_s As的椭圆模拟原型滤波器可以通过两种方法实现：
**方法1：**MATLAB提供了函数[z,p,k] = ellipap(N,Rp,As)来设计一个阶数为N ，通带波动为 R p R_p Rp ，阻带衰减为 A s A_s As的归一化（截止频率为1）椭圆模拟原型滤波器，它返回数阻z和p（零点和积点）以及增益k。将归一化椭圆滤波器系统函数的极点和零点倍乘以 Ω c \Omega_c Ωc，增益倍乘以s=0时非归一化的有理函数和归一化的有理函数的比值，即可实现由归一化椭圆滤波器到非归一化椭圆滤波器的变换。
**方法2：**在MATLAB中还可以直接利用ellip函数设计椭圆模拟原型滤波器。ellip函数既可用于设计模拟椭圆滤波器也可用于设计数字椭圆滤波器，用于设计模拟椭圆波器时，其调用格式为:
[b,a] = ellip(N,Rp,As,Wn,‘s’) N表示滤波器的阶数，Rp表示通带波动，As表示阻带衰减，Wn表示截止频率，'s’表示设计模拟滤波器。
3.模拟滤波器到数字滤波器的变换
从模拟滤波器到数字滤波器的变换就是要由 H a ( s ) H_a(s) Ha(s)进一步求得H(z)，也就是由s平面到z平面的变换，这种变换要求数字滤波器能模仿模拟滤波器的特性。我们最常用的从模拟滤波器到数字滤波器的变换方法有两种：脉冲响应不变法和双线性变换法。
3.1脉冲响应不变法
从时域响应出发，使数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的单位脉冲响应 h a ( t ) h_a(t) ha(t)， h ( n ) h(n) h(n)等于 h a ( t ) h_a(t) ha(t)的取样值。
变换方法：
-
将 H a ( s ) H_a(s) Ha(s)进行部分展开：
-
对 H a ( s ) H_a(s) Ha(s)进行拉氏变换：
-
对 h a ( t ) h_a(t) ha(t)时域采样得到 h ( n ) h(n) h(n)
-
对 h ( n ) h(n) h(n)进行z变换
设计步骤：
-
确定数字滤波器性能指标 ω p , ω s t , R p , A s \omega_p,\omega_{st},R_p,A_s ωp,ωst,Rp,As。
-
将数字滤波器频率指标转换成相应的模拟滤波器频率指标
Ω p = ω p / T , Ω s t = ω s t / T \Omega_p=\omega_p/T,\Omega_{st}=\omega_{st}/T Ωp=ωp/T,Ωst=ωst/T -
根据指标设计模拟滤波器
-
将 H a ( s ) H_a(s) Ha(s)展成部分分式形式
-
把模拟极点 p k p_k pk转换成数字极点 e p k T e^{p_kT} epkT，得到数字滤波器
3.2双线性变换法
从频率响应出发，直接使数字滤波器的频率响应逼近模拟滤波器的频率响应，进而求得数字滤波器的系统函数H(z)。
-
通过正切变换 Ω = t a n ( Ω 1 T / 2 ) \Omega=tan(\Omega_1T/2) Ω=tan(Ω1T/2)将s平面的整个虚轴 j Ω j\Omega jΩ压缩到s1平面 j Ω 1 j\Omega_1 jΩ1轴上的 − j π / T -j\pi/T −jπ/T到 j π / T j\pi/T jπ/T段上。
-
通过 z = e s 1 T z=e^{s_1T} z=es1T变换将 s 1 s_1 s1平面映射到z平面上，得到s平面和z平面的映射关系为：
双线性变换的映射关系如图：
一般来说为了调节模拟频带与数字频带之间的关系，可引入常数c，使映射关系为：
设计时可采用使模拟频率特性与数字频率特性在低频处有确切对应关系的常数 c = 2 / T c=2/T c=2/T。
设计步骤：
-
确定数字滤波器性能指标 ω p , ω s t , R p , A s \omega_p,\omega_{st},R_p,A_s ωp,ωst,Rp,As。
-
根据频率的非线性关系，确定预期的模拟滤波器频率指标
Ω p = 2 / T t a n ( ω p / 2 ) , Ω s t = 2 / T t a n ( ω s t / 2 ) \Omega_p=2/Ttan(\omega_p/2),\Omega_{st}=2/Ttan(\omega_{st}/2) Ωp=2/Ttan(ωp/2),Ωst=2/Ttan(ωst/2) -
根据指标设计模拟滤波器
-
将双线性关系带入，得到
三、实验内容
设采样频率为 f s = 10 k H z f_s=10kHz fs=10kHz，设计数字低通滤波器，满足如下指标
- 通带截止频率： f p = 1 k H z f_p=1kHz fp=1kHz
- 阻带截止频率： f s t = 1.5 k H z f_{st}=1.5kHz fst=1.5kHz
- 通带波动： R p = 1 d B R_p=1dB Rp=1dB
- 阻带衰减： A s = 15 d B A_s=15dB As=15dB
要求分别设计巴特沃斯、切比雪夫Ⅰ型、切比雪夫Ⅱ型和椭圆模拟原型滤波器，并分别结合脉冲响应不变法和双线性变换法进行设计。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及适用范围。
代码中提供了脉冲响应不变法和双线性变换法两种方法，通过注释其中一个来选择另外一个。
巴特沃斯滤波器
clc;
clear;
close all
%巴特沃斯滤波器
fs = 10*1000;%采样率
fp = 1*1000;%通带截止频率
fst = 1.5*1000;%阻带截止频率
Rp = 1;%通带波动
As = 15;%阻带衰减
Wp = 2*pi*fp/fs;%归一化通带截止频率
Ws = 2*pi*fst/fs;%归一化阻带截止频率
N = ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws)));%阶数
Wc = Wp/((10^(Rp/10)-1)^(1/(2*N)));
[b,a] = butter(N,Wc,'s');%巴特沃斯滤波器系数
%[bz,az] = impinvar(b,a,1);%脉冲响应不变法转化为数字滤波器
[bz,az] = bilinear(b,a,1);%双线性变换法转化为数字滤波器
w=[0:500]*pi/500;%频率范围
[H,w]=freqz(bz,az);%频率响应
subplot(221);
plot(w/pi,abs(H));
grid on;
xlabel('\omega(\pi)');
ylabel('|H(e^j^\omega)|');
subplot(222);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
grid on;
xlabel('\omega(\pi)');
ylabel('|H(e^j^\omega)|(dB)');
subplot(223);
plot(w/pi,angle(H)/pi);
grid on;
xlabel('\omega(\pi)');
ylabel('Phase of H(e^j^\omega)(\pi)');
grd = grpdelay(bz,az,w);
subplot(224);
plot(w/pi,grd);
grid on;
xlabel('\omega(\pi)');
ylabel('Group delay');
disp('bz=');
disp(bz);
disp('az=')
disp(az)
脉冲响应不变法
bz=
-0.0000 0.0006 0.0101 0.0161 0.0041 0.0001 0
az=
1.0000 -3.3635 5.0684 -4.2759 2.1066 -0.5706 0.0661
>>
观察可知，满足题目要求的通带波动和阻带衰减。
双线性变换法
bz=
0.0005 0.0030 0.0074 0.0099 0.0074 0.0030 0.0005
az=
1.0000 -3.3960 5.1623 -4.3843 2.1721 -0.5911 0.0687
>>
观察可知，满足题目要求的通带波动和阻带衰减。
切比雪夫Ⅰ型滤波器
clc;
clear;
close all
%切比雪夫1型滤波器
fs = 10*1000;%采样率
fp = 1*1000;%通带截止频率
fst = 1.5*1000;%阻带截止频率
Rp = 1;%通带波动
As = 15;%阻带衰减
Wp = 2*pi*fp/fs;%归一化通带截止频率
Ws = 2*pi*fst/fs;%归一化阻带截止频率
e = sqrt(10^(Rp/10)-1);
A = 10^(As/20);
N = ceil(acosh(sqrt(A^2-1)/e)/(acosh(Ws/Wp)));%阶数
[b,a] = cheby1(N,Rp,Wp,'s');%切比雪夫1型滤波器系数
%[bz,az] = impinvar(b,a,1);%脉冲响应不变法转化为数字滤波器
[bz,az] = bilinear(b,a,1);%双线性变换法转化为数字滤波器
w=[0:500]*pi/500;%频率范围
[H,w]=freqz(bz,az);%频率响应
subplot(221);
plot(w/pi,abs(H));
grid on;
xlabel('\omega(\pi)');
ylabel('|H(e^j^\omega)|');
subplot(222);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
grid on;
xlabel('\omega(\pi)');
ylabel('|H(e^j^\omega)|(dB)');
subplot(223);
plot(w/pi,angle(H)/pi);
grid on;
xlabel('\omega(\pi)');
ylabel('Phase of H(e^j^\omega)(\pi)');
grd = grpdelay(bz,az,w);
subplot(224);
plot(w/pi,grd);
grid on;
xlabel('\omega(\pi)');
ylabel('Group delay');
disp('bz=');
disp(bz);
disp('az=')
disp(az)
脉冲响应不变法
bz=
-0.0000 0.0054 0.0181 0.0040 0
az=
1.0000 -3.0591 3.8323 -2.2919 0.5495
>>
观察可知，满足题目要求的通带波动和阻带衰减。
双线性变换法
bz=
0.0016 0.0065 0.0098 0.0065 0.0016
az=
1.0000 -3.0928 3.9012 -2.3402 0.5610
>>
切比雪夫Ⅱ型滤波器
clc;
clear;
close all
%切比雪夫2型滤波器
fs = 10*1000;%采样率
fp = 1*1000;%通带截止频率
fst = 1.5*1000;%阻带截止频率
Rp = 1;%通带波动
As = 15;%阻带衰减
Wp = 2*pi*fp/fs;%归一化通带截止频率
Ws = 2*pi*fst/fs;%归一化阻带截止频率
e = sqrt(10^(Rp/10)-1);
A = 10^(As/20);
N = ceil(acosh(sqrt(A^2-1)/e)/(acosh(Ws/Wp)));%阶数
[b,a] = cheby2(N,As,Ws,'s');%切比雪夫2型滤波器系数
[bz,az] = impinvar(b,a,1);%脉冲响应不变法转化为数字滤波器
%[bz,az] = bilinear(b,a,1);%双线性变换法转化为数字滤波器
w=[0:500]*pi/500;%频率范围
[H,w]=freqz(bz,az);%频率响应
subplot(221);
plot(w/pi,abs(H));
grid on;
xlabel('\omega(\pi)');
ylabel('|H(e^j^\omega)|');
subplot(222);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
grid on;
xlabel('\omega(\pi)');
ylabel('|H(e^j^\omega)|(dB)');
subplot(223);
plot(w/pi,angle(H)/pi);
grid on;
xlabel('\omega(\pi)');
ylabel('Phase of H(e^j^\omega)(\pi)');
grd = grpdelay(bz,az,w);
subplot(224);
plot(w/pi,grd);
grid on;
xlabel(' omega(\pi)');
ylabel('Group delay');
disp('bz=');
disp(bz);
disp('az=')
disp(az)
脉冲响应不变法
bz=
-0.2435 0.6762 -0.5578 0.3282 0.0166
az=
1.0000 -1.6658 1.4289 -0.5193 0.0935
>>
观察图像可知，通带波动不满足题目要求。
双线性变换法
bz=
0.1729 -0.1320 0.2624 -0.1320 0.1729
az=
1.0000 -1.7139 1.5036 -0.5665 0.1209
>>
观察频率响应可知，通带波动和阻带衰减均满足题目要求。
椭圆模拟原型滤波器
clc;
clear;
close all
%椭圆模拟原型滤波器
fs = 10*1000;%采样率
fp = 1*1000;%通带截止频率
fst = 1.5*1000;%阻带截止频率
Rp = 1;%通带波动
As = 15;%阻带衰减
Wp = 2*pi*fp/fs;%归一化通带截止频率
Ws = 2*pi*fst/fs;%归一化阻带截止频率
e = sqrt(10^(Rp/10)-1);
A = 10^(As/20);
N = ceil(acosh(sqrt(A^2-1)/e)/(acosh(Ws/Wp)));%阶数
[b,a] = ellip(N,Rp,As,Wp,'s');%椭圆模拟原型滤波器系数
[bz,az] = impinvar(b,a,1);%脉冲响应不变法转化为数字滤波器
%[bz,az] = bilinear(b,a,1);%双线性变换法转化为数字滤波器
w=[0:500]*pi/500;%频率范围
[H,w]=freqz(bz,az);%频率响应
subplot(221);
plot(w/pi,abs(H));
grid on;
xlabel('\omega(\pi)');
ylabel('|H(e^j^\omega)|');
subplot(222);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
grid on;
xlabel('\omega(\pi)');
ylabel('|H(e^j^\omega)|(dB)');
subplot(223);
plot(w/pi,angle(H)/pi);
grid on;
xlabel('\omega(\pi)');
ylabel('Phase of H(e^j^\omega)(\pi)');
grd = grpdelay(bz,az,w);
subplot(224);
plot(w/pi,grd);
grid on;
xlabel(' omega(\pi)');
ylabel('Group delay');
disp('bz=');
disp(bz);
disp('az=')
disp(az)
脉冲响应不变法
bz=
0.0809 -0.1552 0.2296 -0.1916 0.1031
az=
1.0000 -2.9706 3.7605 -2.2886 0.5799
>>
观察图像可知，通带波动不满足题目要求。
双线性变换法
bz=
0.1752 -0.4635 0.6438 -0.4635 0.1752
az=
1.0000 -3.0140 3.8430 -2.3499 0.5963
>>
观察图像可知，通带波动不满足题目要求。
脉冲响应不变法的优点是频率坐标的变换是线性的，因此如果模拟滤波器的频响是限带于折叠频率以内的话，通过变换后的数字滤波器的频响可以不失真地反映原响应与频率之间的关系。但是，其最大的缺点是有频谱的周期延拓效应，会产生混叠失真。
因此，脉冲响应不变法只能用于限带的频响特性，如衰减特性较好的低通或带通，而且高频衰减越大，频响的混叠效应就越小。
双线性变换法的优点是消除了脉冲响应不变法的固有的混叠失真，但缺点是频率变换之间的非常严重的非线性，会使变换后时域上的图像产生非常严重的失真。
因此，双线性变换法适用于低通数字滤波器。
四、实验心得体会
通过实验,学习了利用脉冲响应不变法设计IIR数字滤波器的基本原理，即从时域响应出发使数字滤波器的单位脉冲响应 h ( n ) h(n) h(n)模仿模拟滤波器的单位冲激响应 h a ( t ) h_a(t) ha(t)， h ( n ) h(n) h(n)等于 h a ( t ) h_a(t) ha(t)的取样值。所以设计中，应先根据数字指标转换来的模拟指标，设计模拟滤波器，再利用impinvar函数实现脉冲响应不变法模拟滤波器到数字滤波器的变换。脉冲响应不变法和双线性法是设计IIR数字滤波器设计的基础，掌握了这俩种方法能够使我们对数字滤波器设计有更加深刻的认识。
实验四 FIR数字滤波器设计
一、实验目的
掌握窗函数法和频率取样法设计FIR数字滤波器的原理及具体方法。
二、实验原理
1.线性相位FIR数字滤波器
线性相位条件：
线性相位的FIR数字滤波器的单位脉冲相应h(n)具有如下特性
h
(
n
)
=
h
(
N
−
1
−
n
)
，偶对称
h(n)=h(N-1-n)，偶对称
h(n)=h(N−1−n)，偶对称
h ( n ) = − h ( N − 1 − n ) , 奇对称 h(n)=-h(N-1-n),奇对称 h(n)=−h(N−1−n),奇对称
频率响应：
FIR数字滤波器的频率响应为
我们将 H ( e j ω ) H(e^{j\omega}) H(ejω)表示为
其中 ∣ H ( e j ω ) ∣ |H(e^{j\omega})| ∣H(ejω)∣是幅度响应， H ( ω ) H(\omega) H(ω)为幅度函数， θ ( ω ) \theta(\omega) θ(ω)为相位函数。幅度函数是一个纯实数，可为正值和负值。
h(n)为偶对称时，幅度函数和相位函数为
h(n)奇对称时，幅度函数和相位函数为
幅度函数：
h(n)的上述两种情况又可分别分为N为偶数和奇数两种情况，所以一共可以分为四种类型的线性相位FIR滤波器。
-
h(n)偶对称，N为奇数（1型）
-
h(n)偶对称，N为偶数（2型）
-
h(n)奇对称，N为奇数（3型
-
h(n)奇对称，N为偶数（4型）
2.窗函数法设计FIR数字滤波器
窗函数设计法的基本思想为，首先选择一个适当的理想的滤波器 H d ( e j ω ) H_d(e^{j\omega}) Hd(ejω)，然后用窗函数截取它的单位脉冲响应 H d ( n ) H_d(n) Hd(n)，得到线性相位和因果的FIR滤波器。这种方法的重点是选择一个合适的的窗函数和理想滤波器，使设计的滤波器的单位脉冲响应逼近理想滤波器的单位脉冲响应。
设计步骤：
-
给定理想滤波器的频率响应 H d ( e j ω ) H_d(e^{j\omega}) Hd(ejω)，在通带上具有单位增益和线性相位，在阻带上具有零响应。一个带宽为 ω c ( ω c < π ) \omega_c(\omega_c<\pi) ωc(ωc<π)的低通滤波器由下式给定
其中 α \alpha α为采样延迟，其作用是为了得到因果的系统。
-
确定这个滤波器的单位脉冲响应
为了得到一个h(n)长度为N 的因果的线性相位 FIR 滤波器，我们令 α = ( N − 1 ) / 2 \alpha=(N-1)/2 α=(N−1)/2
-
用窗函数截取 h d ( n ) h_d(n) hd(n)得到所设计FIR数字滤波器h(n)
h ( n ) = h d ( n ) ω ( n ) h(n)=h_d(n)\omega(n) h(n)=hd(n)ω(n)
窗函数的选择：
常用的窗函数有矩形（Rectangular）窗、汉宁（Hanning）窗、海明（Hamming）窗、布莱克曼（Blackman）窗、凯瑟（Kaiser）窗等。MATLAB 提供了一些函数用于产生窗函数，如下：
在设计过程中我们需要根据给定的滤波器技术指标，选择滤波器长度N 和窗函数 ω ( n ) \omega(n) ω(n)。
3.频率取样法设计FIR数字滤波器
频率取样法从频域出发，把理想的滤波器 H d ( e j ω ) H_d(e^{j\omega}) Hd(ejω)等间隔采样得到 H d ( k ) H_d(k) Hd(k)，将 H d ( k ) H_d(k) Hd(k)作为实际设计滤波器的H(k)
得到H(k)以后可以由H(k)来唯一确定滤波器的单位脉冲响应h(n)， H ( e j ω ) H(e^{j\omega}) H(ejω)也可以由H(k)求得
其中 ϕ ( x ) \phi(x) ϕ(x)为内插函数
由H(k)求得的频率响应 H ( e j ω ) H(e^{j\omega}) H(ejω)将逼近 H d ( e j ω ) H_d(e^{j\omega}) Hd(ejω) 。
将H(k)表示成如下形式
当h(n)为实数，则 H ( k ) = H ∗ ( N − k ) H(k)=H^*(N-k) H(k)=H∗(N−k)，由此得到
即 H r ( k ) H_r(k) Hr(k)以 k = N / 2 k=N/2 k=N/2为中心呈偶对称。再利用线性相位条件可知，对于1型和2型线性相位滤波器
对于3型和4型线性相位滤波器
其中，$\left\lfloor x \right\rfloor $表示取小于该数的最大的整数。
4.利用MATLAB函数设计FIR数字滤波器
- MATLAB 提供了函数fir1采用窗函数法设计线性相位FIR 数字滤波器
- fir2函数结合频率取样法和窗函数法综合设计具有任意频率响应的FIR数字滤波器
三、实验内容
1.设计一个数字低通FIR滤波器，其技术指标如下
- ω p = 0.2 π , R p = 0.25 d B {\omega}_p = 0.2\pi,R_p=0.25dB ωp=0.2π,Rp=0.25dB
- ω s t = 0.3 π , A s = 50 d B {\omega}_{st} = 0.3\pi,A_s=50dB ωst=0.3π,As=50dB
分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。
矩形窗
clc;
clear;
close all
%矩形窗
wp = 0.2*pi;
wst = 0.3*pi;
tr_width = wst-wp;
%窗函数长度
N = ceil(1.8*pi/tr_width)+1;
n = 0:N-1;
wc = (wp+wst)/2;%中频
alpha = (N-1)/2;%离散频率向量的第一个频率索引
hd = (wc/pi)*sinc((wc/pi)*(n-alpha));
%产生窗函数，对hd(n)加窗得到h(n)
w_boxcar = boxcar(N)';
h = hd.*w_boxcar;
%绘图
subplot(221);
stem(n,hd,'filled');
axis tight;
xlabel('n');
ylabel('hd(n)');
[Hr,w1] = zerophase(h);%窗函数零点
subplot(222);
plot(w1/pi, Hr);
axis tight;
xlabel('\omega/\pi');
ylabel('H(\omega)');
subplot(223);
stem(n,h,'filled');
axis tight;
xlabel('n');
ylabel('h(n)');
[H,w]= freqz(h,1);
subplot(224);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
xlabel('\omega/\pi');
ylabel('dB');
grid on;
读图可知，阻带衰减小于50dB，因此不符合指标。
汉宁窗
clc;
clear;
close all
%汉宁窗
wp = 0.2*pi;
wst = 0.3*pi;
tr_width = wst-wp;
%窗函数长度
N = ceil(6.2*pi/tr_width)+1;
n = 0:N-1;
wc = (wp+wst)/2;%中频
alpha = (N-1)/2;%离散频率向量的第一个频率索引
hd = (wc/pi)*sinc((wc/pi)*(n-alpha));
%产生窗函数，对hd(n)加窗得到h(n)
w_hanning = hanning(N)';
h = hd.*w_hanning;
%绘图
subplot(221);
stem(n,hd,'filled');
axis tight;
xlabel('n');
ylabel('hd(n)');
[Hr,w1] = zerophase(h);%窗函数零点
subplot(222);
plot(w1/pi, Hr);
axis tight;
xlabel('\omega/\pi');
ylabel('H(\omega)');
subplot(223);
stem(n,h,'filled');
axis tight;
xlabel('n');
ylabel('h(n)');
[H,w]= freqz(h,1);
subplot(224);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
xlabel('\omega/\pi');
ylabel('dB');
grid on;
读图可知，阻带衰减小于50dB，因此不符合要求。
海明窗
clc;
clear;
close all
%海明窗
wp = 0.2*pi;
wst = 0.3*pi;
tr_width = wst-wp;
%窗函数长度
N = ceil(6.6*pi/tr_width)+1;
n = 0:N-1;
wc = (wp+wst)/2;%中频
alpha = (N-1)/2;%离散频率向量的第一个频率索引
hd = (wc/pi)*sinc((wc/pi)*(n-alpha));
%产生窗函数，对hd(n)加窗得到h(n)
w_hamming = hamming(N)';
h = hd.*w_hamming;
%绘图
subplot(221);
stem(n,hd,'filled');
axis tight;
xlabel('n');
ylabel('hd(n)');
[Hr,w1] = zerophase(h);%窗函数零点
subplot(222);
plot(w1/pi, Hr);
axis tight;
xlabel('\omega/\pi');
ylabel('H(\omega)');
subplot(223);
stem(n,h,'filled');
axis tight;
xlabel('n');
ylabel('h(n)');
[H,w]= freqz(h,1);
subplot(224);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
xlabel('\omega/\pi');
ylabel('dB');
grid on;
读图可知，阻带衰减大于50dB，因此符合要求。
布莱克曼窗
clc;
clear;
close all
%布莱克曼窗
wp = 0.2*pi;
wst = 0.3*pi;
tr_width = wst-wp;
%窗函数长度
N = ceil(11*pi/tr_width)+1;
n = 0:N-1;
wc = (wp+wst)/2;%中频
alpha = (N-1)/2;%离散频率向量的第一个频率索引
hd = (wc/pi)*sinc((wc/pi)*(n-alpha));
%产生窗函数，对hd(n)加窗得到h(n)
w_blackman = blackman(N)';
h = hd.*w_blackman;
%绘图
subplot(221);
stem(n,hd,'filled');
axis tight;
xlabel('n');
ylabel('hd(n)');
[Hr,w1] = zerophase(h);%窗函数零点
subplot(222);
plot(w1/pi, Hr);
axis tight;
xlabel('\omega/\pi');
ylabel('H(\omega)');
subplot(223);
stem(n,h,'filled');
axis tight;
xlabel('n');
ylabel('h(n)');
[H,w]= freqz(h,1);
subplot(224);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
xlabel('\omega/\pi');
ylabel('dB');
grid on;
读图可知，阻带衰减大于50dB，因此符合要求。
凯瑟窗
clc;
clear;
close all
%凯瑟窗
wp = 0.2*pi;
wst = 0.3*pi;
tr_width = wst-wp;
As = 50;
%窗函数长度
N = ceil((As-7.95)/(2.285*tr_width))+1;
beta = 0.1102*(As-8.7);
n = 0:N-1;
wc = (wp+wst)/2;%中频
alpha = (N-1)/2;%离散频率向量的第一个频率索引
hd = (wc/pi)*sinc((wc/pi)*(n-alpha));
%产生窗函数，对hd(n)加窗得到h(n)
w_kaiser = kaiser(N,beta)';
h = hd.*w_kaiser;
%绘图
subplot(221);
stem(n,hd,'filled');
axis tight;
xlabel('n');
ylabel('hd(n)');
[Hr,w1] = zerophase(h);%窗函数零点
subplot(222);
plot(w1/pi, Hr);
axis tight;
xlabel('\omega/\pi');
ylabel('H(\omega)');
subplot(223);
stem(n,h,'filled');
axis tight;
xlabel('n');
ylabel('h(n)');
[H,w]= freqz(h,1);
subplot(224);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
xlabel('\omega/\pi');
ylabel('dB');
grid on;
读图可知，阻带衰减大于50dB，因此符合要求。
2.设计一个数字带通FIR滤波器，其技术指标如下
- 下阻带边缘： ω s t 1 = 0.2 π , A s = 60 d B \omega_{st1}=0.2\pi,A_s=60dB ωst1=0.2π,As=60dB
- 下通带边缘： ω p 1 = 0.35 π , R p = 1 d B \omega_{p1}=0.35\pi,R_p=1dB ωp1=0.35π,Rp=1dB
- 上通带边缘： ω p 2 = 0.65 π , R p = 1 d B \omega_{p2}=0.65\pi,R_p=1dB ωp2=0.65π,Rp=1dB
- 上阻带边缘： ω s t 2 = 0.8 π , A s = 60 d B \omega_{st2}=0.8\pi,A_s=60dB ωst2=0.8π,As=60dB
clc;
clear;
close all
%参数设置
ws1 = 0.2*pi;
ws2 = 0.8*pi;
wp1 = 0.35*pi;
wp2 = 0.65*pi;
tr_width = wp1-ws1;
wc1 = (ws1+wp1)/2;
wc2 = (ws2+wp2)/2;
N = ceil(11*pi/tr_width)+1;
n = 0:N-1;
alpha =(N-1)/2;
hd = (wc2/pi)*sinc((wc2/pi)*(n-alpha))-(wc1/pi)*sinc((wc1/pi)*(n-alpha));
w_blackman = blackman(N)';%布莱克曼窗
h = hd.*w_blackman;
%绘图
subplot(221);
stem(n,hd,'filled');
axis tight;
xlabel('n');
ylabel('hd(n)');
[Hr,w1] = zerophase(h);
subplot(222);
plot(w1/pi,Hr);
axis tight;
xlabel('\omega/\pi');
ylabel('H(\omega)');
subplot(223);
stem(n,h,'filled');
axis tight;
xlabel('n');
ylabel('h(n)');
[H,w] = freqz(h,1);
subplot(224);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
xlabel('\omega/\pi');
ylabel('dB');
grid on;
读图可知，符合设计要求。
3.采取频率取样法设计FIR数字低通滤波器，满足以下指标
- ω p = 0.2 π , R p = 0.25 d B {\omega}_p = 0.2\pi,R_p=0.25dB ωp=0.2π,Rp=0.25dB
- ω s t = 0.3 π , A s = 50 d B {\omega}_{st} = 0.3\pi,A_s=50dB ωst=0.3π,As=50dB
（1）取N=20，过渡带没有样本。
clc;
clear;
close all
%参数设置
N = 20;
alpha = (N-1)/2;%希尔伯特变换参数
l= 0:N-1;%频率点数
wl = (2*pi/N)*l;%频率向量
%希尔伯特变换矩阵
Hrs = [ones(1,3),zeros(1,15),ones(1,2)];
Hdr = [1,1,0,0];
%权重向量
wdl = [0,0.25,0.25,1];
%希尔伯特变换索引
k1 = 0:floor((N-1)/2);
k2 = floor((N-1)/2)+1:N-1;
%希尔伯特变换角
angH = [-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)];
H = Hrs.*exp(1i* angH);
%对希尔伯特变换矩阵进行快速傅里叶变换
h = ifft(H,N);
%计算频率响应
w = [0:500]*pi/500;
H = freqz(h,1,w);
[Hr,wr] = zerophase(h);
%绘图
subplot(221);
plot(wdl,Hdr,wl(1:11)/pi,Hrs(1:11),'o');
axis([0, 1,-0.1, 1.1]);
xlabel('\omega(\pi)');
ylabel('Hr(k)');
subplot(222);
stem(l,h,'filled');
axis([0,N-1,-0.1,0.3]);
xlabel('n');
ylabel('h(n)');
subplot(223);
plot(wr/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o');
axis([0, 1,-0.2, 1.2]);
xlabel('\omega(\pi)');
ylabel('Hr(\omega)');
subplot(224);
plot(w/pi,20*log10((abs(H)/max(abs(H)))));
axis([0,1,-50,5]);
grid;
xlabel('\omega(\pi)');
ylabel('dB');
（2）取N=40，过渡带有一个样本，T=0.39。
clc;
clear;
close all
%参数设置
N = 40;
alpha = (N-1)/2;%希尔伯特变换参数
l= 0:N-1;%频率点数
wl = (2*pi/N)*l;%频率向量
%希尔伯特变换矩阵
Hrs = [ones(1,5),0.39,zeros(1,29),0.39,ones(1,4)];
Hdr = [1,1,0.39,0,0];
%权重向量
wdl = [0,0.2,0.25,0.3,1];
%希尔伯特变换索引
k1 = 0:floor((N-1)/2);
k2 = floor((N-1)/2)+1:N-1;
%希尔伯特变换角
angH = [-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)];
H = Hrs.*exp(1i* angH);
%对希尔伯特变换矩阵进行快速傅里叶变换
h = ifft(H,N);
%计算频率响应
w = [0:500]*pi/500;
H = freqz(h,1,w);
[Hr,wr] = zerophase(h);
%绘图
subplot(221);
plot(wdl,Hdr,wl(1:21)/pi,Hrs(1:21),'o');
axis([0, 1,-0.1, 1.1]);
xlabel('\omega(\pi)');
ylabel('Hr(k)');
subplot(222);
stem(l,h,'filled');
axis([0,N-1,-0.1,0.3]);
xlabel('n');
ylabel('h(n)');
subplot(223);
plot(wr/pi,Hr,wl(1:21)/pi,Hrs(1:21),'o');
axis([0, 1,-0.2, 1.2]);
xlabel('\omega(\pi)');
ylabel('Hr(\omega)');
subplot(224);
plot(w/pi,20*log10((abs(H)/max(abs(H)))));
axis([0,1,-50,5]);
grid;
xlabel('\omega(\pi)');
ylabel('dB');
（3）取N=60，过渡带有两个样本，T1=0.5925，T2=0.1099。
clc;
clear;
close all
%参数设置
N = 60;
alpha = (N-1)/2;%希尔伯特变换参数
l= 0:N-1;%频率点数
wl = (2*pi/N)*l;%频率向量
%希尔伯特变换矩阵
Hrs = [ones(1,7),0.5925,0.1099,zeros(1,43),0.1099,0.5925,ones(1,6)];
Hdr = [1,1,0.5925,0.1099,0,0];
%权重向量
wdl = [0,0.2,7/30,8/30,0.3,1];
%希尔伯特变换索引
k1 = 0:floor((N-1)/2);
k2 = floor((N-1)/2)+1:N-1;
%希尔伯特变换角
angH = [-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)];
H = Hrs.*exp(1i* angH);
%对希尔伯特变换矩阵进行快速傅里叶变换
h = ifft(H,N);
%计算频率响应
w = [0:500]*pi/500;
H = freqz(h,1,w);
Hr = abs(H);
%绘图
subplot(221);
plot(wdl,Hdr,wl(1:31)/pi,Hrs(1:31),'o');
axis([0, 1,-0.1, 1.1]);
xlabel('\omega(\pi)');
ylabel('Hr(k)');
subplot(222);
stem(l,h,'filled');
axis([0,N-1,-0.1,0.3]);
xlabel('n');
ylabel('h(n)');
subplot(223);
plot(w/pi,Hr,wl(1:31)/pi,Hrs(1:31),'o');
axis([0, 1,-0.2, 1.2]);
xlabel('\omega(\pi)');
ylabel('Hr(\omega)');
subplot(224);
plot(w/pi,20*log10((abs(H)/max(abs(H)))));
axis([0,1,-50,5]);
grid;
xlabel('\omega(\pi)');
ylabel('dB');
（4）分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。
此数字滤波器阻带衰减、通带波动要求均满足。
第一个由实验结果可以看出滤波器通带波动约为1.09dB和阻带衰减约为16.2dB，阻带衰减明显不满足滤波器设计要求而通带波动也远达不到要求。需要进一步增大取样点数N，并在过渡带内取样。
第二个由实验结果可以看出滤波器通带波动约为0.306dB和阻带衰减约为63.5dB，阻带满足滤波器设计要求而通带波动不满足要求。第三个数字滤波器阻带衰减、通带波动要求均满足。
第三个数字滤波器阻带衰减、通带波动均满足要求。
4.采用频率取样技术设计下面的高通滤波器
- ω p = 0.8 π , R p = 1 d B {\omega}_p = 0.8\pi,R_p=1dB ωp=0.8π,Rp=1dB
- ω s t = 0.6 π , A s = 50 d B {\omega}_{st} = 0.6\pi,A_s=50dB ωst=0.6π,As=50dB
对于高通滤波器，N必须为奇数（或1型滤波器）。选择N=33，过渡带有两个样本，过渡带样本最优值为T1=0.1095，T2=0.598。
clc;
clear;
close all
%参数设置
N = 33;
alpha = (N-1)/2;%希尔伯特变换参数
l= 0:N-1;%频率点数
wl = (2*pi/N)*l;%频率向量
%希尔伯特变换矩阵
Hrs = [zeros(1,11),0.1095,0.598,ones(1,7),0.598,0.1095,zeros(1,11)];
Hdr = [0,0,0.1095,0.598,1,1];
%权重向量
wdl = [0,0.6,2/3,24/33,0.8,1];
%希尔伯特变换索引
k1 = 0:floor((N-1)/2);
k2 = floor((N-1)/2)+1:N-1;
%希尔伯特变换角
angH = [-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N-k2)];
H = Hrs.*exp(1i* angH);
%对希尔伯特变换矩阵进行快速傅里叶变换
h = ifft(H,N);
%计算频率响应
w = [0:500]*pi/500;
Hr = h*cos((alpha-l)'*w);
%绘图
subplot(221);
plot(wdl,Hdr,wl(1:17)/pi,Hrs(1:17),'o');
axis([0, 1,-0.1, 1.1]);
xlabel('\omega(\pi)');
ylabel('Hr(k)');
subplot(222);
stem(l,h,'filled');
axis([0,N-1,-0.1,0.3]);
xlabel('n');
ylabel('h(n)');
subplot(223);
plot(w/pi,Hr,wl(1:17)/pi,Hrs(1:17),'o');
axis([0, 1,-0.2, 1.2]);
xlabel('\omega(\pi)');
ylabel('Hr(\omega)');
subplot(224);
plot(w/pi,20*log10((abs(Hr)/max(abs(Hr)))));
axis([0,1,-50,5]);
grid;
xlabel('\omega(\pi)');
ylabel('dB');
读图可知，该滤波器满足阻带衰减大于50dB。
四、实验心得体会
本次实验要求我们掌握频率取样法设计FIR数字滤波器的原理及具体方法。
通过查阅相关资料，可以知道利用频率取样法设计FIR数字滤波器设计的基本原理是:从频域出发，把理想的滤波器
H
d
(
e
j
ω
)
H_d(e^{j\omega})
Hd(ejω)等间隔采样得到
H
d
(
k
)
H_d(k)
Hd(k)，将
H
d
(
k
)
H_d(k)
Hd(k)作为实际设计滤波器的
H
(
k
)
H(k)
H(k)，然后通过内插公式恢复出实际的滤波器的频响。
通过此次实验，我充分了解了频率取样法在设计FIR滤波器时的原理和使用方法对FIR滤波器的性能加深了理解，并且对设计滤波器的过程更加熟悉了，这对我以后的学习和将来的工作都大有裨益。
扫一扫
5万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
