基于MATLAB的频谱分析与信号去噪汇总

学生毕业设计报告

基于MATLAB的频谱分析与信号去噪

作 者 专 业

系（院） 年 级

学 号 日 期

指导教师

学生诚信承诺书

本人郑重承诺：所呈交的设计报告是我个人在导师指导下进行的研究工作及取得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，报告中不包含其他人已经发表或撰写的研究成果，也不包含为获得安阳师范学院或其他教育机构的学位或证书所使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在报告中作了明确的说明并表示了谢意。

签名： 日期：

报告使用授权说明

本人完全了解 有关保留、使用学位报告的规定，即：学校有权保留送交报告的复印件，允许报告被查阅和借阅；学校可以公布报告的全部或部分内容，可以采用影印、缩印或其他复制手段保存报告。

签名： 导师签名： 日期：

基于MATLAB的频谱分析与信号去噪

摘要：本课题是基于数字信号处理的理论知识对语音信号、图像信号等的频谱分析以及对加噪声语音信号进行时域、频域分析和滤波设计。然后利用 MATLAB软件进行编程，调试并完善程序，最终在计算机上得以实现。滤波器设计在数字信号处理中占有极其重要的地位，数字滤波器又有FIR数字滤波器和IIR滤波器两种。利用MATLAB信号处理工具箱可以快速有效地设计各种数字滤波器。在设计实现的过程中，使用窗函数法来设计FIR数字滤波器，用巴特沃斯、切比雪夫和双线性变法设计IIR数字滤波器，并利用MATLAB 作为辅助工具完成设计中的计算与图形的绘制。通过对对所设计滤波器的仿真和频率特性分析，可知利用MATLAB信号处理工具箱可以有效快捷地设计FIR和IIR数字滤波器，过程简单方便，结果的各项性能指标均符合指定要求。

关键词 频域分析， 滤波器 ， MATLAB 1 引言 1.1 课题背景

随着信息时代和数字世界的到来，数字信号处理己成为当今一门极其重要的学科和技术领域，数字信号处理在通信、语音、图像、自动控制、医疗和家用电器等众多领域得到了广泛的应用。任意一个信号都具有时域与频域特性，信号的频谱完全代表了信号，因而研究信号的频谱就等于研究信号本身。通常从频域角度对信号进行分析与处理，容易对信号的特性获得深入的了解。因此，信号的频谱分析是数字信号处理技术中的一种较为重要的工具。

数字滤波器, 是数字信号处理中及其重要的一部分。随着数字技术的发展，受到人们越来越多的重视。数字滤波器可以通过数值运算实现滤波，所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活、不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊功能。数字滤波器种类很多，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种，即有限冲激响应( FIR，Finite Impulse Response)滤波器和无限冲激响应( IIR，Infinite Impulse Response)滤波器。

在工程领域中，MATLAB是一种倍受程序开发人员青睐的语言，对于一些需要做大量数据运算处理的复杂应用以及某些复杂的频谱分析算法MATLAB显得游刃有余。 1.2 研究意义

信号处理几乎涉及到所有的工程技术领域，而频谱分析又是信号处理中一个非常重要的分析手段。一般的频谱分析都依靠传统频谱分析仪来完成，价格昂贵，体积庞大，不便

第 1 页

于工程技术人员的携带。虚拟频谱分析仪改变了原有频谱分析仪的整体设计思路，用软件代替了硬件，使工程技术人员可以用一部笔记本电脑到现场就可轻松完成信号的采集、处理及频谱分析。

信号去噪被用于从一堆波音资料中提取有用信息，去除干扰，提高波音资料信噪比。为了提高信噪比，人们根据信号和噪声的各种特征差异，设计了许多去噪方法，并在应用中并取得了很好的成果。信号去噪的很多方法都是利用短时傅立叶变换来滤波去噪，但是短时傅立叶变换不能同时兼顾时间分辨率和频率分辨率。 用不同种滤波器滤波也是一种有效可行的方法。 1.3 本文研究内容

信号的频谱分析就是利用傅里叶分析的方法，求出与时域描述相对应的频域描述,从中找出信号频谱的变化规律，以达到特征提取的目的。不同信号的傅里叶分析理论与方法，在有关专业书中都有介绍。

但实际的待分析信号一般没有解析式，直接利用公式进行傅里叶分析非常困难.。DFT是一种时域和频域均离散化的傅里叶变换，适合数值计算且有快速算法，是分析信号的有力工具。

DFT及FFT是数字信号处理的重要内容。DFT是FFT的基础，FFT是DFT的快速算法，在 MAT LAB中可以利用函数 FFT来计算序列的离散傅里叶变换DFT。基于此首先介绍了MAT LAB软件；其次给出了基于MAT LAB软件的DFT和 FFT频谱分析的方法，利用 MAT LAB软件方法，使得设计方便、快捷，大大减轻了工作量；再而我们提取一段语音信号，加噪，再通过设计的不同滤波器进行滤波后的频谱分析。本文将重点介绍基于MATLAB的频谱分析设计，包括：

(1) 音频信号频谱分析； (2) 图像信号频谱分析； (3) 离散信号/序列频谱分析； (4)语音信号提取，分析和加噪； (5)设计IIR和FIR的各种滤波器； (6)用设计的滤波器进行滤波； (7)分析滤波前后信号特征。 2 频谱分析技术及MATLAB简介 2.1 时域抽样定理

时域抽样定理给出了连续信号抽样过程中信号不失真的约束条件：对于基带信号，信号抽样频率fsam大于等于2倍的信号最高频率fm，即fsam?2fm。时域抽样是把连续信号

X?t?变成适于数字系统处理的离散信号X[k]。对连续信号X?t?以间隔T抽样，则可得到

的离散序列为X[k]?X?kt?t?kT，如图 2-1所示。

第 2 页

x(t)0T2Tx[k]t012k 图2-1 连续信号抽样的离散序列

若X[k]?X?kt?t?kT，则信号X?t?与X[k]的频谱之间存在：

1?X(e)??X?j(??n?sam)?

Tn???j?其中，X(ej?)的频谱为X(j?)，X[k]的频谱为X(ej?)。

可见，信号时域抽样导致信号频谱的周期化。?sam?2?/T(rad/s)为抽样角频率，fsam?1T为抽

样频率。数字角频率Ω与模拟角频率ω的关系为：Ω=ωT。

2.2 离散傅立叶变换(DFT)

有限长序列x(n)的离散傅立叶变换(DFT)为

?knX(k)?DFT[x(n)]??x(n)WN,0?n?N?1

n?0N?1逆变换为

1N?1?knx(n)?IDFT[X(k)]??X(k)WN,0?n?N?1

Nn?02.3 快速傅立叶变换(FFT)

在各种信号序列中，有限长序列占重要地位。对有限长序列可以利用离散傅立叶变换(DFT)进行分析。DFT不但可以很好的反映序列的频谱特性，而且易于用快速算法(FFT)在计算机上进行分析。

有限长序列的DFT是其z变换在单位圆上的等距离采样，或者说是序列傅立叶的等距离采样，因此可以用于序列的谱分析。FFT是DFT的一种快速算法，它是对变换式进行一次次分解，使其成为若干小数据点的组合，从而减少运算量。

MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有fft、ifft、fft2 、ifft2, fftn、ifftn和fftshift、ifftshift等。当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。

fft函数调用方式： (1)Y=fft(X) (2)Y=fft(X,N)

(3)Y＝fft(X,[],dim)或Y＝fft(X,N,dim)。 函数ifft的参数应用与函数fft完全相同。

第 3 页

图3-3 纹理频谱图

4 有噪声的语音信号分析与去噪 4.1 有噪语音信号提取 4.1.1 语音信号的采集

利用PC 机上的声卡和WINDOWS 操作系统可以进行数字信号的采集。将话筒输入计算机的语音输入插口上,启动录音机。按下录音按钮,接着对话筒说话“语音信号处理”,说完后停止录音,屏幕左侧将显示所录声音的长度。点击放音按钮,可以实现所录音的重现。以文件名“speech”保存入C:\\MATLAB6p5\\work 中。可以看到,文件存储器的后缀默认为. wav ,这是WINDOWS 操作系统规定的声音文件存的标准。 4.1.2 语音信号的时频分析

MATLAB软件平台下，利用wavread函数对语音信号进行采样，记住采样频率和采样点数

wavread 函数调用格式

y=wavread(file)%读取file所规定的wav文件，返回采样值放在向量y中。 [y,fs,nbits]=wavread(file) %采样值放在向量y中，fs表示采样频率(hz)，

nbits表示采样位数。

y=wavread(file，N)%读取钱N点的采样值放在向量y中。

y=wavread(file，[N1,N2])%读取从N1到N2点的采样值放在向量y中。 对语音信号speech.wav进行采样其程序如下：

[y,fs,nbits]=wavered ('speech'); %把语音信号进行加载入MATLAB 仿真软

件平台中

fs =

第 9 页

44100 nbits = 16

首先画出语音信号的时域波形，然后对语音信号进行频谱分析。在MATLAB中利用fft对信号进行快速傅里叶变换，得到信号的频谱特性。

其程序如下：

[y,fs,nbits]=wavread ('speech'); sound(y,fs,nbits); %回放语音信号 n = length (y) ; %求出语音信号的长度 Y=fft(y,n); %傅里叶变换

subplot(2,1,1);plot(y);title('原始信号波形'); subplot(2,1,2);plot(abs(Y));title('原始信号频谱')

程序结果如图4-1所示。

原始信号波形10.50-0.5-10x 1045原始信号频谱1015x 10521.510.50051015x 105

图4-1 原始信号特征

4.1.3 语音信号加噪与频谱分析

利用MATLAB中的随机函数(rand或randn)产生噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。

其程序如下：

[y,fs,nbits]=wavread ('speech');

第 10 页

n = length (y) ; %求出语音信号的长度 noise=0.21*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 subplot(2,1,1);

plot(s);title('加噪语音信号的时域波形'); S=fft(s); %傅里叶变换 subplot(2,1,2);

plot(abs(S));title('加噪语音信号的频域波形') 程序结果如图4-2所示。

加噪语音信号的时域波形210-1-20x 10451015x 105加噪语音信号的频域波形21.510.50051015x 105

图4-2 加噪信号特征

4.2 设计FIR和IIR数字滤波器

IIR滤波器和FIR滤波器的设计方法完全不同。IIR滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行的。其设计步骤是：先设计过渡模拟滤波器得到系统函数H(s)，然后将H(s)按某种方法转换成数字滤波器的系统函数H(z)。FIR滤波器采用间接法，常用的方法有窗函数法、频率采样发和切比雪夫等波纹逼近法。对于线性相位滤波器，经常采用FIR滤波器。

对于数字带通滤波器的设计，通用方法为双线性变换法。可以借助于模拟滤波器的频率转换设计一个所需类型的过渡模拟滤波器，再经过双线性变换将其转换策划那个所需的数字滤波器。具体设计步骤如下：

(1)确定所需类型数字滤波器的技术指标。

第 11 页

(2)将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为Ω=2/T tan(0.5ω)

(3)将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。 (4)设计模拟低通滤波器。

(5)通过频率变换将模拟低通转换成相应类型的过渡模拟滤波器。

(6)采用双线性变换法将相应类型的过渡模拟滤波器转换成所需类型的数字滤波器。 我们知道，脉冲响应不变法的主要缺点是会产生频谱混叠现象，使数字滤波器的频响偏离模拟滤波器的频响特性。为了克服之一缺点，可以采用双线性变换法。

下面我们总结一下利用模拟滤波器设计IIR数字低通滤波器的步骤：

(1)确定数字低通滤波器的技术指标：通带边界频率、通带最大衰减，阻带截止频率、阻带最小衰减。

(2)将数字低通滤波器的技术指标转换成相应的模拟低通滤波器的技术指标。 (3)按照模拟低通滤波器的技术指标设计及过渡模拟低通滤波器。

(4)用双线性变换法，模拟滤波器系统函数转换成数字低通滤波器系统函数。 如前所述，IIR滤波器和FIR滤波器的设计方法有很大的区别。下面我们着重介绍用窗函数法设计FIR滤波器的步骤。如下：

(1)根据对阻带衰减及过渡带的指标要求，选择串窗数类型(矩形窗、三角窗、汉宁窗、哈明窗、凯塞窗等)，并估计窗口长度N。先按照阻带衰减选择窗函数类型。原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。

(2)构造希望逼近的频率响应函数。 (3)计算h(n)。

(4)加窗得到设计结果。

接下来，我们根据语音信号的特点给出低通滤波器的性能指标：

fp=1000Hz，fc=1200Hz，As=50db ,Ap=1dB

在MATLAB中，可以利用函数fir1设计FIR滤波器，利用函数butter,cheby1和ellip设计IIR滤波器，利用MATLAB中的函数freqz画出各步步器的频率响应。

hn=fir1(M，wc，window)，可以指定窗函数向量window。如果缺省window参数，则fir1默认为哈明窗。其中可选的窗函数有Rectangular Barlrtt Hamming Hann Blackman窗，其相应的都有实现函数。

MATLAB信号处理工具箱函数buttp buttor butter是巴特沃斯滤波器设计函数，其有5种调用格式，本课程设计中用到的是[N,wc]=butter(N,wc,Rp,As,’s’),该格式用于计算巴特沃斯模拟滤波器的阶数N和3dB截止频率wc。

MATLAB信号处理工具箱函数cheblap,cheblord和cheeby1是切比雪夫I型滤波器设计函数。我们用到的是cheeby1函数，其调用格式如下：

[B,A]=cheby1(N,Rp,wpo,’ftypr’) [B,A]=cheby1(N,Rp,wpo,’ftypr’,’s’)

函数butter,cheby1和ellip设计IIR滤波器时都是默认的双线性变换法，所以在设计滤波器时只需要代入相应的实现函数即可。

第 12 页

下面我们将给出FIR和IIR数字滤波器的主要程序。 IIR低通滤波器程序： Ft=8000; Fp=1000; Fs=1200; wp=2*pi*Fp/Ft; ws=2*pi*Fs/Ft; fp=2*Ft*tan(wp/2); fs=2*Fs*tan(wp/2);

[n11,wn11]=buttord(wp,ws,1,50,'s'); %求低通滤波器的阶数和截止频率 [b11,a11]=butter(n11,wn11,'s'); %求S域的频率响应的参数

[num11,den11]=bilinear(b11,a11,0.5); %利用双线性变换实现频率响应S域到Z域的变换

[h,w]=freqz(num11,den11); %根据参数求出频率响应 plot(w*8000*0.5/pi,abs(h)); legend('用butter设计'); grid

生成相应图片如图4-3所示。

1.4用butter设计1.2 10.80.60.40.20 05001000150020002500300035004000

图4-3 IIR低通滤波器

第 13 页

相应的IIR带通滤波器程序见附录 1 生成相应图片见图4-4所示。

1.5用butter设计 10.50 05001000150020002500300035004000

图4-4 IIR带通滤波器

FIR低通用窗函数设计低通滤波器的程序： Ft=8000; Fp=1000; Fs=1200; wp=2*Fp/Ft; ws=2*Fs/Ft; rp=1; rs=50;

p=1-10.^(-rp/20); %通带阻带波纹 s=10.^(-rs/20); fpts=[wp ws]; mag=[1 0]; dev=[p s];

[n21,wn21,beta,ftype]=kaiserord(fpts,mag,dev);%由kaiserord求滤波器的阶数和截止频率

b21=fir1(n21,wn21,Kaiser(n21+1,beta)); %由fir1设计滤波器

第 14 页

[h,w]=freqz(b21,1); %得到频率响应 plot(w/pi,abs(h)); title('FIR低通滤波器'); grid;

生成相应图片见图4-5所示。

FIR低通滤波器1.41.210.80.60.40.2000.10.20.30.40.50.60.70.80.91

图4-5 FIR低通滤波器

FIR带通滤波器程序见附录 2 生成相应图片见图4-6所示。

第 15 页

1.5用butter设计 10.50 05001000150020002500300035004000

图4-6 FIR带通滤波器

4.3 用滤波器对加噪语音信号进行滤波和分析 4.3.1 滤波器滤波

用自己设计的各滤波器分别对加噪的语音信号进行滤波，在MATLAB中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。

1)双线性变换法： i.低通滤波器

z11=filter(num11,den11,s);

ii.带通滤波器

z12=filter(num12,den12,s);

2)窗函数法 i.低通滤波器 z21=fftfilt(b21,s);

ii.带通滤波器 z22=fftfilt(b22,s);

4.3.2 比较滤波前后语音信号的波形及频谱

(1)双线性变换法

第 16 页

21.510.50x 10滤波前信号的频谱4双线性变4换法低通滤波后信号的频谱x 1021.510.5051015x 10500510155滤波前信号的波形210-1-2051015x 105x 10双线性变换法低通滤波后的信号波形210-1-2051015x 105

图4-7 双线性变换法低通滤波

1)用低通滤波器比较程序： [y,fs,nbits]=wavread ('speech');

n = length (y) ; %求出语音信号的长度 noise=0.21*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 S=fft(s); %傅里叶变换 z11=filter(num11,den11,s);

m11=fft(z11); %求滤波后的信号 subplot(2,2,1); plot(abs(S),'g');

title('滤波前信号的频谱'); grid;

subplot(2,2,2); plot(abs(m11),'r');

title('双线性变换法低通滤波后信号的频谱'); grid;

subplot(2,2,3); plot(s);

title('滤波前信号的波形');

第 17 页

grid;

subplot(2,2,4); plot(z11);

title('双线性变换法低通滤波后的信号波形'); grid;

生成相应的滤波前后信号特征比较见图4-7所示。 2)用带通滤波器比较程序见附录3

生成相应的滤波前后信号特征比较见图4-8所示。

x 10滤波前信号的频谱421.510.50双线性变换法带通滤波后信号的频谱150010005000051015x 1050510155滤波前信号的波形210-1-2051015x 105x 10双线性变换法带通滤波后的信号波形210-1-2051015x 105

图4-8 双线性变换法带通滤波

(2)窗函数法

1)用低通滤波器比较程序： [y,fs,nbits]=wavread ('speech');

n = length (y) ; %求出语音信号的长度 noise=0.21*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 S=fft(s); %傅里叶变换 z21=fftfilt(b21,s);

m21=fft(z21); %求滤波后的信号 subplot(2,2,1);

第 18 页

plot(abs(S),'g');

title('滤波前信号的频谱'); grid;

subplot(2,2,2); plot(abs(m21),'r');

title('窗函数法低通滤波后信号的频谱'); grid;

subplot(2,2,3); plot(s);

title('滤波前信号的波形'); grid;

subplot(2,2,4); plot(z21);

title('窗函数法低通滤波后的信号波形'); grid;

生成相应的滤波前后信号特征比较见图4-9所示。

21.510.50051015x 10滤波前信号的波形210-1-2051015x 1055x 10滤波前信号的频谱4窗函数4法低通滤波后信号的频谱x 1021.510.500510155x 10窗函数法低通滤波后的信号波形210-1-2051015x 105

图4-9 窗函数法低通滤波

2)用带通滤波器比较程序见附录 4

生成相应的滤波前后信号特征比较见图4-10所示。

第 19 页

21.510.50x 10滤波前信号的频谱4窗函数4法带通滤波后信号的频谱x 104321051015x 10500510155滤波前信号的波形210-1-2051015x 105x 10窗函数法带通滤波后的信号波形420-2-4051015x 105

图4-10 窗函数法带通滤波

分析：加入噪声后回放的声音与原始的语音信号有明显的不同,其伴随较尖锐的干扰啸叫声。从含噪语音信号的频谱图中可以看出是含噪的语音信号的频谱,在整个频域范围内分布均匀,这正是干扰所造成的。通过滤波前后的对比，低通滤波后效果最好，带通滤波后的效果相对较差。由此可见，语音信号主要分布在低频，而噪声主要分布在高频。 6 结论

我感受到只有在了解课本知识的前提下，才能更好的应用这个工具；并且熟练的应用MATLAB也可以很好的加深我对课程的理解，为我带来方便。这次设计使我了解了MATLAB的使用方法，学会分析滤波器的优劣和性能，提高了分析和动手实践能力。同时我相信，进一步加强对MATLAB的学习与研究对我今后的学习将会起到很大的帮助。

通过这次设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能真正为社会服务，从而提高自己的实际动手能力和独立思考的能力。在设计的过程中遇到问题，可以说得是困难重重，这毕竟第一次做的，难免会遇到过各种各样的问题，同时在设计的过程中发现了自己的不足之处，对以前所学过的知识理解得不够深刻，掌握得不够牢固。

第 20 页

参考文献

[1]黄植功. 基于MATLAB的连续信号频谱分析[J]. 广西物理, 2009, (03) ，36-52 [2]刘小群,周云波. 基于Matlab的DFT及FFT频谱分析[J]. 山西电子技术, 2010, (04) ，47-58

[3]董静薇,于广艳. 基于Matlab的FFT频谱分析及IIR数字滤波器设计[J]. 软件导刊, 2008, (10) ，36-40

[4]赵淑敏. 基于MATLAB实现对语音信号频谱分析[J]. 江西通信科技, 2010, (01) ，59-63

[5]赵淑敏. 基于MATLAB实现对语音信号频谱分析[J]. 信息通信,2010, (04) ，58-62 [6]曾尚璀,沈华,俞振利. 基于Matlab系统的信号FFT频谱分析与显示[J]. 科技通报, 2000, (04)：103-110

[7]曾秀莲,程保胜,潘雪峰. 连续时间信号频谱分析及MATLAB实现[J]. 中国水运(学术版), 2006, (05) ，93-98

[8]杨峰,苏玉萍,余冬菊. 用MATLAB模拟实现数字信号的调制与频谱分析[J]. 电脑学习, 2008, (04) ，77-81

[9]王占丽. 频谱分析技术 [J]. 黑龙江科技信息, 2009, (27)，23-29

[10]李媛媛,徐岩,王靖岳. 对MATLAB实现数字信号的频谱分析[J]. 通信技术, 2008, (01) ，69-75

[11]魏克新. MATLAB语言与自动控制系统设计. 北京: 机械工业出版社.2002，1-150 [12]高西全，丁玉美．数字信号处理．第3版．北京：西安电子科技大学出版社，2008,177-189

[13]刘泉，阙大顺．数字信号处理原理．北京：电子工业出版社，2005，88-92

Based on the MATLAB spectrum analysis and signal

denoising

Abstract:Subject is based on MATLAB speech signal processing, the design and implementation of

integrated use of digital signal processing theory knowledge of speech signal with noise adding time domain and frequency domain analysis and filtering. Filter design in digital signal processing in occupies an extremely important status, FIR digital filters and filter IIR filter is the important component. MATLAB signal processing toolbox can effectively design of digital filter. Through theoretical derivation corresponding conclusion as a programming tool, using MATLAB for computer implementation. In the design and implementation of process, USES window function method to design FIR digital filters, with bart leavenworth, chebyshev and bilinear anshi IIR digital filter, and design MATLAB as the auxiliary tool complete design of computing and graphics

第 21 页

drawing. Through the design of the filter selvesor simulation and frequency characteristics were analyzed, using MATLAB signal processing toolbox that can effectively quickly design FIR and IIR digital filters, process is simple and convenient, results of various performance indicators are as specified requirements.

Key words:frequency domain analysis, filter, Wavelet denoising,MATLAB

第 22 页

附录一

Fp1=1200; Fp2=3000; Fs1=1000; Fs2=3200; Ft=8000;

wp1=tan(pi*Fp1/Ft); %带通到低通滤波器的转换 wp2=tan(pi*Fp2/Ft); ws1=tan(pi*Fs1/Ft); ws2=tan(pi*Fs2/Ft); w=wp1*wp2/ws2; bw=wp2-wp1; wp=1;

ws=(wp1*wp2-w.^2)/(bw*w);

[n12,wn12]=buttord(wp,ws,1,50,'s'); %求低通滤波器阶数和截止频率 [b12,a12]=butter(n12,wn12,'s'); %求S域的频率响应参数

[num2,den2]=lp2bp(b12,a12,sqrt(wp1*wp2),bw);%将S域低通参数转为带通的 [num12,den12]=bilinear(num2,den2,0.5);%利用双线性变换实现频率响应S域到Z域的转换

[h,w]=freqz(num12,den12);%根据参数求出频率响应 plot(w*8000*0.5/pi,abs(h)); axis([0 4000 0 1.5]); legend('用butter设计'); grid;

附录二

Fp1=1200; Fp2=3000; Fs1=1000; Fs2=3200; Ft=8000;

wp1=tan(pi*Fp1/Ft); %带通到低通滤波器参数转换 wp2=tan(pi*Fp2/Ft); ws1=tan(pi*Fs1/Ft); ws2=tan(pi*Fs2/Ft); w=wp1*wp2/ws2; bw=wp2-wp1;

第 23 页

wp=1;

ws=(wp*wp2-w.^2)/(bw*w);

[n22,wn22]=buttord(wp,ws,1,50,'s'); %求低通滤波器阶数和截止频率 [b22,a22]=butter(n22,wn22,'s'); %求S域的频率响应的参数

[num2,den2]=lp2bp(b22,a22,sqrt(wp1*wp2),bw); %将S域低通参数转为带通的 [num22,den22]=bilinear(num2,den2,0.5);%利用双线性变换实现S域到Z域的转换 [h,w]=freqz(num22,den22); %根据参数求出频率响应 plot(w*8000*0.5/pi,abs(h)); axis([0 4000 0 1.5]); legend('用butter设计'); grid

附录三

[y,fs,nbits]=wavread ('speech');

n = length (y) ; %求出语音信号的长度 noise=0.21*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 S=fft(s); %傅里叶变换 z12=filter(num12,den12,s);

m12=fft(z12); %求滤波后的信号 subplot(2,2,1); plot(abs(S),'g');

title('滤波前信号的频谱'); grid;

subplot(2,2,2); plot(abs(m12),'r');

title('双线性变换法带通滤波后信号的频谱'); grid;

subplot(2,2,3); plot(s);

title('滤波前信号的波形'); grid;

subplot(2,2,4); plot(z12);

title('双线性变换法带通滤波后的信号波形'); grid;

第 24 页

附录四

[y,fs,nbits]=wavread ('speech');

n = length (y) ; %求出语音信号的长度 noise=0.21*randn(n,2); %随机函数产生噪声 s=y+noise; %语音信号加入噪声 S=fft(s); %傅里叶变换 z22=fftfilt(b22,s);

m22=fft(z22); %求滤波后的信号 subplot(2,2,1); plot(abs(S),'g');

title('滤波前信号的频谱'); grid;

subplot(2,2,2); plot(abs(m22),'r');

title('窗函数法带通滤波后信号的频谱'); grid;

subplot(2,2,3); plot(s);

title('滤波前信号的波形'); grid;

subplot(2,2,4); plot(z22);

title('窗函数法带通滤波后的信号波形'); grid;

第 25 页

本文来源：https://www.bwwdw.com/article/0mxr.html