北理工数字信号处理实验报告.docx

本科实验报告
实验名称： 数字信号处理实验

课程名称：
数字信号处理实验
实验时间：
周二8~10节
任课教师：
范哲意
实验地点：
4—423
实验教师：
范哲意
实验类型：
□ 原理验证
□ 综合设计
□ 自主创新
学生姓名：
任科飞
学号/班级：
1120121159/05111251
组 号：

学 院：
信息与电子
同组搭档：

专 业：
信息工程
成 绩：

实验1 利用DFT分析信号频谱
一、实验目的
1.加深对DFT原理的理解。
2.应用DFT分析信号频谱。
3.深刻理解利用DFT分析信号频谱的原理，分析现实过程现象及解决办法。

二、实验原理
1、DFT和DTFT的关系
有限长序列的离散时间傅里叶变换在频率区间的N个等分点上的N个取样值可以由下式表示：

由上式可知，序列的N点DFT，实际上就是序列的DTFT在N个等间隔频率点上样本。
2、利用DFT求DTFT
方法1：由恢复出的方法如图2.1所示：

图 2.1.由 N点DFT恢复频谱DTFT的流程
由图2.1所示流程图可知：

由式2-2可以得到

其中为内插函数

方法2：然而在实际MATLAB计算中，上诉插值公式不见得是最好的方法。由于DFT是DTFT的取样值，其相邻的两个频率样本点的间距为,所以如果我们增加数据的长度N，使得得到的DFT谱线就更加精细，其包络就越接近DTFT的结果，这样可以利用DFT来近似计算DTFT。如果没有更多的数据，可以通过补零来增加数据长度。
3、利用DFT分析连续时间信号的频谱
采用计算机分析连续时间信号的频谱，第一步就是把连续时间信号离散化，这里需要进行连个操作：一是采样，二是截断。
对于连续非周期信号，按采样间隔T进行采样，截取长度为M，那么

对进行N点的频率采样，得到

因此，可以将利用DFT分析连续非周期信号频谱的步骤归纳如下：
（1）确定时域采样间隔T，得到离散序列；
（2）确定截取长度M，得到M点离散序列，这里的为窗函数。
（3）确定频域采样点数N，要求。
（4）利用FFT计算离散序列的N点DFT，得到。
（5）根据式（2-6）由计算采样点的近似值。
采用上诉方法计算的频谱，需要注意如下三点问题：
（1）频谱混叠。如果不满足采样定理的条件，频谱会很出现混叠误差。对于频谱无限宽的信号，应考虑覆盖大部分主要频率的范围。
（2）栅栏效应和频谱分辨率。使用DFT计算频谱，得到的结果只是N个频谱样本值，样本值之间的频谱是未知的，就像通过一个栅栏观察频谱，称为“栅栏效应”。频谱分辨率与记录长度成正比，提高频谱分辨率，就要增加记录时间。
（3）频谱泄露。对于信号截断会把窗函数的频谱会引入到信号频谱中，造成频谱泄露。解决这问题的主要办法是采用旁瓣小的窗函数，频谱泄露和窗函数均会引起误差。
因此，要合理选取采样间隔和截取长度，必要时还需考虑适当的窗。
对于连续周期信号，我们在采用计算机进行计算时，也总是要进行截断，序列总是有限长的，仍然可以采用上诉方法近似计算。
4、可能用到MATLAB函数与代码
实验中的DFT运算可以采用MATLAB中提供的FFT来实现。
DTFT可以利用MATLAB矩阵运算的方法进行计算。
三、实验内容
1. ，完成如下要求：
（1）计算其DTFT，并画出区间的波形。
（2）计算4 点DFT，并把结果显示在（1）所画的图形中。
（3）对补零，计算64 点DFT，并显示结果。
（4）是否可以由DFT 计算DTFT，如果可以，请编程实现。
2. 考察序列

（1）时，用DFT 估计的频谱；将补零加长到长度为100点序列用DFT
估计的频谱。要求画出相应波形。
（2）时，用DFT 估计x(n)的频谱，并画出波形。
3. 已知信号 ，其中，，。从的表达式可以看出，它包含三个频率的正弦波，但是，从其时域波形来看，似乎是一个正弦信号，利用DFT做频谱分析，确定适合的参数，使得到的频谱的频率分辨率符合需要。
4.利用DFT近似分析连续时间信号xt=e-0.1tu(t)的频谱（幅度谱）。分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定合适的参数。

四、实验代码及实验结果
第一题
实验代码
x = [2, -1, 1, 1];
n = 0 : 3;
n1 = 0 : 63;
z = zeros(1,60);
x1 = [x z];
y = fft(x);
y1 = fft(x1);
w = -pi : 0.01 * pi : pi;
X = x * exp(-j * n' * w);
subplot(211);
plot(w, abs(X));
xlabel('\Omega/\pi');
title('Magnitude');
hold on;
stem(n, y, 'filled');
axis tight;
subplot(212);
plot(w, angle(X)/pi);
xlabel('\Omega/\pi');
title('Phase');
axis tight;
figure;
stem(n1, y1, 'filled');

实验结果

分析
可以由DFT计算DTFT。增加补零的个数，提高采样密度，可以由DFT近似计算DTFT。

第二题
实验代码
n = 0 : 10;
x = cos(0.48 * pi * n) + cos(0.52 * pi * n);
n1 = 0 : 99;
n2 = 0 : 100;
z = zeros(1, 89);
x1 = [x z];
x2 = cos(0.48 * pi * n2) + cos(0.52 * pi * n2);
y = fft(x);
y1 = fft(x1);
y2 = fft(x2);
subplot(211);
stem(n, y, 'filled');
title('0<=n<=10');
subplot(212);
stem(n1, y1, 'filled');
title('0<=n<99');
figure;
stem(n2, y2, 'filled');
title('0<=n<=100');

实验结果

(3)分析
可以通过增大截取长度和增加补零的个数来提高频谱分辨率，但是补零不能够增加分辨力。

第三题
实验代码
>> n=0:39;
>> x1=0.15*sin(2*pi*n)+sin(4*pi*n)-0.1*sin(6*pi*n);
>> X1=fft(x1,40);
>> subplot(311)
>> stem(n,abs(X1));
>> x2=0.15*sin(2*pi*(1/6)*n)+sin(4*pi*(1/6)*n)-0.1*sin(6*pi*(1/6)*n);
>> X2=fft(x2,40);
>> subplot(312)
>> stem(n,abs(X2));
>> x3=0.15*sin(2*pi*0.1*n)+sin(4*pi*0.1*n)-0.1*sin(6*pi*0.1*n);
>> X3=fft(x3,40);
>> subplot(313)
>> stem(n,abs(X3));
实验结果

分析：减小采样周期，可以使频谱分辨率提高。

第四题
实验代码
n=0:100;
x=exp(-0.1*n);
X1=fft(x);
k=0:100;
stem(k,abs(X1));
xlabel('k');
实验结果

实验代码
n=0:2:100;
x=exp(-0.1*n);
X2=fft(x);
k=0:50;
stem(k,abs(X2));
xlabel('k');
实验结果

实验代码
>> n=0:50;
x=exp(-0.1*n);
X3=fft(x);
k=0:50;
stem(k,abs(X3));
xlabel('k');
实验结果

实验代码
>> n=0:2:50;
x=exp(-0.1*n);
X4=fft(x);
k=0:25;
stem(k,abs(X4));
xlabel('k');
实验结果

实验代码
>> n=0:150;
x=exp(-0.1*n);
X5=fft(x);
k=0:150;
stem(k,abs(X5));
xlabel('k');
实验结果

实验代码
>> n=0:2:150;
x=exp(-0.1*n);
X5=fft(x);
k=0:75;
stem(k,abs(X5));xlabel('k');
实验结果

分析：采样区间可以为0到100，采样间隔为2。

五、心得与体会
通过本次实验，我们掌握并加深对DFT原理的理解并且学会应用DFT分析信号频谱，在此基础上利用DFT分析信号频谱的原理，掌握了分析现实现象的步骤及办法。
并且，通过这次的实验，我对DFT的原理确实有了更深的了解，对DFT和DTFT之间的关系也有了更多的认识：序列x(n)的N点DFT X(K)，实际上就是x(n)序列的DTFT在N个等间隔频率点上的样本X(K)。所以，我们可以通过增加数据的长度N或通过补零来使DFT更加接近DTFT的结果。这样就可以利用DFT近似计算DTFT。

实验2 利用FFT计算线性卷积（选作）
一、实验目的
1.掌握利用FFT计算线性卷积的原理及具体实现方法。
2.加深理解重叠相加法和重叠保留法。
3.考察利用FFT计算线性卷积各种方法的适用范围。
二、实验原理
1.线性卷积与圆周卷积
设x(n)为L点序列，h(n)为M点序列，x(n)和h(n)的线性卷积为
（3-1）
的长度为L+M-1
x(n)和h(n)的圆周卷积为
（3-2）

圆周卷积与线性卷积相等而不产生交叠的必要条件为
N≥L+M+1 (3-3)
圆周卷积定理：根据DFT性质，x(n)和h(n)的N点圆周卷积的DFT等于它们的DFT的乘积：
(3-4)

2.快速卷积
快速卷积发运用圆周卷积实现线性卷积，根据圆周卷积定理利用FFT算法实现圆周卷积。可将快速卷积运算的步骤归纳如下：
(1)必须选择；为了能使用基-2算法，要求。采用补零的办法使得x(n)和h(n)的长度均为N。
(2)计算x(n)和h(n)的N点FFT。

(3)组成乘积

(4)利用IFFT计算Y(k)的IDFT，得到线性卷积y(n)

3.分段卷积
我们考察单位取样响应为h(n)的线性系统，输入为x(n)，输出为y(n)，则
yn=xn*h(n)
当输入序列x(n)极长时，如果要等x(n)全部集齐时再开始进行卷积，会使输出有较大延时；如果序列太长，需要大量存储单元。为此，我们把x(n)分段，为别求出每段的卷积，合在一起得到最后的总输出。这称为分段卷积。分段卷积可以细分为重叠保留法和重叠相加法。
重叠保留法：设x(n)的长度为，h(n)的长度为M。把序列x(n)分成多段N点序列，每段雨前一段重写M-1个样本。并在第一个输入段前面补M-1个零。计算每一段与h(n)的圆周卷积，其结果中前M-1个不等与线性卷积，应当舍去，只保留后面N-M+1个正确的输出样本，把它们合起来得到总的输出。
利用FFT实现重叠保留法的步骤如下：
(1)在x(n)前面填充M-1个零，扩大以后的序列为

(2)将x(n)分为若干段N点子段，设L=N-M+1为每一段的有效长度，则第i段的数据为：

(3)计算每一段与h(n)的N点圆周卷积，利用FFT计算圆周卷积

(4)舍去每一段卷积结果的前M-1个样本，连接剩下的样本得到卷积结果y(n)。
重叠相加法：设h(n)长度为M，将信号x(n)分解成长为L的子段。以表示没断信号，则：

每一段卷积的长度为L+M-1，所以在做求和时，相邻两段序列由M-1个样本重叠，即前一段的最后M-1个样本和下一段前M-1个样本序列重叠，这个重叠部分相加，再与不重叠的部分共同组成y(n)。
利用FFT实现重叠保留法的步骤如下：
(1)将x(n)分为若干L点子段。
(2)计算每一段与h(n)的卷积，根据快速卷积法利用FFT计算卷积。
(3)将各段相加，得到输出y(n)。

4、可能得到的MATLAB函数
实验中FFT运算可采用MATLAB中提供的函数fft来实现。

三、实验内容
假设要计算序列x(n)=u(n)-u(n-L),0≤n≤L和h(n)=cos(0.2πn),0≤n≤M的线性卷积完成以下实验内容。
1.设L=M,根据线性卷积的表达式和快速卷积的原理分别编程实现计算两个序列线性卷积的方法,比较当序列长度分别为8,16,32,64,256,512,1024时两种方法计算线性卷积所需时间。
2．当L=2048且M=256时比较直接计算线性卷积和快速卷积所需的时间,进一步考察当L=4096且M=256时两种算法所需的时间。
3. 编程实现利用重叠相加法计算两个序列的线性卷积,考察L=2048且M=256时计算线性卷积的时间,与2题的结果进行比较。
4. 编程实现利用重叠保留法计算两个序列的线性卷积,考察L=2048且M=256时计算线性卷积的时间,与2题的结果进行比较。

四、实验代码及实验结果
第一题
实验代码
L = input('length=');
n1 = 0 : 2 * L - 1;
n2 = 0 : L - 1;
x = ones(1, L);
h = cos(0.2 * pi * n2);
z = zeros(1, L - 1);
x1 = [x z];
h1 = [h z];
tic;
y1 = conv(x, h);
toc;
tic;
X = fft(x1);
H = fft(h1);
Y = X .* H;
y = ifft(Y);
toc;
实验结果（第一行为线性卷积，第二行为FFT）
L=8
M=8
时间已过 0.232076 秒。
时间已过 0.087391 秒。
L=16
M=16
时间已过 0.000332 秒。
时间已过 0.000098 秒。
L=32
M=32
时间已过 0.000150 秒。
时间已过 0.000051 秒。
L=64
M=64
时间已过 0.289400 秒。
时间已过 0.000758 秒。
L=256
M=256
时间已过 0.000238 秒。
时间已过 0.000169 秒。
L=512
M=512
时间已过 0.000455 秒。
时间已过 0.000235 秒。
L=1024
M=1024
时间已过 0.001132 秒。
时间已过 0.000922 秒。
分析：可见在在相同长度下，快速卷积比线性卷积快差不多一倍的时间。

第二题
实验代码

实验结果：（第一行为线性卷积，第二行为FFT）
L=2048
M=256
时间已过0.012061秒。
时间已过0.000581秒。
分析：快速卷积比线性卷积快。

第三题
Lx=input('L=');
M=input('M=');
N=M+1;
x=1;
y=1;
for i=1:L-1
x=[1 x];
y=[y cos(0.2*pi*i)];
end
x=[x,zeros(1,N-1)];
t=zeros(1,M-1);
Y=zeros(1,Lx+M-1);
a=floor(Lx/N);
tic
for k=0:a
A=x(k*N+1:k*N+N);
y1=fft(A,Lx+M-1);
y2=fft(h,Lx+M-1);
y3=y1.*y2;
q=ifft(y3,Lx+M-1);
Y(k*N+1:k*N+M-1)=q(1:M-1)+t(1:M-1);
Y(k*N+M:k*N+N)=q(M:N);
t(1:M-1)=q(N+1:N+M-1);
end
toc
实验结果：
L=2048
M=256
时间已过0.007656秒。
分析：较第2题结果快。

第四题
实验代码：
Lx=input('L=');
M=input('M=');
N=M+1;
x=1;
y=1;
for i=1:L-1
x=[1 x];
y=[y cos(0.2*pi*i)];
end
M1=M-1;
L=N+M1;
h=[h,zeros(1,N-M)];
x=[zeros(1,M1),x,zeros(1,N-1)];
a=floor ((Lx+M1-1)/(L))+1;
Y=zeros(1,N);
tic
for k=0:a-1
xk=x(k*L+1:k*L+N);
b=fft(xk,N);
C=fft(h,N);
Z=b.*C;
Y(k+1,:)=ifft(Z,N);
end
toc
L=2048
M=256
实验结果：
L=2048
M=256
时间已过0.002848秒。
分析：较第2题结果快。

五、心得与体会
本次实验要求我们掌握利用FFT计算线性卷积的原理及具体实现方法，通过实验加深理解重叠相加法和重叠保留法并考察利用FFT计算线性卷积各种方法的适用范围。
本次实验让我切实看到了FFT算法的高效性及对于庞大的数据的处理实时性，同时对重叠保留法、重叠相加法也有了更深的认识，对以后的实验有了很好的理论基础，受益颇多。

实验3 IIR数字滤波器设计
一、实验目的
1.掌握利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方法。
2.加深理解数字滤波器和模拟滤波器之间的技术指标转化。
3.掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及适用范围。

二、实验内容
1.设采样频率为fs =4kHz，采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3dB截止频率为fc=1kHz。
2.设采样频率为fs=10kHz，设计数字低通滤波器，满足如下指标
通带截止频率：fp=1kHz，通带波动：Rp=1dB
阻带截止频率：fst=1.5kHz，阻带衰减：As=15dB
要求分别采用巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。结合实验结果，分别讨论采用上述设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。

四、实验代码及实验结果
实验代码：
巴特沃斯：
fp = 1000;
fs = 1500;
f = 10000;
Wp = 2 * pi * fp;
Ws = 2 * pi * fs;
wp = 2 * f * tan(2 * pi * fp / (2 * f));
ws = 2 * f * tan(2 * pi * fs / (2 * f));
Rp = 1;
As = 15;
N1 = ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws)));
N2 = ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(wp/ws)));
OmegaC1 = Wp / ((10 ^ (Rp / 10) - 1) ^ (1 / (2 * N1)));
OmegaC2 = wp / ((10 ^ (Rp / 10) - 1) ^ (1 / (2 * N2)));
[z, p, k] = buttap(N1);
b1 = k * OmegaC1 ^ N1;
a1 = poly(p * OmegaC1);
[H, w] = freqs(b1, a1);
subplot(221);
plot(w/pi, abs(H));
grid on;
xlabel('\Omega(\pi)');ylabel('|H(j\Omega)|');title('°ÍÌØÎÖ˹Âö³åÂö³åÏìÓ¦²»±ä');
subplot(223);
plot(w/pi, angle(H));
grid on;
xlabel('\Omega(\pi)');ylabel('Phase of H(j\Omega)|');title('°ÍÌØÎÖ˹Âö³åÂö³åÏìÓ¦²»±ä');
[b2, a2] = butter(N2, OmegaC2, 's');
[H, w] = freqs(b2, a2);
subplot(222);
plot(w/pi, abs(H));
axis tight;
grid on;
xlabel('\Omega(\pi)');ylabel('|H(j\Omega)|');title('°ÍÌØÎÖ˹˫ÏßÐԱ任');
subplot(224);
plot(w/pi, angle(H));
grid on;
xlabel('\Omega(\pi)');ylabel('Phase of H(j\Omega)|');title('°ÍÌØÎÖ˹˫ÏßÐԱ任');
实验结果：

切比雪夫Ⅰ型：
fp = 1000;
fs = 1500;
f = 10000;
Wp = 2 * pi * fp;
Ws = 2 * pi * fs;
wp = 2 * f * tan(2 * pi * fp / (2 * f));
ws = 2 * f * tan(2 * pi * fs / (2 * f));
Rp = 1;
As = 15;
e = (10 ^ (Rp / 10) - 1) ^ 0.5;
N1 = ceil(acos(((As * As) - 1) ^ 0.5 / e) / acos(Ws / Wp));
N2 = ceil(acos(((As * As) - 1) ^ 0.5 / e) / acos(ws / wp));
[b1, a1] = cheby1(N1, Rp, Ws, 's');
[H, w] = freqs(b1, a1);
subplot(221);
plot(w/pi, abs(H));
grid on;
xlabel('\Omega(\pi)');ylabel('|H(j\Omega)|');title('ÇбÈÑ©·ò1Âö³åÂö³åÏìÓ¦²»±ä');
subplot(223);
plot(w/pi, angle(H));
grid on;
xlabel('\Omega(\pi)');ylabel('Phase of H(j\Omega)|');title('ÇбÈÑ©·ò1Âö³åÂö³åÏìÓ¦²»±ä');
[b2, a2] = cheby1(N2, Rp, ws, 's');
[H, w] = freqs(b2, a2);
subplot(222);
plot(w/pi, abs(H));
grid on;
xlabel('\Omega(\pi)');ylabel('|H(j\Omega)|');title('ÇбÈÑ©·ò1Ë«ÏßÐԱ任');
subplot(224);
plot(w/pi, angle(H));
grid on;
xlabel('\Omega(\pi)');ylabel('Phase of H(j\Omega)|');title('ÇбÈÑ©·ò1Ë«ÏßÐԱ任');
实验结果：

切比雪夫Ⅱ型：
fp = 1000;
fs = 1500;
f = 10000;
Wp = 2 * pi * fp;
Ws = 2 * pi * fs;
wp = 2 * f * tan(2 * pi * fp / (2 * f));
ws = 2 * f * tan(2 * pi * fs / (2 * f));
Rp = 1;
As = 15;
e = (10 ^ (Rp / 10) - 1) ^ 0.5;
K = Wp / Ws;
k1 = e / ((As * As) - 1) ^ 0.5;
k = wp / ws;
N1 = ceil(ellipke(K) * ellipke((1 - k1) ^ 0.5) / (ellipke(k1) * ellipke((1 - K) ^ 0.5)));
N2 = ceil(ellipke(k) * ellipke((1 - k1) ^ 0.5) / (ellipke(k1) * ellipke((1 - k) ^ 0.5)));
[b1, a1] = ellip(N1, Rp, As, Ws, 's');
[H, w] = freqs(b1, a1);
subplot(221);
plot(w/pi, abs(H));
grid on;
xlabel('\Omega(\pi)');ylabel('|H(j\Omega)|');title('ÍÖÔ²Âö³åÂö³åÏìÓ¦²»±ä');
subplot(223);
plot(w/pi, angle(H));
grid on;
xlabel('\Omega(\pi)');ylabel('Phase of H(j\Omega)|');title('ÍÖÔ²Âö³åÂö³åÏìÓ¦²»±ä');
[b2, a2] = ellip(N2, Rp, As, ws, 's');
[H, w] = freqs(b2, a2);
subplot(222);
plot(w/pi, abs(H));
grid on;
xlabel('\Omega(\pi)');ylabel('|H(j\Omega)|');title('ÍÖԲ˫ÏßÐԱ任');
subplot(224);
plot(w/pi, angle(H));
grid on;
xlabel('\Omega(\pi)');ylabel('Phase of H(j\Omega)|');title('ÍÖԲ˫ÏßÐԱ任');
实验结果：

椭圆型：
fp = 1000;
fs = 1500;
f = 10000;
Wp = 2 * pi * fp;
Ws = 2 * pi * fs;
wp = 2 * f * tan(2 * pi * fp / (2 * f));
ws = 2 * f * tan(2 * pi * fs / (2 * f));
Rp = 1;
As = 15;
e = (10 ^ (Rp / 10) - 1) ^ 0.5;
K = Wp / Ws;
k1 = e / ((As * As) - 1) ^ 0.5;
k = wp / ws;
N1 = ceil(ellipke(K) * ellipke((1 - k1) ^ 0.5) / (ellipke(k1) * ellipke((1 - K) ^ 0.5)));
N2 = ceil(ellipke(k) * ellipke((1 - k1) ^ 0.5) / (ellipke(k1) * ellipke((1 - k) ^ 0.5)));
[b1, a1] = ellip(N1, Rp, As, Ws, 's');
[H, w] = freqs(b1, a1);
subplot(221);
plot(w/pi, abs(H));
grid on;
xlabel('\Omega(\pi)');ylabel('|H(j\Omega)|');title('ÍÖÔ²Âö³åÂö³åÏìÓ¦²»±ä');
subplot(223);
plot(w/pi, angle(H));
grid on;
xlabel('\Omega(\pi)');ylabel('Phase of H(j\Omega)|');title('ÍÖÔ²Âö³åÂö³åÏìÓ¦²»±ä');
[b2, a2] = ellip(N2, Rp, As, ws, 's');
[H, w] = freqs(b2, a2);
subplot(222);
plot(w/pi, abs(H));
grid on;
xlabel('\Omega(\pi)');ylabel('|H(j\Omega)|');title('ÍÖԲ˫ÏßÐԱ任');
subplot(224);
plot(w/pi, angle(H));
grid on;
xlabel('\Omega(\pi)');ylabel('Phase of H(j\Omega)|');title('ÍÖԲ˫ÏßÐԱ任');
实验结果：

分析：脉冲响应不变法的优点是频率坐标的变换是线性的，因此如果模拟滤波器的频响是限带于折叠频率以内的话，通过变换后的数字滤波器的频响可以不失真地反映原响应与频率之间的关系。但是，其最大的缺点是有频谱的周期延拓效应，会产生混叠失真。故，脉冲响应不变法只能用于限带的频响特性，如衰减特性较好的低通，或带通，而且高频衰减越大，频响的混叠效应就越小。双线性变换法的优点是消除了脉冲响应不变法的固有的混叠失真，但缺点是频率变换之间的非常严重的非线性，会使变换后时域上的图像产生非常严重的失真。

五、心得与体会
通过实验，我学习了利用脉冲响应不变法设计IIR数字滤波器的基本原理为：从时域响应出发，使数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的单位冲激响应ha(t)，h(n)等于ha(t)的取样值。所以设计中，应先根据数字指标转换来的模拟指标，设计模拟滤波器，再利用impinvar函数实现脉冲响应不变法模拟滤波器到数字滤波器的变换。
脉冲响应不变法和双线性法设计IIR数字滤波器是非常实用的技能，也会时我对DSP这门课程有更深刻的理解和认识。
本次实验最大的收获是了解了多种模拟滤波器到数字滤波器的转换方法。通过实验，确实巩固了课本知识。完成实验报告的过程中，查阅了相关资料，不仅开拓了知识面，更加深了对本部分知识的理解。

实验4 频率取样法设计FIR数字滤波器
一、实验目的
掌握频率取样法设计FIR数字滤波器的原理及具体方法

二、实验原理
1、基本原理
频率取样法从频域出发，把理想的滤波器等间隔采样得到,将作为实际设计滤波器的。
（8-1）
得到以后可以由来确定唯一确定滤波器的单位脉冲响应,可以由求得：
（8-2、3）
其中为内插函数
（8-4）
由求得的频率响应将逼近
如果我们设计的是线性相位FIR滤波器，则的幅度和相位满足线性相位滤波器的约束条件。
我们将表示为如下形式：
（8-5）
当为实数，则
由此得到
（8-6）
即为中心偶对称。在利用线性相位条件可知，对于1型和2型线性相位滤波器
（8-7）
对于3型和4型线性相位滤波器
（8-8）
其中，x表示取小于该数的最大的整数。
设计步骤
（1）由给定的理想滤波器给出和。
（2）由求得
（3) 根据求得或

三、实验内容
1. 设计一个数字低通FIR滤波器，其技术指标如下：
ω_p=0.2π，R_p=0.25dB
ω_st=0.3π，A_s=50dB
分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。结合实验结果分别讨论上述方法设计的数字滤波器是否符合指标。
2. 设计一个数字带通FIR滤波器，其技术指标如下：
3. 采用频率取样设计法设计FIR数字低通滤波器，满足以下指标
ω_p=0.2π，R_p=0.25dB
ω_st=0.3π，A_s=50dB
（1）取N=20，过渡带没有样本。
（2）取N=40，过渡带有一个样本，T=0.39。
（3）取N=60，过渡带有两个样本，T1=0.5925, T2=0.1099。
（4）分别采用上述方法设计的数字滤波器是否都能满足给定的指标要求。
4. 采用频率取样技术设计下面的高通滤波器
ω_st=0.6π，A_s=50dB
ω_p=0.8π，R_p=1dB
对于高通滤波器，N必须为奇数（或1型滤波器）。选择N=33，过渡带有两个样本，过渡带样本最优值为T1=0.1095，T2=0.598。

四、实验代码及实验结果
第一题
矩形窗：
实验代码：
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));
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，所以，不符合指标。
汉宁窗：
实验代码：
wp = 0.2 * pi; wst = 0.3 * pi;
tr_width = wst - wp;
N = ceil