c语言中butte函数是什么意思

C语言没有这个函数

MATLAB中有这个函数

说明

巴特沃斯滤波器设计

语法

[z,p,k]=butter(n,Wn)

[z,p,k] = butter(n,Wn,'ftype')

[b,a]=butter(n,Wn)

[b,a]=butter(n,Wn,'ftype')

[A,B,C,D]=butter(n,Wn)

[A,B,C,D] = butter(n,Wn,'ftype')

[z,p,k]=butter(n,Wn,'s')

[z,p,k] = butter(n,Wn,'ftype','s')

[b,a]=butter(n,Wn,'s')

[b,a]=butter(n,Wn,'ftype','s')

[A,B,C,D]=butter(n,Wn,'s')

[A,B,C,D] = butter(n,Wn,'ftype','s')

描述

butter 用来设计低通、带通、高通、和带阻数字和模拟的巴特沃斯滤波器。巴特沃斯滤波器的特征是通带内幅度响应最大平坦，且整体上是单调的。

巴特沃斯滤波器牺牲了在通带和阻带内的单调衰减陡度。除非需要巴特沃斯滤波器的平滑性，椭圆或切比雪夫滤波器可以用更小的滤波器阶数获得更陡峭的衰减特性。

数字域

[z,p,k] = butter(n,Wn) 设计一个阶数为n，归一化截止频率为Wn的低通数字巴特沃斯滤波器。此函数用n列的向量z和p返回零点和极点，以及用标量k返回增益。

[z,p,k] = butter(n,Wn,'ftype') 设计一个高通、低通或带阻滤波器，字符串'ftype'取值是：

'high' 用于设计归一化截止频率为Wn的高通数字滤波器

'low' 用于设计归一化截止频率为Wn的低通数字滤波器

'stop' 用于设计阶数为2n的带阻数字滤波器，Wn应该是有两个元素的向量Wn=[w1 w2]。阻带是w1 < ω <w2

截止频率 是幅度响应为处的的频率。对巴特沃斯滤波器，归一化截止频率Wn必须是介于0和1之间的数，这里的1对应于尼奎斯特频率，即每秒π弧度(π rad/s)。

如果Wn是含有两个元素的向量，Wn=[w1 w2]，butter 返回阶数为 2n的数字带通滤波器，通频带为w1 < ω <w2

通过返回不同数量的输出参数，butter 直接地得到其它的滤波器实现。要获得传输函数形式，使用如下所示的两个输出参数。

注意 参考下面的限制 了解关于影响形成传输函数的数值问题。

[b,a] = butter(n,Wn) 设计一个阶为n，归一化截止频率为Wn的数字低通巴特沃斯滤波器。它返回滤波器系数在长度为n+1的行向量b和a中，这两个向量包含z的降幂系数。

[b,a] = butter(n,Wn,'ftype') 设计一个高通、低通或带阻滤波器，字符串'ftype' 是上面描述的'high'、'low'、或 'stop'。

要获得状态空间形式，使用下面所示的4个输出参数：

[A,B,C,D] = butter(n,Wn) 或

[A,B,C,D] = butter(n,Wn,'ftype') 其中 A、 B,、C,、和D 是

并且u是输入, x是状态向量, y 是输出。

模拟域

[z,p,k] = butter(n,Wn,'s') 设计一个阶n，截止角频率为Wn rad/s的模拟低通巴特沃斯滤波器。它返回零点和极点在长n或2n的列向量z和p中，标量k返回增益。butter的截止角频率Wn必须大于0 rad/s。

如果Wn是有两个元素w1<w2的向量， butter(n,Wn,'s') 返回阶 2n 带通模拟滤波器，其通带是w1 < ω < w2。

[z,p,k] = butter(n,Wn,'ftype','s') 通过使用上面描述的ftype 值可以设计一个高通、低通或带阻滤波器。

只要返回不同数量的输出参数，butter 可以直接地获得其它的模拟滤波器实现。要获得传输函数形式，使用如下所示的两个输出参数：

[b,a] = butter(n,Wn,'s') 设计一个阶n、截止角频率为Wn rad/s的模拟低通巴特沃斯滤波器。它返回滤波器的系数在长n+1的行向量b和a中，这两个向量包含下面这个传输函数中s的降幂系数：

[b,a] = butter(n,Wn,'ftype','s') 通过设置上面描述的ftype 值，可以设计一个高通、低通或带阻滤波器。

要获得状态空间形式，使用下面的四个参数：

[A,B,C,D] = butter(n,Wn,'s') 或

[A,B,C,D] = butter(n,Wn,'ftype','s') 其中A、 B、 C、和D 是

并且u 作为输入, x 是状态向量, y 是输出。

举例

高通滤波器

对于1000Hz的采样，设计一个9阶高通巴特沃斯滤波器，截止频率300Hz，相应的归一化值为06：

[z,p,k] = butter(9,300/500,'high');

[sos,g] = zp2sos(z,p,k); % 转换为二次分式表示形式

Hd = dfiltdf2tsos(sos,g); % 创建dfilt对象

h = fvtool(Hd); % 绘制幅度响应

set(h,'Analysis','freq') % 显示频率响应

11 实验目的

1．了解数字信号处理系统的一般构成；

2．掌握奈奎斯特抽样定理。

12 实验仪器

1．YBLD智能综合信号源测试仪 1台

2．双踪示波器 1台

3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台

4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台

13 实验原理

一个典型的DSP系统除了数字信号处理部分外，还包括A/D和D/A两部分。这是因为自然界的信号，如声音、图像等大多是模拟信号，因此需要将其数字化后进行数字信号处理，模拟信号的数字化即称为A/D转换。数字信号处理后的数据可能需还原为模拟信号，这就需要进行D/A转换。一个仅包括A/D和D/A两部分的简化数字信号处理系统功能如图1所示。

A/D转换包括三个紧密相关的过程，即抽样、量化和编码。A/D转换中需解决的以下几个重要问题：抽样后输出信号中还有没有原始信号的信息？如果有能不能把它取出来？抽样频率应该如何选择？

奈奎斯特抽样定理（即低通信号的均匀抽样定理）告诉我们，一个频带限制在0至fx以内的低通信号x(t)，如果以fs≥2fx的抽样速率进行均匀抽样，则x(t)可以由抽样后的信号xs(t)完全地确定，即xs(t)包含有x(t)的成分，可以通过适当的低通滤波器不失真地恢复出x(t)。最小抽样速率fs=2fx称为奈奎斯特速率。

低通

译码

编码

量化

抽样

输入信号 样点输出 滤波输出

A/D（模数转换） D/A（数模转换）

图1 低通采样定理演示

为方便实现，实验中更换了一种表现形式，即抽样频率固定（10KHz），通过改变输入模拟信号的频率来展示低通抽样定理。我们可以通过研究抽样频率和模拟信号最高频率分量的频率之间的关系，来验证低通抽样定理。

14 实验内容

1．软件仿真实验：编写并调试MATLAB程序，分析有关参数，记录有关波形。

2．硬件实验：输入不同频率的正弦信号，观察采样时钟波形、输入信号波形、样点输出波形和滤波输出波形。

15 MATLAB参考程序和仿真内容

%%

%f—余弦信号的频率

% M—基2 FFT幂次数 N=2^M为采样点数，这样取值是为了便于作基2的FFT分析

%2 采样频率Fs

%%

function samples(f,Fs,M)

N=2^M; % fft点数=取样总点数

Ts=1/Fs; % 取样时间间隔

T=NTs; % 取样总时间=取样总点数取样时间间隔

n=0:N-1;

t=nTs;

Xn=cos(2fpit);

subplot(2,1,1);

stem(t,Xn);

axis([0 T 11min(Xn) 11max(Xn)]);

xlabel('t -->');

ylabel('Xn');

Xk=abs(fft(Xn,N));

subplot(2,1,2);

stem(n,Xk);

axis([0 N 11min(Xk) 11max(Xk)]);

xlabel('frequency -->');

ylabel('!Xk!');

%%

假如有一个1Hz的余弦信号y=cos(2πt)，对其用4Hz的采样频率进行采样，共采样32点，只需执行samples(1,4,5)，即可得到仿真结果。

软件仿真实验内容如下表所示：

仿真参数

f

Fs

Wo(计算)

Xn（图形）

Xk（图形）

（1，4，5）

另外记录图形，并标图号

（1，8，5）

（2，8，6）

自 选

16 硬件实验步骤

本实验箱采样频率fs固定为10KHz，低通滤波器的截止频率约为45KHz。

1、用低频信号源产生正弦信号，正弦信号源频率f自定，并将其接至2TP2（模拟输入）端，将示波器通道一探头接至2TP6（采样时钟）端观察采样时钟波形，示波器通道二探头接至2TP2观察并记录输入信号波形。

2、将示波器通道二探头接至2TP3观察并记录样点输出波形。

3、将示波器通道二探头接至2TP4观察并记录滤波输出波形。

4、根据采样定理，分f=fs /8、f=fs/4、f=fs/2等3种情况更改正弦信号频率，重复步骤2至步骤3。

5、用低频信号源产生方波信号，重复步骤1至步骤4。

17 思考题

1、 讨论在仿真实验中所计算的数字域频率Wo和Xk的图形中非零谱线位置之间的对应关系。

2、 讨论在仿真实验中自选参数的意义。

3、将在2TP2端加方波信号后的恢复波形，与相同频率的正弦信号的恢复波形相比，能够得出哪些结论？

2 FFT频谱分析实验

21 实验目的

1．通过实验加深对快速傅立叶变换（FFT）基本原理的理解。

2．了解FFT点数与频谱分辨率的关系，以及两种加长序列FFT与原序列FFT的关系。

22 实验仪器

1．YBLD智能综合信号源测试仪 1台

2．双踪示波器 1台

3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台

4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台

23 实验原理

离散傅里叶变换（DFT）和卷积是信号处理中两个最基本也是最常用的运算，它们涉及到信号与系统的分析与综合这一广泛的信号处理领域。实际上卷积与DFT之间有着互通的联系：卷积可化为DFT来实现，其它的许多算法，如相关、滤波和谱估计等都可化为DFT来实现，DFT也可化为卷积来实现。

对N点序列x(n)，其DFT变换对定义为：

在DFT运算中包含大量的重复运算。FFT算法利用了蝶形因子WN的周期性和对称性，从而加快了运算的速度。FFT算法将长序列的DFT分解为短序列的DFT。N点的DFT先分解为2个N/2点的DFT，每个N/2点的DFT又分解为2个N/4点的DFT。按照此规律，最小变换的点数即所谓的“基数（radix）。”因此，基数为2的FFT算法的最小变换（或称蝶形）是2点DFT。一般地，对N点FFT，对应于N个输入样值，有N个频域样值与之对应。一般而言，FFT算法可以分为时间抽取（DIT）FFT和频率抽取（DIF）两大类。

在实际计算中，可以采用在原来序列后面补0的加长方法来提高FFT的分辨率；可以采用在原来序列后面重复的加长方法来增加FFT的幅度。

24 实验内容

1．软件仿真实验：分别观察并记录正弦序列、方波序列及改变FFT的点数后的频谱；分别观察并记录正弦序列、方波序列及2种加长序列等信号的频谱。

2．硬件实验：分别观察并记录正弦信号、方波信号及改变FFT的点数后的频谱。

25 MATLAB参考程序和仿真内容

%%

function[x]=ffts(mode,M)

Nfft=2^M;

x=zeros(1,Nfft); %定义一个长度为Nfft的一维全0数组

if mode= =1 for n=0:Nfft-1 x(n+1)=sin(2pin/Nfft); end

end %定义一个长度为Nfft的单周期正弦序列

if mode= =2 for n=0:Nfft-1 x(n+1)=sin(4pin/Nfft); end

end %定义一个长度为Nfft的双周期正弦序列

if mode= =3 for n=0:Nfft/2-1 x(n+1)=sin(4pin/Nfft); end

end %定义一个长度为Nfft/2的正弦序列，后面一半为0序列。

if mode= =4 for n=0:Nfft-1 x(n+1)=square(2pin/Nfft); end

end

if mode= =5 for n=0:Nfft-1 x(n+1)=square(2pin/Nfft); end

end

if mode= =6 for n=0:Nfft/2-1 x(n+1)=square(4pin/Nfft); end

end

n=0:Nfft-1;

subplot(2,1,1);

stem(n,x);

axis([0 Nfft-1 11min(x) 11max(x)]);

xlabel('Points-->');

ylabel('x(n)');

X=abs(fft(x,Nfft));

subplot(2,1,2);

stem(n,X);

axis([0 Nfft-1 11min(X) 11max(X)]);

xlabel('frequency-->');

ylabel('!X(k)!');

%%

假设需观察方波信号的频谱，对一个周期的方波信号作32点的FFT，则只需在MATLAB的命令窗口下键入：[x]=ffts(21,5) ，程序进行模拟，并且输出FFT的结果。

关于软件仿真实验内容，建议在完成大量仿真例子的基础上，选择能够体现实验要求的4个以上的例子进行记录。例如要观察后面补0的加长方法来提高FFT的分辨率的现象，可以仿真ffts(4,5)和ffts(6,6)两个例子。

26 硬件实验步骤

1．将低频信号源输出加到实验箱模拟通道1输入端，将示波器探头接至模拟通道1输出端。

2．在保证实验箱正确加电且串口电缆连接正常的情况下，运行数字信号处理与DSP应用实验开发软件，在“数字信号处理实验”菜单下选择“FFT频谱分析”子菜单，出现显示FFT频谱分析功能提示信息的窗口。

3．用低频信号产生器产生一个1KHz的正弦信号。

4．选择FFT频谱分析与显示的点数为64点，开始进行FFT运算。此后，计算机将周期性地取回DSP运算后的FFT数据并绘图显示

5．改信号源频率，观察并记录频谱图的变化。

6．选择FFT的点数为128点，观察并记录频谱图的变化。

7．更改正弦信号的频率，重复步骤4 ~步骤6。

8．用低频信号产生器产生一个1KHz的方波信号，重复步骤4 ~步骤7。注意：应根据实验箱采样频率fs为10KHz和方波信号的频带宽度选择方波信号的频率。

本硬件实验要进行两种信号，每个信号两种频率，每个信号两种点数等共8次具体实验内容，性质能够体现实验要求的4个以上的例子进行记录。

27 思考题

1．对同一个信号，不同点数FFT观察到的频谱图有何区别？

2．序列加长后FFT与原序列FFT的关系是什么，试推导其中一种关系。

3．用傅立叶级数理论，试说明正弦信号频谱和方波信号频谱之间的关系。

3 IIR滤波器设计实验

31 实验目的

1．通过实验加深对IIR滤波器基本原理的理解。

2．学习编写IIR滤波器的MATLAB仿真程序。

32 实验仪器

1．YBLD智能综合信号源测试仪 1台

2．双踪示波器 1台

3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台

4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台

33 实验原理

IIR滤波器有以下几个特点：

1．IIR数字滤波器的系统函数可以写成封闭函数的形式。

2．IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。

3．IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。

4．IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络。

在MATLAB下设计IIR滤波器可使用Butterworth函数设计出巴特沃斯滤波器，使用Cheby1函数设计出契比雪夫I型滤波器，使用Cheby2设计出契比雪夫II型滤波器，使用ellipord函数设计出椭圆滤波器。下面主要介绍前两个函数的使用。

与FIR滤波器的设计不同，IIR滤波器设计时的阶数不是由设计者指定，而是根据设计者输入的各个滤波器参数（截止频率、通带滤纹、阻带衰减等），由软件设计出满足这些参数的最低滤波器阶数。在MATLAB下设计不同类型IIR滤波器均有与之对应的函数用于阶数的选择。

一、巴特沃斯IIR滤波器的设计

在MATLAB下，设计巴特沃斯IIR滤波器可使用butter函数。

Butter函数可设计低通、高通、带通和带阻的数字和模拟IIR滤波器，其特性为使通带内的幅度响应最大限度地平坦，但同时损失截止频率处的下降斜度。在期望通带平滑的情况下，可使用butter函数。

butter函数的用法为：

[b,a]=butter(n,Wn,/ftype/)

其中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n，同时给出对应的截止频率Wn。buttord函数的用法为：

[n,Wn]= buttord(Wp,Ws,Rp,Rs)

其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。

不同类型（高通、低通、带通和带阻）滤波器对应的Wp和Ws值遵循以下规则：

1．高通滤波器：Wp和Ws为一元矢量且Wp>Ws；

2．低通滤波器：Wp和Ws为一元矢量且Wp<Ws；

3．带通滤波器：Wp和Ws为二元矢量且Wp<Ws，如Wp=[02,07],Ws=[01,08];

4．带阻滤波器：Wp和Ws为二元矢量且Wp>Ws，如Wp=[01,08],Ws=[02,07]。

二、契比雪夫I型IIR滤波器的设计

在期望通带下降斜率大的场合，应使用椭圆滤波器或契比雪夫滤波器。在MATLAB下可使用cheby1函数设计出契比雪夫I型IIR滤波器。

cheby1函数可设计低通、高通、带通和带阻契比雪夫I型滤IIR波器，其通带内为等波纹，阻带内为单调。契比雪夫I型的下降斜度比II型大，但其代价是通带内波纹较大。

cheby1函数的用法为：

[b,a]=cheby1(n,Rp,Wn,/ftype/)

在使用cheby1函数设计IIR滤波器之前，可使用cheblord函数求出滤波器阶数n和截止频率Wn。cheblord函数可在给定滤波器性能的情况下，选择契比雪夫I型滤波器的最小阶和截止频率Wn。

cheblord函数的用法为：

[n,Wn]=cheblord(Wp,Ws,Rp,Rs)

其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。

34 实验内容

1．软件仿真实验：编写并调试MATLAB程序，选择不同形式，不同类型的4种滤波器进行仿真，记录幅频和相频特性，对比巴特沃斯滤波器和契比雪夫滤波器。

2．硬件实验：设计IIR滤波器，在计算机上观察冲激响应、幅频特性和相频特性，然后下载到实验箱。用示波器观察输入输出波形，测试滤波器的幅频响应特性。

35 MATLAB参考程序和仿真内容

%%

%mode: 1--巴特沃斯低通；2--巴特沃斯高通；3--巴特沃斯带通；4--巴特沃斯带阻

% 5--契比雪夫低通；6--契比雪夫高通；7--契比雪夫带通；8--契比雪夫带阻

%fp1,fp2： 通带截止频率，当高通或低通时只有fp1有效

%fs1, fs2： 阻带截止频率，当高通或低通时只有fs1有效

%rp: 通带波纹系数

%as: 阻带衰减系数

%sample: 采样率

%h: 返回设计好的滤波器系数

%%

function[b,a]=iirfilt(mode,fp1,fp2,fs1,fs2,rp,as,sample)

wp1=2fp1/sample;wp2=2fp2/sample;

ws1=2fs1/sample;ws2=2fs2/sample;

%得到巴特沃斯滤波器的最小阶数N和3bd频率wn

if mode<3[N,wn]=buttord(wp1,ws1,rp,as);

elseif mode<5[N,wn]=buttord([wp1 wp2],[ws1 ws2],rp,as);

%得到契比雪夫滤波器的最小阶数N和3bd频率wn

elseif mode<7[N,wn]=cheb1ord(wp1,ws1,rp,as);

else[N,wn]=cheblord([wp1 wp2],[ws1 ws2],rp,as);

end

%得到滤波器系数的分子b和分母a

if mode= =1[b,a]=butter(N,wn);end

if mode= =2[b,a]=butter(N,wn,/high/);end

if mode= =3[b,a]=butter(N,wn);end

if mode= =4[b,a]=butter(N,wn,/stop/);end

if mode= =5[b,a]=cheby1(N,rp,wn);end

if mode= =6[b,a]=cheby1(N,rp,wn,/high/);end

if mode= =7[b,a]=cheby1(N,rp,wn);end

if mode= =8[b,a]=cheby1(N,rp,wn,/stop/);end

set(gcf,/menubar/,menubar);

freq_response=freqz(b,a);

magnitude=20log10(abs(freq_response));

m=0:511;

f=msample/(2511);

subplot(3,1,1);plot(f,magnitude);grid; %幅频特性

axis([0 sample/2 11min(magnitude) 11max(magnitude)]);

ylabel('Magnitude');xlabel('Frequency-->');

phase=angle(freq_response);

subplot(3,1,2);plot(f,phase);grid; %相频特性

axis([0 sample/2 11min(phase) 11max(phase)]);

ylabel('Phase');xlabel('Frequency-->');

h=impz(b,a,32); %32点的单位函数响应

t=1:32;

subplot(3,1,3);stem(t,h);grid;

axis([0 32 12min(h) 11max(h)]);

ylabel('h(n)');xlabel('n-->');

%%

假设需设计一个巴特沃斯低通IIR滤波器，通带截止频率为2KHz，阻带截止频率为3KHz，通带波纹系数为1，阻带衰减系数为20，采样频率为10KHz，则只需在MATLAB的命令窗口下键入：

[b,a]=iirfilt(1,2000,3000,2400,2600,1,20,10000)

程序进行模拟，并且按照如下顺序输出数字滤波器系统函数

的系数

b= b0 b1 ……bn

a= a0 a1 ……an

关于软件仿真实验内容，建议在完成大量仿真例子的基础上，选择能够体现实验要求的4个例子进行记录，系统函数只要记录系统的阶数。

36 硬件实验步骤

1．根据实验箱采样频率fs为10KHz的条件，用低频信号发生器产生一个频率合适的低频正弦信号，将其加到实验箱模拟通道1输入端，将示波器通道1探头接至模拟通道1输入端，通道2探头接至模拟通道2输出端。

2．在保证实验箱正确加电且串口电缆连接正常的情况下，运行数字信号处理与DSP应用实验开发软件，在“数字信号处理实验”菜单下选择“IIR滤波器”子菜单，出现提示信息。

3．输入滤波器类型、滤波器截止频率等参数后，分别点击“幅频特性”和“相频特性”按钮，在窗口右侧观察IIR滤波器的幅频特性和相频特性。此时提示信息将消失，如需查看提示信息，可点击“设计说明”按钮。

4．点击“下载实现”按钮，IIR滤波器开始工作，此时窗口右侧将显示IIR滤波器的幅频特性。

5．根据输入滤波器类型，更改低频信号源的频率，观察示波器上输入输出波形幅度的变化情况，测量IIR滤波器的幅频响应特性，看其是否与设计的幅频特性一致。

6．更改滤波器类型、滤波器截止频率等参数（共4种），重复步骤3至步骤5。所选择的例子参数最好和MATLAB仿真程序的例子一样。

7．用低频信号产生器产生一个500Hz的方波信号，分别设计3种滤波器，完成如下表要求的功能，并且记录参数和波形。

功 能

滤波器类型

参 数

输出波形

fp1

fp2

fs1

fs2

通过3次及以下次数的谐波

另外记录图形，并标图号

滤除5次及以下次数的谐波

通过3次到5次的谐波

37 思考题

1．在实验箱采样频率fs固定为10KHz的条件下，要观察方波信号频带宽度内的各个谐波分量，方波信号的频率最高不能超过多少，为什么？

2．硬件实验内容7中输出信号各个谐波分量，与原来方波信号同样谐波分量相比，有没有发生失真？主要发生了什么类型的失真？为什么？

4 窗函数法FIR滤波器设计实验

41 实验目的

1．通过实验加深对FIR滤波器基本原理的理解。

2．学习使用窗函数法设计FIR滤波器，了解窗函数的形式和长度对滤波器性能的影响。

42 实验仪器

1．YBLD智能综合信号源测试仪 1台

2．双踪示波器 1台

3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台

4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台

43 实验原理

数字滤波器的设计是数字信号处理中的一个重要内容。数字滤波器设计包括FIR（有限单位脉冲响应）滤波器与IIR（无限单位脉冲响应）滤波器两种。

与IIR滤波器相比，FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到严格的线性相位特性。设FIR滤波器单位脉冲响应h(n)长度为N，其系统函数H（z）为：

H（z）是z－1的N－1次多项式，它在z平面上有N－1个零点，原点z=0是N－1阶重极点，因此H（z）是永远稳定的。稳定和线性相位特性是FIR滤波器突出的优点。

FIR滤波器的设计任务是选择有限长度的h(n)。使传输函数H( )满足技术要求。FIR滤波器的设计方法有多种，如窗函数法、频率采样法及其它各种优化设计方法，本实验介绍窗函数法的FIR滤波器设计。

窗函数法是使用矩形窗、三角窗、巴特利特窗、汉明窗、汉宁窗和布莱克曼窗等设计出标准响应的高通、低通、带通和带阻FIR滤波器。

一、firl函数的使用

在MATLAB下设计标准响应FIR滤波器可使用firl函数。firl函数以经典方法实现加窗线性相位FIR滤波器设计，它可以设计出标准的低通、带通、高通和带阻滤波器。firl函数的用法为：

b=firl(n,Wn,/ftype/,Window)

各个参数的含义如下：

b—滤波器系数。对于一个n阶的FIR滤波器，其n+1个滤波器系数可表示为：b(z)=b(1)+b(2)z－1+…+b(n+1)z－n。

n—滤波器阶数。

Wn—截止频率，0≤Wn≤1，Wn=1对应于采样频率的一半。当设计带通和带阻滤波器时，Wn=[W1 W2],W1≤ω≤W2。

ftype—当指定ftype时，可设计高通和带阻滤波器。Ftype=high时，设计高通FIR滤波器；ftype=stop时设计带阻FIR滤波器。低通和带通FIR滤波器无需输入ftype参数。

Window—窗函数。窗函数的长度应等于FIR滤波器系数个数，即阶数n+1。

二、窗函数的使用

在MATLAB下，这些窗函数分别为：

1．矩形窗：w=boxcar(n)，产生一个n点的矩形窗函数。

2．三角窗：w=triang(n)，产生一个n点的三角窗函数。

当n为奇数时，三角窗系数为w(k)=

当n为偶数时，三角窗系数为w(k)=

3．巴特利特窗：w=Bartlett(n)，产生一个n点的巴特利特窗函数。

巴特利特窗系数为w(k)=

巴特利特窗与三角窗非常相似。巴特利特窗在取样点1和n上总以零结束，而三角窗在这些点上并不为零。实际上，当n为奇数时bartlett(n)的中心n－2个点等效于triang(n－2)。

4．汉明窗：w=hamming(n)，产生一个n点的汉明窗函数。

汉明窗系数为w(k+1)=054－046cos（ ） k=0,…,n－1

5．汉宁窗：w=hanning(n)，产生一个n点的汉宁窗函数。

汉宁窗系数为w(k)=05[1－cos( )] k=1,…,n

6．布莱克曼窗：w=Blackman(n),产生一个n点的布莱克曼窗函数。

布莱克曼窗系数为w(k)=042－05cos(2π )+08cos（4π ）] k=1,…,n

与等长度的汉明窗和汉宁窗相比，布莱克曼窗的主瓣稍宽，旁瓣稍低。

7．凯泽窗：w=Kaiser(n,beta)，产生一个n点的凯泽窗数，其中beta为影响窗函数旁瓣的β参数，其最小的旁瓣抑制α与β的关系为:

01102(α－087) α>50

β= 05842(α－21)04+007886(α－21) 21≤α≤50

0 α<21

增加β可使主瓣变宽，旁瓣的幅度降低。

8．契比雪夫窗：w=chebwin(n,r)产生一个n点的契比雪夫窗函数。其傅里叶变换后的旁瓣波纹低于主瓣r个db数。

44 实验内容

1．软件仿真实验：编写并调试MATLAB程序，观察不同窗，不同类型滤波器不同点数等共4种FIR滤波器的h(n)，并记录幅频特性和相频特性。

2．硬件实验：用窗函数法设计标准响应的FIR滤波器，在计算机上观察窗函数幅频特性、幅频特性和相频特性，然后下载到实验箱。用示波器观察输入输出波形，测试滤波器的幅频响应特性。

45 MATLAB参考程序和仿真内容

%%

%mode: 模式（1--高通；2--低通；3--带通；4--带阻）

%n: 阶数，加窗的点数为阶数加1

%fp: 高通和低通时指示截止频率，带通和带阻时指示下限频率

%fs: 带通和带阻时指示上限频率

%window:加窗（1--矩形窗；2--三角窗；3--巴特利特窗；4--汉明窗；

% 5--汉宁窗；6--布莱克曼窗；7--凯泽窗；8--契比雪夫窗）

%r: 代表加chebyshev窗的r值和加kaiser窗时的beta值

%sample: 采样率

%h: 返回设计好的FIR滤波器系数

%%

%mode: 模式（1--高通；2--低通；3--带通；4--带阻）

%n: 阶数，加窗的点数为阶数加1

%fp: 高通和低通时指示截止频率，带通和带阻时指示下限频率

%fs:

1 buttord

（1）[N,wc]=buttord(wp，ws，αp，αs)

用于计算巴特沃斯数字滤波器的阶数N和3dB截止频率wc。

调用参数wp，ws分别为数字滤波器的通带、阻带截止频率的归一化值，要求：0≤wp≤1，0≤ws≤1。1表示数字频率pi。

αp，αs分别为通带最大衰减和组带最小衰减（dB）。

当ws≤wp时，为高通滤波器；

当wp和ws为二元矢量时，为带通或带阻滤波器，这时wc也是二元向量。

N，wc作为butter函数的调用参数。

（2）[N,Ωc]=buttord(Ωp，Ωs，αp，αs，‘s’)

用于计算巴特沃斯模拟滤波器的阶数N和3dB截止频率Ωc。

Ωp，Ωs，Ωc均为实际模拟角频率。

说明：buttord函数使用阻带指标计算3dB截止频率，这样阻带会刚好满足要求，而通带会有富余。

2buttap（N）

[z0，p0，k0]=buttap(N)

用于计算N阶巴特沃斯归一化（3dB截止频率Ωc=1）模拟低通原型滤波器系统函数的零、极点和增益因子。

说明：如果要从零、极点模型得到系统函数的分子、分母多项式系数向量ba、aa，可调用

[B，A]=zp2tf（z0，p0，k0）

3butter

（1）[b，a]=butter（N，wc，‘ftype’）

计算N阶巴特沃斯数字滤波器系统函数分子、分母多项式的系数向量b、a。

调用参数N和wc分别为巴特沃斯数字滤波器的阶数和3dB截止频率的归一化值（关于pi归一化），一般是调用buttord（1）格式计算N和wc。

系数b、a是按照z-1的升幂排列。

（2）[B，A]=butter（N，Ωc，‘ftype’，‘s’）

计算巴特沃斯模拟滤波器系统函数的分子、分母多项式系数向量ba、aa。

调用参数N和Ωc分别为巴特沃斯模拟滤波器的阶数和3dB截止频率(实际角频率)，可调用buttord（2）格式计算N和Ωc。

系数B、A按s的正降幂排列。

tfype为滤波器的类型：

◇ftype=high时，高通；Ωc只有1个值。

◇ftype=stop时，带阻阻；此时Ωc=[Ωcl,Ωcu]，分别为带阻滤波器的通带3dB下截止频率和上截止频率。

◇ ftype缺省时：

若Ωc只有1个值，则默认为低通；

若Ωc有2个值，则默认为带通；其通带频率区间Ωcl < Ω < Ωcu。

注意：所设计的带通和带阻滤波器系统函数是2N阶。因为带通滤波器相当于N阶低通滤波器与N阶高通滤波器级联。

相关文章:

数学思想及理论

均匀乱数

MDSC很可能,我已经用尽了现有原始数据的所有潜能

是可以用butter函数直接设计带通的，

[b,a]=BUTTER(N,wn)只是此时的wn=[w1,w2];

对N和wn可用[N,wn]=buttord(wp,ws,Rp,As)此时的wp=[wp1 wp2] 和ws=[ws1 ws2]

PAUSE Wait for user response PAUSE等待用户的相应。

PAUSE(n) pauses for n seconds before continuing, where n can also be a fraction The resolution（分辨率） of the clock is platform specific Fractional pauses of 001 seconds should be supported on most platforms

PAUSE(n)函数是程序停止n秒后继续，n可以使分数或小数。函数的分辨率根据平台而定，大多数平台都能满足001秒的精度。

PAUSE causes a procedure to stop and wait for the user to strike any key before continuing

可以使用PAUSE让程序停止运行直到用户敲击键盘。

PAUSE OFF indicates that any subsequent PAUSE or PAUSE(n) commands should not actually pause This allows normally interactive scripts to run unattended

PAUSE OFF 关闭之后所有的PAUSE函数的暂停功能。

PAUSE ON indicates that subsequent PAUSE commands should pause

PAUSE ON 打开之后的所有PAUSE函数的功能。

MATLAB提供的用于图形控制的函数和命令：

axis: 人工选择坐标轴尺寸

clf：清图形窗口

ginput： 利用鼠标的十字准线输入

hold： 保持图形

shg：显示图形窗口

subplot： 将图形窗口分成N块子窗口。

1．图形窗口(figure window)

(1) 图形窗口的创建和选择(Creating and selecting of figure window)

figure(n)：用于为当前的绘图创建图形窗口，每运行一次figure就会创建一个新的图形窗口,n表示第n个

窗口,如果窗口定义了句柄，也可以用figure(h)将句柄h的窗口作为当前窗口。

clf ：用于清除当前图形窗口中的内容。

shg ：用于显示当前图形窗口。

(2) 在一个图形窗口中绘制多个子图形(Drawing several subfigures in a single window)

subplot(m,n,p)：把窗口分成m×n个小窗口，并把第p个窗口当作当前窗口。

例：将4 个图形显示在同一个图形窗口中。

t=0:pi/20:2pi; [x,y]=meshgrid(t);

subplot(2,2,1); plot(sin(t),cos(t)); axis equal

subplot(2,2,2); z=sin(x)+cos(y); plot(t,z); axis([0 2pi –2 2])

subplot(2,2,3); z=sin(x)cos(y); plot(t,z); axis([0 2pi –1 1])

subplot(2,2,4); z=sin(x)^2-cos(y)^2; plot(t,z); axis([0 2pi –1 1])

(3) 在一个已有的图形上绘图(Drawing a figure on the figure was existed)

hold on ：在一个已有的图形上继续绘图；

hold off： 命令结束继续绘图。

例：将peaks函数的等高线图与伪彩色画在一起。

[x,y,z]=peaks; %产生双变量数组

contour(x,y,z,20,'k') %绘制等高线

hold on

pcolor(x,y,z) %绘制伪彩色图

shading interp %表面色彩渲染

hold off

2．坐标轴控制命令(Axis control commands)

控制坐标性质的axis函数的多种调用格式：

axis(xmin xmax ymin ymax) ：指定二维图形x和y轴的刻度范围,

axis auto ：设置坐标轴为自动刻度（缺省值）

axis manual（或axis(axis)）：保持刻度不随数据的大小而变化

axis tight ：以数据的大小为坐标轴的范围

axis ij ：设置坐标轴的原点在左上角，i为纵坐标，j为横坐标

axis xy ：使坐标轴回到直角坐标系

axis equal ：使坐标轴刻度增量相同

axis square ：使各坐标轴长度相同，但刻度增量未必相同

axis normal ：自动调节轴与数据的外表比例，使其他设置失效

axis off ：使坐标轴消隐

axis on ：显现坐标轴

(1) 坐标轴的范围(Domain of coordinates axis)

二维图形坐标轴范围在缺省状态下是根据数据的大小自动设置的，如欲改变，可利用axis(xmin xmax ymin ymax),函数来定义。

例： 定义坐标轴范围对观察图形的影响。

x=0:01:pi/2; figure(1); plot(x,tan(x),'-ro') %ymax=tan(157),而其他数据都很小，结果将

%使图形难于进行观察和判断。

figure(2); plot(x,tan(x),'-ro'); axis([0, pi/2,0,5]) %对坐标轴的范围进行控制就可得到较满意的绘图结果

(2) 显示比例对绘图结果的影响(Effect of display scaling on plotting results)

例：比较(Default, axis square, axis equal, axis tight)几种不同的显示方式的显示效果。

t=0:pi/20:2pi; figure(1);

subplot(2,1,1); plot(sin(t),2cos(t)); grid on %缺省状态下的图形比例

subplot(2,1,2); plot(sin(t),2cos(t)); axis square; grid on %正方形的显示比例

figure(2)

subplot(1,2,1); plot(sin(t),2cos(t)) ; axis equal; grid on %具有相等的刻度比例

subplot(1,2,2); plot(sin(t),2cos(t)); axis tight ; grid on %紧缩形式

3．图形标注(Marking on the figure)：MATLAB的图形标注方法（表 6—7）

title ：标题,

xlabel ：x轴标注,

ylabel ：y轴标注,

text ：任意定位的标注

gtext ：鼠标定位标注,

legent ：标注图例

图形标注可以使用字母，数字，汉字或按规定的方法表示希腊字母。如：\pi表示π，\leq表示≤，\rm表示后面的字恢复为正体字，\it表示斜体字，FontSize表示字体的大小， FontName表示字体的类型等。

可以使用图形窗口的Insert菜单，也可以使用属性编辑器，还可以使用函数输入的方法加标注,以下介绍相关函数的使用方法。

(1) 加注坐标轴标识和图形标题(Add axis labels and title of figure)

加注坐标轴标识：xlabel(‘s’), ylabel(‘s’)

图形标题: title(‘s’)

例：加注坐标轴标示和图形标题。

t=0:pi/100:2pi;y=sin(t);

plot(t,y)

axis([0 2pi,-1 1])

xlabel('0 \leq \itt \rm \leq \pi','FontSize',16)

ylabel('sin(t)','FontSize',20)

title('正弦函数图形','FontName','隶书','FontSize',20)

(2) 图中加注文本(Add text in the figure)

text(x,y,’字符串’)

例：在上图中加语句。

t=0:pi/100:2pi;

y=sin(t);

plot(t,y)

axis([0 2pi,-1 1])

xlabel('0 \leq \itt \rm \leq \pi','FontSize',16)

ylabel('sin(t)','FontSize',20)

title('正弦函数图形','FontName','隶书','FontSize',20)

text(3pi/4,sin(3pi/4),'\leftarrowsin(t)=0707', 'FontSize',16)

text(pi,sin(pi),'\leftarrowsin(t)=0', 'FontSize',16)

text(5pi/4,sin(5pi/4),'sin(t)=-0707\rightarrow','FontSize',16,

'HorizontalAlignment','right')

句中：

leftarrow 表示加一个向左的箭头

rightarrow 表示加一个向右的箭头

HorizontalAlignment 表示右对齐水平排列

gtext('字符串'): 在图形窗口上用鼠标直接在指定的位置上加注文本。

例：

t=0:pi/100:2pi;

y=sin(t);

plot(t,y)

axis([0 2pi,-1 1])

xlabel('0 \leq \itt \rm \leq \pi','FontSize',16)

ylabel('sin(t)','FontSize',20)

title('正弦函数图形','FontName','隶书','FontSize',20)

gtext('MATLAB')

(3) 指定TeX字符

例：在标题中指定TeX字符

t=0:pi/100:2pi;

alpha=-08;

beta=15;

y=sin(betat)exp(alphat);

plot(t,y)

title('{\itAe}^{-\it\alpha\itt}sin\it\beta{\itt}\it\alpha<<\it\beta')

xlabel('时间\mus'),

ylabel('幅值')

在title中的字符串表现的是 Aeαt sinβt α<<β

斜体Ae 上标斜体αt 斜体βt 斜体α 斜体β

(4) 在图形中添加图例框(Add legend in the figure)

legend(字符串1，字符串2，…)

例：在当前图形中添加图例说明。

x=0:pi/10:2pi;

y1=sin(x);

y2=06sin(x);

y3=03sin(x);

plot(x,y1,x,y2,'-o',x,y3,'-')

legend( '曲线1','曲线2','曲线3')

legend('boxoff')

legend函数的其他功能见（表 6—8）

4 图线形式和颜色(Style and color of plot)

(1) 图线的形式: (style of plot)

四种线形: 实线'-'，虚线'--'， 点线':'，点划线'-'

标记点类型：点'', 圆'o', 加号'+', 星号'', x符号'x', 方形's', 菱形'd', hexagram 'h'

上三角△ '^', 下三角▽ 'v', 左三角'<', 右三角'>', 正五边形'p'

命令：plot(x,y,'—'), plot(x1,y1,':’,x2,y2,'')

例1：选择不同的线形绘图。

t=0:pi/100:2pi; y=sin(t); y2=sin(t-025); y3=sin(t-05);

plot(t,y,'-',t,y2,'-',t,y3,':')

例2：选择不同的标记点绘图。

t=0:pi/20:2pi; x=t^3; y=sin(t); plot(x,y,'o')

(2) 线的颜色(color of plot)

可选颜色： 红r，绿g, 蓝b, 黄y, 粉红m, 青c, 黑k

例：t=0:pi/20:2pi;

y=sin(t); plot(x,y,'r'), plot(x,y,'g+')

(3) 图线的其他属性(other characters of plot)

设置图线的宽度 : 'LineWidth'

标记点的边缘颜色: 'MarkerEdgeColor'

填充颜色 : 'MarkerFaceColor'

标记点的大小 : 'MarkerSize'

例： 设置图线的线形、颜色、宽度、标记点的颜色及大小。

t=0:pi/20:pi; y=sin(4t)sin(t)/2;

plot(t,y,'-bs','LineWidth',2,'MarkerEdgeColor','k', 'MarkerFaceColor', 'y','MarkerSize',10);

首先要明白，由你给的代码知道你是要设计数字低通滤波器，而你给的频率是模拟的，你需要将模拟的频率转换到数字，butter函数是用在数字中的函数，Wn应在0~1之间，即你写的2pi500是不可以的。

具体程序可以这样写（去采样周期为T=1s,采样频率FS=1/T=1hz）

（1）设通带频率为400hz，截止频率为500hz，通带衰减Rp=1db，阻带衰减Rs=20db

（2）用双线性变换设计时，模拟和数字的频率对应关系为W（模拟）=2/Ttan(1/2)w(数字)，得到数字的技术指标为Wpz，Wsz。

（3）设计滤波器

[Nd,Wdc]=buttor(Wpz,Wsz,Rp,Rs);

[B,A]=butter(Nd,Wdz);

(4)滤波

z=filter(B,A,要滤波的函数)