《数字信号处理》实验报告

2017年5月

实验一 离散时间信号和系统响应

一. 实验目的

1. 掌握连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解

2. 掌握求系统响应的方法

3. 掌握时域离散系统的时域特性

4. 利用卷积方法观察分析系统的时域特性

二. 实验原理与方法

采样是连续信号数字化处理的第一个环节。对采样过程的研究不仅可以了解采样前后信号时域和频域特性的变化以及信号信息不丢失的条件，而且可以加深对离散傅里叶变换、Z变换和序列傅里叶变换之间关系式的理解。

对连续信号以T为采样间隔进行时域等间隔理想采样，形成采样信号：

式中为周期冲激脉冲，

的傅里叶变换为：

上式表明将连续信号xa(t)采样后其频谱将变为周期的，周期为Ωs=2π/T。也即采样信号的频谱是原连续信号xa(t)的频谱Xa(jΩ)在频率轴上以Ωs为周期，周期延拓而成的。

在时域中，描述系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。一个时域离散线性时不变系统的输出与输入间的关系为：

这里，为系统的输出序列，为输入序列。可以是无限长，也可以是有限长。

三. 实验内容

1. 时域采样定理的验证

给定模拟信号：

，

式中。

对其进行采样，可得到采样序列

。

选择三种采样频率Fs=1kHz, 300Hz, 200Hz, 观测时间选Tp=64ms，分别得到三个采样序列。编写程序计算这三个序列的幅度特性，并绘图显示其幅频特性曲线。观察分析频谱混叠现象。

2. 给定一个低通滤波器的差分方程为：

输入信号

a. 分别求出和的系统响应，并画出其波形

b. 求出系统的单位脉冲响应，画出其波形

3. 给定系统的单位脉冲响应为

用线性卷积法求分别对系统和的输出响应，并画出波形

A=444.128;

a=50*sqrt(2)*pi;

w=50*sqrt(2)*pi;

n=0:63;

fs = 1000;

c=A*exp((-a)*n/fs).*sin(w*n/fs);

subplot(3,2,1);

stem(n,c,'.');

xlabel('n');

ylabel('xa(n)');

title('xa(n)的时域序列');

N=64;

k=0:length(c)-1;

w=k*pi/100;

X=fft(c,N);

subplot(3,2,2);

plot(w/pi,abs(X));

xlabel('w/pi');

ylabel('|X(jw)|');

title('xa(n)的傅氏变换|X(jw)|');

A=444.128;

a=50*sqrt(2)*pi;

w=50*sqrt(2)*pi;

n=0:63;

fs = 300;

c=A*exp((-a)*n/fs).*sin(w*n/fs);

subplot(3,2,3);

stem(n,c,'.');

xlabel('n');

ylabel('xa(n)');

title('xa(n)的时域序列');

N=64;

k=0:length(c)-1;

w=k*pi/100;

X=fft(c,N);

subplot(3,2,4);

plot(w/pi,abs(X));

xlabel('w/pi');

ylabel('|X(jw)|');

title('xa(n)的傅氏变换|X(jw)|');

A=444.128;

a=50*sqrt(2)*pi;

w=50*sqrt(2)*pi;

n=0:63;

fs = 200;

c=A*exp((-a)*n/fs).*sin(w*n/fs);

subplot(3,2,5);

stem(n,c,'.');

xlabel('n');

ylabel('xa(n)');

title('xa(n)的时域序列');

N=64;

k=0:length(c)-1;

w=k*pi/100;

X=fft(c,N);

subplot(3,2,6);

plot(w/pi,abs(X));

xlabel('w/pi');

ylabel('|X(jw)|');

title('xa(n)的傅氏变换|X(jw)|');

内容一：调用filter解差分方程，由系统对u(n)的响应判断稳定性

A=[1,-0.9];

B=[0.05,0.05];

x1n=[1 1 1 1 1 1 1 1 zeros(1,122) ];

x2n=ones(1,122);

hn=impz(B,A,58);

subplot(2,2,1);

y='hn';

stem(hn,'r','.');

xlabel('n');

ylabel('h(n)');

title('(a)系统单位脉冲响应h(n)');

y1n=filter(B,A,x1n);

subplot(2,2,2);

y='y1n';

stem(y1n,'g','.');

xlabel('n');

ylabel('h(n)');

title('(b)系统对R8(n)的响应y1(n)');

y2n=filter(B,A,x2n);

subplot(2,2,4);

y='y1n';

stem(y2n,'b','.');

xlabel('n');

ylabel('y2(n)');

title('(c)系统对u(n)的响应y2(n)');

内容一:系统的单位冲响应的波形如下图(a)所示，系统 和 的响应序列的波形如下图(b)和图(c)

%内容2:调用conv函数计算卷积

x1n=[1 1 1 1 1 1 1 1];

h1n=[ones(1,10) zeros(1,22)];

h2n=[1 2.5 2.5 1 zeros(1,22)];

y21n=conv(h1n,x1n);

y22n=conv(h2n,x1n);

figure(2);

subplot(2,2,1);

y='h1(n)';

stem(h1n,'b','.');

axis([0 22 0 1.2]);

title('(d)系统单位脉冲响应h1(n)');

subplot(2,2,2);

y='y21n';

stem(y21n,'r','.');

axis([0 52 0 8.2]);

title('(e)h1(n)与R8(n)的卷积y21(n)');

subplot(2,2,3);

y='h2(n)';

stem(h2n,'b','.');

axis([0 22 0 3.2]);

title('(f)系统对单位脉冲响应h2(n)');

subplot(2,2,4);

y='y22n';

stem(y22n,'r','.');

grid on;

title('(g)h2(n)与R8(n)的卷积y22(n)');

四. 实验思考

1. 在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？它们所对应的模拟频率是否相同？为什么？；

对于生活中的连续信号（图像，语音等），采样是十分重要的，采样频率必须满足采样定理（f>2fs）才能安全重构原来的信号，当采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是不同的，它们所对应的模拟频率也是不同的。因为采样的频率不同，其傅里叶变换所能包含的范围也不同。而我们是在抽样的过程中考虑采样频率，在FFT中其实可以不考虑，归一化都可以，所以说提高采样率我们的重点是能不能够重构原来的信号，如果原连续信号的最高频率分量或者带宽很宽的话，那么按照采样定理就要提高采样率了。

2. 如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应？如何求？

直接进行线性卷积是不可行的，因为输入信号为无限长序列。但在实验中函数fftt的调用格式为 y=fftfilt(b,x) 该格式是利用基于FFT的重叠相加法对数据进行滤波,这种频域滤波技术只对FIR滤波器有效.该函数是通过向量b描述的滤波器对x数据进行滤波.

x是等待滤波的信号；

b是FIR滤波器的H(z）的分子多项式系数

3. 如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号会有何变化？用前面第二个实验结果进行分析说明

如果信号经过低通滤波器,把信号的高频分量滤掉,时域信号将变得更加平滑。并且由实验内容（2）可见,经过系统低通滤波使输入信号、和的阶跃变化变得缓慢上升与下降.

实验总结

(1) 时域求系统响应的方法有两种，第一种是通过解差分方程求得系统的响应输出，注意合理的选择系统的初始条件，第二种是已知系统的单位冲激响应，通过求系统输入与单位脉冲响应的线性卷积求得系统的输出。

(2) 实验中要检查系统的稳定性，其方法是在输入端加上单位阶跃序列，观察系统的输出波形，如果波形的幅值稳定在一个常数值附近，则系统稳定，否则系统不稳定。上面第三个实验是稳定的。

实验二 用FFT对信号作频谱分析

一. 实验目的

1. 进一步加深DFT算法原理和基本性质的理解

2. 掌握用FFT对连续信号和时域离散信号进行频谱分析的方法

3. 了解用FFT进行频谱分析时可能出现的分析误差及其原因，以便正确应用FFT。

二. 实验原理

（补充离散傅里叶变换的性质）

用FFT对信号作频谱分析是学习数字信号处理的重要内容，经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率F和分析误差。频谱分辨率直接和FFT的变换区间N有关，FFT能够实现的频率分辨率是2/N，因此要求2/NF。可以根据此式选择FFT的变换区间N。误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时，离散谱的包络才能逼近于连续谱，因此N要适当选择大一些。

周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察世间长一些。

对模拟信号进行谱分析时，首先要按照采样定理将其变为时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。

三. 实验内容

1. 对以下给出的各序列进行谱分析：

选择FFT的变换区间N为8和16两种情况进行频谱分析。

2. 对以下各周期序列进行频谱分析

选择FFT的变换区间N为8和16两种情况分别对以上序列进行频谱分析。

3. 对模拟周期信号进行频谱分析

选择采样频率Fs=64Hz，对变换区间N=16，32，64三种情况进行谱分析

1. 对以下序列进行谱分析:

选择FFT的变换区间N为8和16 两种情况进行频谱分析。分别打印其幅频特性曲线, 并进行对比、分析和讨论。

代码如下：

x1n=[ones(1,4)];%产生R4（n）序列向量

X1k8=fft(x1n,8);

X1k16=fft(x1n,16);

N=8;

f=2/N*(0:N-1);

figure(1);

subplot(1,2,1);

stem(f,abs(X1k8),'.');

title('(1a)8点DFT[x_1(n)]');

xlabel('w/pai');

ylabel('幅度');

N=16;

f=2/N*(0:N-1);

figure(1);

subplot(1,2,2);

stem(f,abs(X1k16),'.');

title('(1a)16点DFT[x_1(n)]');

xlabel('w/pai');

ylabel('幅度');

M=8;

xa=1:(M/2);

xb=(M/2):-1:1;

x2n=[xa,xb];

x3n=[xb,xa];

X2k8=fft(x2n,8);

X2k16=fft(x2n,16);

X3k8=fft(x3n,8);

X3k16=fft(x3n,16);

figure(2);

N=8;

f=2/N*(0:N-1);

subplot(2,2,1);

stem(f,abs(X2k8),'.');

title('(2a)8点DFT[x_1(n)]');

xlabel('w/pai');

ylabel('幅度');

subplot(2,2,3);

stem(f,abs(X3k8),'.');

title('(3a)8点DFT[x_1(n)]');

xlabel('w/pai');

ylabel('幅度');

N=16;

f=2/N*(0:N-1);

subplot(2,2,2);

stem(f,abs(X2k16),'.');

title('(2a)16点DFT[x_1(n)]');

xlabel('w/pai');

ylabel('幅度');

subplot(2,2,4);

stem(f,abs(X3k16),'.');

title('(3a)16点DFT[x_1(n)]');

xlabel('w/pai');

ylabel('幅度');

N=8;

n=0:N-1;

x4n=cos(pi*n/4);

x5n=cos(pi*n/4)+cos(pi*n/8);

X4k8=fft(x4n,8);

X4k16=fft(x4n,16);

X5k8=fft(x5n,8);

X5k16=fft(x5n,16);

figure(3);

N=8;

f=2/N*(0:N-1);

subplot(2,2,1);

stem(f,abs(X4k8),'.');

title('(4a)8点DFT[x_4(n)]');

xlabel('w/pai');

ylabel('幅度');

subplot(2,2,3);

stem(f,abs(X5k8),'.');

title('(5a)16点DFT[x_5(n)]');

xlabel('w/pai');

ylabel('幅度');

N=16;

f=2/N*(0:N-1);

subplot(2,2,2);

stem(f,abs(X4k16),'.');

title('(4a)16点DFT[x_4(n)]');

xlabel('w/pai');

ylabel('幅度');

subplot(2,2,4);

stem(f,abs(X5k16),'.');

title('(5a)16点DFT[x_5(n)]');

xlabel('w/pai');

ylabel('幅度');

四. 实验思考

1. 在N=8时，和的幅频特性会相同吗？为什么？N=16时呢？

答：当N=8时，幅频特性相同，由时域循环移位定理，x2(n)和x3(n)满足时域循环移位定理。

当N=16时，它们的幅频特性不会相同，因为，x2(n)和x3(n)后都要补零，x2(n)进行循环移位后的序列跟x3(n)不相同。

2. 对于周期序列，如果周期不知道，如何用FFT进行谱分析？

答：周期信号的周期预先不知道时,可先截取M点进行DFT,再将截取长度扩大1倍截取,比较结果,如果二者的差别满足分析误差要求,则可以近似表示该信号的频谱,如果不满足误差要求就继续将截取长度加倍,重复比较,直到结果满足要求

3. 如何选择FFT的变换区间？（包括非周期信号和周期信号）

答：一、对于非周期信号：有频谱分辨率F，而频谱分辨率直接和FFT的变换区间有关，因为FFT能够实现的频率分辨率是2π/N...因此有最小的N>2π/F。就可以根据此式选择FFT的变换区间。

二、对于周期信号，周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。

答：周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。

对模拟信号进行频谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的普分析进行。

实验三 IIR数字滤波器的设计

一. 实验目的

1. 熟悉用双线性变换法设计IIR数字滤波器的原理和方法

2. 掌握数字滤波器的Matlab实现方法

3. 通过观察对实际心电图信号的滤波作用，获得数字滤波的感性认识

二. 实验原理

设计IIR数字滤波器一般采用间接设计法——脉冲响应不变法和双线性变换法，应用最广泛的是双线性变换法。

脉冲响应不变法的基本思想是：使数字滤波器的单位脉冲响应h(n)近似于模拟滤波器的单位脉冲响应ha(t)，即使

其S平面和Z平面的映射关系为：

双线性变换法的基本思想是：使描述数字滤波器的差分方程近似描述模拟滤波器的微分方程

S平面和Z平面的映射关系为：

双线性变换法中的频率变换是一种非线性变换，这种非线性引起的幅频特性畸变可通过预变形矫正法而得到校正。

设计IIR 数字滤波器的一般步骤:

（1）确定所需类型数字滤波器的技术指标：通带截止频率ωp、通带衰减αp、阻带截止频率ωs、阻带衰减αs。

（2）将所需类型数字滤波器的技术指标转换成相应类型模拟滤波器的技术指标。

（3）设计该类型模拟滤波器

（4）通过复频率变换将模拟滤波器转换成所需类型的数字滤波器。

三. 实验内容

1. 分别用脉冲响应不变法和双线性变换法设计一个巴特沃斯低通IIR数字滤波器，设计指标参数为：在通带内频率低于0.2时，最大衰减小于1dB，在阻带内[0.3，]频率区间上，最小衰减大于15dB。观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。

Fp=0.1;

Fs=0.15;

Ft=1;

Nn=256;

As=15;

Ap=1

wp=2*pi*Fp;

ws=2*pi*Fs;

[n,wn]=buttord(wp,ws,Ap,As,'s'); %求低通滤波器的阶数和截止频率；n为滤波器最小阶数，wn为其截止频率，

%wp，ws分别为通带频率和截止频率，Ap为通带最大衰减，As阻带最小衰减。

[b,a]=butter(n,wn,'s'); %求S域的频率响应的参数

[bz,az]=impinvar(b,a,Ft); %脉冲响应不变法实现S域到Z域的变换

[h,f]=freqz(bz,az,Nn,Fs); %根据参数求出频率响应

figure(1)

subplot(2,1,1);plot(f,20*log10(abs(h)));

xlabel('频率/Hz');ylabel('振幅/dB');

grid on;

subplot(2,1,2);plot(f,180/pi*unwrap(angle(h)));

xlabel('频率/Hz');ylabel('相位/^o');

grid on;

title('脉冲响应不变法');

[bz,az]=bilinear(b,a,Ft); %双线性变换实现S域到Z域的变换

[h,f]=freqz(bz,az,Nn,Fs); %根据参数求出频率响应

figure(2)

subplot(2,1,1);plot(f,20*log10(abs(h)));

xlabel('频率/Hz');ylabel('振幅/dB');

grid on;

subplot(2,1,2);plot(f,180/pi*unwrap(angle(h)));

xlabel('频率/Hz');ylabel('相位/^o');

grid on;

title('双线性变换法');

3.心电图

xn=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-19,-38,-60,84,-90,-66,-32,-4,-2,-4,8,12,12, 10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];

n=0:55;

Fp=0.1;

Fs=0.15;

Ft=1;

Nn=256;

As=15;

Ap=1

wp=2*pi*Fp;

ws=2*pi*Fs;

[k,wk]=buttord(wp,ws,Ap,As,'s'); %求低通滤波器的阶数和截止频率；k滤波器最小阶数，wk截止频率，

%wp，ws分别为通带频率和截止频率，Ap为通带最大衰减，As阻带最小衰减。

[b,a]=butter(k,wk,'s'); %求S域的频率响应的参数

[bz,az]=bilinear(b,a,Ft); %双线性变换实现S域到Z域的变换

[h,f]=freqz(bz,az,Nn,Fs); %根据参数求出频率响应

figure(1)

stem(n,xn);

title('原始心电信号');

y1=filter(bz,az,xn);

figure(2)

plot(n,y1);

title('滤波后的心电信号')

四. 实验思考

1. 用双线性变换法设计数字滤波器过程中，变换公式

中T的取值，对设计结果有无影响？为什么？

答：无影响，因为数字滤波器的传输函数H（ejw）以2为周期，最高频率在处，因此， <，按照线性关系,那么一定满足,因此T可以任选。在数字滤波器的设计过程中，设计问题起始于数字滤波器上的传输函数的技术要求，当这些技术要求通过变换公式得到模拟滤波器的传输函数HC (s)的技术要求并设计出模拟滤波，然后模拟滤波器的传输函数HC(s)再通过公式变换得到数字滤波器的传输函数，这样取值T对设计的影响就被抵消，因此T的取值对设计结果没有影响。

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

一. 实验目的

1. 掌握用窗函数法设计FIR数字滤波器的原理和方法

2. 熟悉线性相位FIR数字滤波器特性

3. 了解各种窗函数对滤波特性的影响

二. 实验原理

窗函数法设计 FIR 滤波器的步骤为：

（1）构建希望逼近的理想频率响应函数及技术指标

（2）求滤波器的单位脉冲响应

如果复杂，可对从采样M个点，采样值为，则：

（3）根据对过渡带及阻带衰减的要求，选择窗函数的形式，并估计窗口宽度N，设要求的过渡带宽为，则

（4）计算滤波器的单位脉冲响应

（5）求H(ejω)，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。

窗函数傅里叶变换W(ejω)的主瓣决定了H(ejω)过渡带宽，W(ejω)的旁瓣大小和多少决定了H(ejω)在通带和阻带范围内波动幅度，常用的几种窗函数有：

矩形窗；Hanning窗；Hamming窗；Blackmen窗；Kaiser窗

三. 实验内容

1.用升余弦窗设计一线性相位低通FIR 数字滤波器，截止频率。窗口长度N=15，33。要求在两种窗口长度情况下，分别求出h(n)，打印出相应的幅频特性和相频特性曲线，观察3dB带宽和20dB 带宽，总结窗口长度N对滤波特性的影响

2. n=33, , 用四种窗函数设计线性相位低通滤波器，绘制相应的幅频特性曲线，观察3dB带宽和20dB 带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响

实验代码：

主程序：

s = -1;

while(s<0)

clc;

N=input('请输入窗函数长度 N=');

s=1;

end

close all;

i=0;wc=pi/4;

while(s)

n=0:N-1;

hd=ideal(wc,N);

k=input('请选择窗口类型:\n[1](boxcar)\n[2](hamming)\n[3](hanning)\n[4](blackman)\n请选择：','s');

k=str2num(k);

if(k==1)

B=boxcar(N);

string=['Boxcar','N=',num2str(N)];

else

if(k==2)

B=hamming(N);

string=['Hamming','N=',num2str(N)];

else

if(k==3)

B=hanning(N);

string=['Hanning','N=',num2str(N)];

else

if(k==4)

B=blackman(N);

string=['Blackman','N=',num2str(N)];

end

end

end

end

h=hd.*(B)';

[H,f]=freqz(h,[1],1024,'whole');

db=20*log10(abs(H)/max(abs(H)));

pha=angle(H);

i=i+1;

figure(i)

subplot(2,2,1);

n=0:N-1;

stem(n,h,'.');

xlabel('n');

ylabel('h(n)');

title('实际低通滤波器的h(n)');

text((0.3*N),0.27,string);

hold off;

subplot(2,2,2);

plot(f/pi,db);

axis([0 1 -100 0]);

xlabel('w/pi');

ylabel('20log|H(jw)|');

title('衰减特性(dB)');

grid;

subplot(2,2,3);

plot(f,pha);

hold on;

n=0:7;

x=zeros(8);

plot(n,x);

title('相频特性');

xlabel('频率w（rad)');

ylabel('相位（rad）');

axis([0 3.15 -4 4]);

subplot(2,2,4);

plot(f,abs(H));

title('幅频特性');

xlabel('频率W(rad)');

ylabel('幅值|H(jw)|');

axis([0 3.15 0 1.5]);

text(0.9,1.3,string);

s=input('继续实验请选择:\n[1](继续实验)\n[0](退出)\n请选择:');

if(s==1)

N=input('请输入窗函数的长度 N=');

end

End

子程序：

function hd=ideal(w,N);

alpha=(N-1)/2;

n=0:N-1;

m=n-alpha+eps;

hd=sin(w*m)./(pi*m);

（1）N=15时的矩形窗：

N=33时的矩形窗：

(2) N=15时的汉明窗：

N=33时的汉明窗：

(3) N=15时的汉宁窗：

N=33时的汉宁窗：

(4) N=15时的布莱克曼窗：

N=33时的布莱克曼窗：

（1）矩形窗的频率响应主瓣宽度为4π/N，第一副瓣比主瓣低13dB。

（2）哈明窗是改进的升余弦窗，能量更加集中在主瓣中，主瓣的能量约占99.96%，第一旁瓣的峰值比主瓣小40dB,但主瓣宽度和汉宁窗相同，为8π/N。

（3）汉宁窗的幅度函数由三部分相加，使能量更集中在主瓣中，但代价是主瓣宽度加宽到8π/N。

（4）布莱克曼窗的幅度函数由五部分组成，它们都是移位不同，且幅度也不同的函数，使旁瓣再进一步抵消，阻带衰减进一步增加，过渡带是矩形窗过渡带的3倍。

（5）调整窗口长度N可以有效低控制过渡带的宽度，减少带内波动以及加大阻带的衰减只能从窗函数的形状上找解决方法，如果能找到的窗函数形状，使其谱函数的主瓣包含更多的能量，相应旁瓣幅度就减小了，旁瓣的减小可使通带、阻带波动减小，从而加大阻带衰减，但这样总是以加宽过渡带为代价的。

3. 调用信号产生函数xtg产生具有加性噪声的信号xt，显示xt波形及其频谱。设计一低通滤波器，从高频噪声中提取xt中的单频调幅信号，要求信号幅频失真小于0.1dB,噪声频谱衰减不小于60dB。

（1）观察xt的频谱，确定滤波器指标参数。

（2）根据滤波器指标选择合适的窗函数，计算窗函数的长度N，设计一个FIR低通滤波器，绘图显示滤波器的频响特性曲线

（3）用设计的FIR低通滤波器对xt进行滤波，绘图显示滤波器输出信号的时域和频域波形图

实验程序框图如图所示。

1、信号产生函数xtg程序清单：

%xt=xtg(N) 产生一个长度为N,有加性高频噪声的单频调幅信号xt,采样频率Fs=1000Hz

%载波频率fc=Fs/10=100Hz,调制正弦波频率f0=fc/10=10Hz.

function xt=xtg

N=1000;Fs=1000;T=1/Fs;Tp=N*T;

t=0:T:(N-1)*T;

fc=Fs/10;f0=fc/10; %载波频率fc=Fs/10，单频调制信号频率为f0=Fc/10;

mt=cos(2*pi*f0*t); %产生单频正弦波调制信号mt，频率为f0

ct=cos(2*pi*fc*t); %产生载波正弦波信号ct，频率为fc

xt=mt.*ct; %相乘产生单频调制信号xt

nt=2*rand(1,N)-1; %产生随机噪声nt

%=======设计高通滤波器hn,用于滤除噪声nt中的低频成分,生成高通噪声=======

fp=120; fs=150;Rp=0.2;As=60; % 滤波器指标

fb=[fp,fs];m=[0,1]; % 计算remezord函数所需参数f,m,dev

dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];

[n,fo,mo,W]=remezord(fb,m,dev,Fs); % 确定remez函数所需参数

hn=remez(n,fo,mo,W); % 调用remez函数进行设计,用于滤除噪声nt中的低频成分

yt=filter(hn,1,10*nt); %滤除随机噪声中低频成分，生成高通噪声yt

%================================================================

xt=xt+yt; %噪声加信号

fst=fft(xt,N);k=0:N-1;f=k/Tp;

subplot(3,1,1);plot(t,xt);grid;xlabel('t/s');ylabel('x(t)');

axis([0,Tp/5,min(xt),max(xt)]);title('(1) 信号加噪声波形')

subplot(3,1,2);plot(f,abs(fst)/max(abs(fst)));grid;title('(2) 信号加噪声的频谱')

axis([0,Fs/2,0,1.2]);xlabel('f/Hz');ylabel('幅度')

2、主程序/实验程序清单：

clear all;clear all;

%==调用xtg产生信号xt, xt长度N=1000,并显示xt及其频谱,=========

N=1000;xt=xtg;

fp=120; fs=150;Rp=0.2;As=60;Fs=1000; % 输入给定指标

% (1) 用窗函数法设计滤波器

wc=(fp+fs)/Fs; %理想低通滤波器截止频率(关于pi归一化）

B=2*pi*(fs-fp)/Fs; %过渡带宽度指标

Nb=ceil(11*pi/B); %blackman窗的长度N

hn=fir1(Nb-1,wc,blackman(Nb));

Hw=abs(fft(hn,1024)); % 求设计的滤波器频率特性

ywt=fftfilt(hn,xt,N); %调用函数fftfilt对xt滤波

%以下为用窗函数法设计法的绘图部分（滤波器损耗函数，滤波器输出信号波形）

f=[0:1023]*Fs/1024;

figure(2)

subplot(2,1,1)

plot(f,20*log10(Hw/max(Hw)));grid;title('(3) 低通滤波器幅频特性')

axis([0,Fs/2,-120,20]);

xlabel('f/HZ');ylabel('幅度')

t=[0:N-1]/Fs;Tp=N/Fs;

subplot(2,1,2)

plot(t,ywt);grid;

axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_w(t)');

title('(4) 滤波噪声后的信号波形')

% (2) 用等波纹最佳逼近法设计滤波器

fb=[fp,fs];m=[1,0]; % 确定remezord函数所需参数f,m,dev

dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)];

[Ne,fo,mo,W]=remezord(fb,m,dev,Fs); % 确定remez函数所需参数

hn=remez(Ne,fo,mo,W); % 调用remez函数进行设计

Hw=abs(fft(hn,1024)); % 求设计的滤波器频率特性

yet=fftfilt(hn,xt,N); % 调用函数fftfilt对xt滤波

%以下为用等纹波设计法的绘图部分（滤波器损耗函数，滤波器输出信号波形）

f=[0:1023]*Fs/1024;

figure(3)

subplot(2,1,1)

plot(f,20*log10(Hw/max(Hw)));grid;title('(5) 低通滤波器幅频特性')

axis([0,Fs/2,-80,10]);

xlabel('f/HZ');ylabel('幅度')

%t=[0:N-1]/Fs;Tp=N/Fs;

subplot(2,1,2)

plot(t,yet);grid;

axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_e(t)');

title('(6) 滤波噪声后的信号波形')

实验图形如下：

四. 实验思考

1. 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器？

①根据所需设计的数字滤波器类型，确定线性相位数字滤波器类型； 

②选择合适的窗函数； 

③确定理想低通数字滤波器的频率响应函数； 

④计算理想低通数字滤波器的单位脉冲响应； 

⑤加窗得到设计结果； 

2. 如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为和，试求理想带通的单位脉冲响应

答：希望逼近的理想带通滤波器的截止频率分别为：

。

用窗函数法设计FIR滤波器的主要特点：

设＝FT[]为希望逼近的频响特性函数,H() =FT[h(n)]为用窗函数法设计的实际滤波器的频响函数。通常取H()相应的理想频响特性作为。

实验总结

（1）希望逼近的理想滤波器频响函数的表达式。因为数字滤波器一般要求设计成线性相位特性，所以必须满足上述线性相位FIR滤波器的频域特点。

（2）熟悉各种常用窗函数的技术指标和加窗后对滤波特性的影响，这样才能根据设计指标正确选择窗函数类型及其长度N。

（3）检验设计结果：窗函数法的设计结果单位脉冲响应。而检验一般在频域进行。所以要计算检验3dB截止频率和阻带最小衰减，其计算量相当大，必须用计算机进行。

（4）熟悉窗函数设计法的特点：设计过程简单.方便实用。但边界频率不易精确控制。所以设计完以后，必须检验结果。