语音信号的采集与频谱分析(附代码)

《信号与系统》大作业
语音信号的采集与频谱分析
——基于Matlab的语音信号处理
学生姓名：
学号：
专业班级：电子工程学院卓越班
指导老师：
2015年6月22日
摘要
本设计用苹果手机自带的录音设备采集了原始语音，并导入了电脑转成wav格式，然后用MATLAB和Adobe audition对其进行时域分析。

接着利用傅里叶变换进行了频域分析，绘制频谱图，再录制一段加上歌曲的伴奏的语音与原唱进行了对比分析，得出了我与歌星在频域上的差别。

本设计给信号加了两种噪声并通过观察加噪后的频谱和试听回放效果比较加噪前后的差别，
最后，设计了FIR数字低通滤波器和带通滤波器，分析滤波前后的频谱。

再次试听回放效果，得出结论。

关键词：语音、FFT、频谱图、噪声、滤波器
Abstract
This design is based on the general function of Matlab and Adobe edition to deal with Audio signals. The original signals are collected by iPhone’s built-in recording equipment.
First，I compare the file generated by myself with that of thesame song sang by a famous singer.The emphasis is generally laid on analysing the difference in frequncy domain,but time domain will be included too.
After that,two noise signals are added to the original signal respectively and let them pass a filter to analyse it.In the two process mentioned before,I make comparison between the before and after frequency domain.
Sampling Theorem is the base of my design.It is by sampling we can get discrete signals from the original one and draw the image in time domain.Also,fast fourier transform is employed(FFT)to get the signals in frequency domain.The ayalysis of frequency domain is the highlight of this design.
Through this design,I can deepen my comprehension of principles of audio signals and I have learnt how to deal with it.Through met with much hindrance,I improved my skills finally.
Keywords： audio signal、TTT、noise、filter
1 绪论
1.1课题的研究意义
语音信号处理属于信息科学的一个重要分支，它是研究用数字信号处理技术对语音信号进行处理的一门新兴学科，同时又是综合性的多学科领域和涉及面很广的交叉学科，因此我们进行语言信号处理具有时代的意义。

学会运用MATLAB的信号处理功能，采集语音信号，并对语音信号进行滤波及变换处理，观察其时域和频域特性，加深对信号处理理论的理解，这为今后熟练使用MATLAB进行系统的分析仿真和设计奠定基础。

1.2 课题的目的与要求
利用Matlab和Adobe audition对语音信号进行数字信号处理和分析，要求采集语音信号后，在软件平台进行频谱分析；并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，再进行频谱分析和再播放来感受噪声滤除的效果。

2.理论依据
采样频率、采样位数的概念，采样定理；时域信号的FFT分析；数字滤波器设计原理和方法，各种不同类型滤波器的性能比较。

2.1. 时域信号的FFT分析
FFT即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

它利用“一分为二”的思想不断进行，点数越多，运算量的节约就越大，这就是FFT的优越性。

在MATLAB 的信号处理工具箱中函数FFT 和IFFT 用于快速傅立叶变换和逆
变换。

函数FFT 用于序列快速傅立叶变换，其调用格式为y=fft(x)，其中，x 是序列，y 是序列的FFT ，x 可以为一向量或矩阵，若x 为一向量，y 是x 的FFT 且和x 相同长度；若x 为一矩阵，则y 是对矩阵的每一列向量进行FFT 。

如果x 长度是2的幂次方，函数fft 执行高速基－2FFT 算法，否则fft 执行一种混合基的离散傅立叶变换算法，计算速度较慢。

函数FFT 的另一种调用格式为y=fft(x,N)，式中，x ，y 意义同前，N 为正整数。

函数执行N 点的FFT ，若x 为向量且长度小于N ，则函数将x 补零至长度N ；若向量x 的长度大于N ，则函数截短x 使之长度为N ；若x 为矩阵，按相同方法对x 进行处理。

2.1.2采样定理
采样定理，又称香农采样定律、奈奎斯特采样定律，是信息论，特别是通
讯与信号处理学科中的一个重要基本结论。

采样是将一个信号（即时间或空间上的连续函数）转换成一个数值序列
（即时间或空间上的离散函数），这个数值序列即是一个数字信号。

采样定理是将连续的模拟信号转变为离散的数字信号的方法。

对于一个时域连续的信号f(t) 利用开关函数s(t)进行采样，或者冲激函
数序列进行采样（()T t δ）。

两种取样都具有一定的取样间隔，同时要求f(t)为频带带限信号，否则取样后将发生频谱混合重叠，破坏原信号。

采样定理公式如下：
其中Ts 为采样间隔，设f(t)的频谱区间在
(,)m m ωω-之间，即信号频带带限，需满足的采样条件为：2s m f f ≥，其中fs 为采样频率， 这样便完成了对于连续信号的采样形成离散的数字信号方便处理。

2.1.2 数字滤波器
()()()()
2s s s n t f t f nT t nT Sa ωδ∞=-∞=
-*∑2m ωπ
数字滤波器是一个输入输出都是离散时间信号的特定功能装置。

应用数字
滤波器处理模拟信号（对应模拟频率）时，首先须对输入模拟信号进行限带、
抽样和模数转换。

数字滤波器输入信号的数字频率（2π*f/fs,f为模拟信号的
频率，fs为采样频率，注意区别于模拟频率),按照奈奎斯特取样定理，要使抽
样信号的频谱不产生重叠，应小于折叠频率（ws/2=π)，其频率响应具有以
2π为间隔的周期重复特性，且以折叠频率即ω=π点对称。

为得到模拟信号，数字滤波器处理的输出数字信号须经数模转换。

数字滤波器具有高精度、高可
靠性、可程控改变特性或复用、便于集成等优点。

数字滤波器有低通、高通、带通、带阻和全通等类型。

它可以是时不变的
或时变的、因果的或非因果的、线性的或非线性的。

应用最广的是线性、时不
变数字滤波器，以及FIR滤波器。

从单位脉冲响应的角度，可以把数字滤波器分为：IIR滤波器（无限长单
位冲激响应滤波器）和FIR滤波器（有限长单位冲激响应滤波器）
3.设计思想
3.1语音信号的采集
此次我选择的歌曲是赵薇的《情深深雨蒙蒙》，用苹果手机自带的录音功
能（Voice Memos）录音，和windows自带的录音机程序相比，优点是方便，而且可以忽略电脑运行本身的噪声。

然后将其导入到电脑，转换格式存为wave格式，并对从网上下载原唱的语音用同样的方式转换格式，并保存到正确的目录，以便以后能在matlab中打开。

同时手机自带的录音功能还能完成简单的音频剪辑操作，减掉多余的不需要的部分。

3.2 语音处理工具的选择
① Matlab
Matlab是大家都熟悉的一种科学计算软件，在信号与系统中运用很广泛，
优点是可编程性强，且具有很好的可移植性和可扩展性，而且和专业信号与系
统关联更大。

但缺点是数据都以矩阵形式进行处理，对我来说不够直观。

② Adobe audition
Adobe audition是一款专门处理音频的软件，可提供先进的编辑、混音、
控制和效果处理功能，优点是很方便，比如采样频率可以直接选择，处理效果
很直观，缺点是和信号与系统相关性较小，可编程性差。

这次设计的频谱分析主要是采用了matlab进行分析，但对原始语音频谱的分析以及由采样频率不同带来的变化中我同时采用了matlab和adobe
audition来分析。

3.3噪声函数的选择
①rand函数产生随机噪声，具体语句将在下面的列表中列出
②高斯白噪声Gaussian White Noise
此噪声在《随机信号分析课程》中已学过，此次向原始信号中加入此噪声，加强理解。

3.4 滤波器的选择
设计中选择的是FIR有限长单位冲激响应滤波器，相比于IIR滤波器，它
有如下的优点：
⑴ FIR滤波器，由于冲激响应有限长，可以用快速傅里叶变换，这样运算
速度可以快很多；
⑵ FIR滤波器可以做到严格的线性相位，而若要IIR得到严格的线性相位，则要大大增加滤波器的阶数，
⑶ FIR滤波器主要采用非递归结构，因而从理论上以及时性从实际的有限
精度的运算中，都是稳定的。

有限精度运算误差也较小，IIR滤波器必须采用
递归的结构，极点必须在Z平面单位圆内，才能稳定，这种结构，运算中的四
舍五入处理，有时会引起寄生振荡。

但是IIR也有相同技术指标下阶数较少，存储单元少，运算次数少的有点，此处不再加以赘述。

两种滤波器：
①低通巴特沃斯滤波器
巴特沃斯滤波器的特点为的通频带频率响应最为平滑，阻带缓慢衰减。

用窗函数法设计低通滤波器，具体指标为
fp=3000Hz fc=4000Hz Ap=1dB As=100dB
②高通滤波器
具体指标为：fp=35ooHz fc=4000Hz Ap=1dB As=10dB
3.5 Matlab函数的选择
（1）[y,fs,bits]=wavread(‘filename’)
wavread函数可读入语音信号y是采样获得的离散时域信号，fs是返回的采样
频率， bits是采样位数，filename是音频信号的文件名。

（2）wavwrite(y,fs,bits,’filename’)
wavwrite函数，改变采样频率，可以将22.05KHz或者11.025KHz的信号写回
并观察这样的采样频率会产生什么样的结果。

其中y是已经获得的时域离散信号，fs为写回的频率，bits是写回的位数，filename为要写回的文件名。

（3）sound(y,fs,bits)
sound函数，播放存在的离散音频信号其中的三个参数同wavread中的相同
（4）rnt=rand(a,b)
rand函数，为产生随机噪声，从而在信号中加入噪声。

改随机噪声是在（0,1）内的随机数，均值为0.5其中a、b为产生随机数矩阵的行和列。

（5）fy=fft(y)
y为离散时域信号,fft函数是所有程序分析的核心。

fft函数提供了以快速傅
立叶变换算法的离散傅里叶变换的计算，等于对一个离散信号进行DFT离散傅立叶变换。

（6）y1=awgn(y,q,’measured’)
awgn函数，为向一个现有的离散信号加入高斯白噪声的函数， y为被添加噪声的信号，q为信噪比，一般要大于1，measured是一种控制名称。

（7）[n,fo,ao,w]=firpmord(f,a,dev,fs)
firpmord函数，为设计滤波器阶数的函数，f是滤波器的参数，主要为带边沿的参数，低通滤波器则是通带与阻带截止频率。

a是通带与止带的幅频响应的
幅度，也是一个向量。

dev是通带与止带的震荡幅度大小比例（介于0到1），也是一个向量，fs则是采样频率。

左边的4个参数，n是滤波器阶数，f0是频率向量，a0是幅频特性向量，w是权重。

这四个参数是用来下一步制作模拟滤波器的。

（8）h=firpm(n,fo,ao,w)
firpm函数，利用firpmord所产生的参数制作模拟滤波器，后面四个参数
是firpmord所计算出的。

（9）[ht,w]=freqz(h)
freqz函数，对一个设计好的滤波器求冲激响应，h为前一步firmp设计好的
滤波器，ht为频率响应，用abs取绝对值得到幅频特性，或者用log转化成dB 作图。

w是频率范围（归一化的）
（10） yf=fftfilt(h,y)
fftfilt函数，用设计好的滤波器滤波（数字信号）h为设计好的滤波器，y为待滤的信号。

（11）h=fir1(n,w,window)
fir1函数，用于使用做好的窗函数进行带通滤波器设计，做好窗函数如
hanning窗（MATLAB也有这个函数），根据相应其他参数做出带通滤波器。

window是已经做好的窗，n是滤波器阶数，w是通带范围（一定要在0到1
之间），为一个向量（因为是带通滤波器）。

（13） a=mean(y)
mean函数，用于求一个离散序列的均值，在实验中求采样所得的音频信号的均
值y是一个向量（数组），a是该数组的均值。

（14） a=var(y)
var函数，用于求离散序列的方差，在实验中用于求采样所得音频信号的方差， y为一个数组，a为该数组的方差。

4. 课程设计的具体实现及仿真结果分析
4.1 语音信号的采集、录入、与打开
把我的原始语音信号保存为qingshenshen.wav保存在相应路径，此语音信号的长度小于30s。

通过matlab的函数wavread（）可读入一个.wav格式的音频文件，例如下面的语句
[y,fs,bits]=wavread(‘qingshenshen.wav’);
读入wav文件并放入向量y中，对于双声道的音频文件，y有两行，此处
我们仅对一个声道的音频进行分析。

fs为采样频率，.wav格式默认的采样频率为44.1KHz，下文将比较不同采样频率对语音信号的影响。

回放声音：用sound（）函数，例如：
sounf(x,fs,bits);
此处我将我唱的歌和伴奏合起来保存为qingshenshen.wav
（具体程序见附录一）
运行结果如下
原始语音信号时域频域图
此时可在workspace中看到fs的值value为44100
通过调整放大第三章小图，可以看出我声音的频率范围主要集中在200-450hz 之间属于女低音
原始语音信号频域幅值图
用Adobe audition画出的频域图为，分析可得频域范围与Matlab分析所得很接近，只是Adobe audition的横轴显示的标准不同。

然后我将原唱赵薇的（带伴奏）进行同样的分析，（具体程序见附录一）
原唱语音信号时域频域图
原唱语音信号频域幅值图
通过对比可以发现，赵薇的频域更宽，从几十到1000hz都有较大的分布，算是个女高音吧，但大部分集中在几十到600hz这个区间段，看来我在高音部分与原唱有较大差距，怪不得唱她的歌的时候总有种唱不上去的感觉。

特征值的分析:
通过matlab函数 mean（）和var（）的计算，可以得出原唱的均值与方差为：
m1= 1.0e-003 *0.1372
v1=0.8996
而我的语音的均值与方差为：
m2=1.0e-005*0.5838
v2=0.0046
比对可知，我唱的歌在均值和方差上和原唱差距都很差，说明我对歌曲的还原度不高，还需要加强唱歌技巧的联系。

特别要注意的是，在这次作图过程中，有两个横坐标的量要加以特别的控制，一个是时域分析的横坐标时间，所以加入了对横坐标的控制语句t=(0:n-1)/fs;，使它的单位为时间s；二是频域分析的横坐标控制语句
yol=(0:length(X)-1)'*fs/length(X);如果不加这句很可能出来的不是应有的频率范围。

上面的分析都是默认频率44.1KHz下进行的，可以通过wavwrite函数改变采样频率，例如：
wavwrite(x,220500,’qingshenshen.wav’)
以上语句便使得采样频率变成了22.05khz
然后用sound函数试听，发现音乐失真很大，速度变得很快，原本30s的音频几秒便放完了，同时音调也变得很高，这是我最直观的感受。

同理，如果用500hz的频率去采样，声音就变得很低沉，像风的呼啸与野兽的低喘，可在workspace中观察到fs的值变了
4.2噪声的加入
此处，我对我的原始语音分别加入了两个噪声，一个是随机序列噪声，一个是高斯白噪声。

注意，以下的分析全都是基于44.1khz的采样频率下进行的。

⑴随机序列噪声
程序中用rand函数产生与采样获得的信号y维度相同的噪声信号，加入原信号中，同样进行时域与频域的分析，画出时域图像、幅频特性，相频特性用处不大在此忽略不画。

这里rand的维度通过y的维度得到，可用语句whos查看，
（具体程序见附录二）
运行结果如下：
与上文未加噪声的波形对比，可看出时域波形幅值变化较大，波形更加参差不齐，二频域波形变化不大。

用sound回放，可明显感觉到有杂音，就像老式的电视机放出的那种刺啦刺啦的声音。

然后再计算加噪后的均值与方差：
mean=1.0e-004*0.1374
var=0.048
与加噪声前相比，均值变增大了135%，而方差只变化了4.3%，主要因为加入了大小同信号大小幅度差不多的噪声，并且全都在大于0的区域。

此处分析时我仅仅分析了我翻唱的音频信号，对翻唱的信号加噪声分析，对于原唱的分析类似，只要改变文件名即可，因此不再赘述。

同样在画图时有坐标控制，此后的程序也有，不再加以称述。

⑵高斯白噪声
用函数awgn（）产生高斯白噪声，这里取信噪比为10dB,防止噪声太大而淹没原始信号。

（具体程序如下见附录二）
运行结果如下：
从图片上可以看出，通过对比可以明显看出加完噪声后时域幅值变大，频域变
化不大。

加入高斯白噪声后试听的效果同加入随机噪声的结果都有刺啦刺啦的声音但相
对而言高斯白噪声的噪声更大。

同样对于歌星原唱的音频信号不再重复操作。

均值与方差分别为:
mean2=1.0e-001*1.2676
var2=0.0460
与未加噪声前相比，均值变化很大，而方差仍基本没有变化。

4.3滤波器的设计及滤波结果
⑴低通滤波器
具体的指标为：通带截止频率为3000Hz ，阻带截止频为4000Hz，通带最大衰减为1dB，阻带最大衰减为100dB. 通带阻带的幅频相应幅度一个是1一个是0，震荡比例在通带是0.05，阻带为0.01。

并用frez计算频率响应。

（具体程序见附录三）
设计出滤波器后，画出幅频响应和相频响应图为：
通过滤波器前后的噪声的时域和频域图为：
对比可看出，低频部分，尤其是频率小于3000hz的部分的频谱基本没有发生变化，而在蓝色的途中看以明显看出绿色图中的高频部分被滤除了。

试听感觉，声音相比未滤波之前的更加沉闷一些，因为高频成分被滤去的缘故。

⑵高通滤波器
此处采用凯塞窗函数构成高通滤波器，设计指标为：阻带截止频率为
3500Hz，通带截止频率为4000Hz,阻带衰减为1dB,通带衰减为100dB。

（具体程序见附录三）滤波器的频率响应图为：
加噪信号通过滤波器前后频谱的变化如下图：
从视觉上可以很明显地看出，通过滤波器后，低频信号都被滤除了，只剩蓝色的高频分量，同时滤波前后的波形变化也很大，滤波后的波形更为规整。

试听后觉得声音很尖，是因为低频部分的分量被滤除的原因。

5.设计中遇到的问题及总结
5.1问题
遇到的问题首先与矩阵的维度有关，Matlab是以矩阵形式存储和处理数据的，而我经常忽略矩阵的维度而直接对两个矩阵进行加减运算，无疑会出现矩阵维度不匹配的问题，一开始我并不知道数据的维度究竟是多少，也不知道如何解决，但后来我学会了用whos语句查看变量的维度，之后又学会了怎么取出矩阵的某一行某一列，问题便得到了解决。

然后遇到的问题是如何改变采样频率，一开始我会直接在程序开头定义采样频率为22050Hz或者5500Hz，但之后发现这样做出来的试听效果是一样的，后来查阅了才知道wave格式默认的采样频率为44.1kHz，要用wavwrite函数改变信号的采样频率，问题迎刃而解。

接着还遇到了滤波器怎么设计的问题，这方面老师上课虽然有提到过，但是并没有深入地讲，但我记得老师说过滤波器是一个很专业的东西，有人为此设计了一辈子。

然后我通过查阅资料了解了IIR滤波器和FIR滤波器，并知道了通过Matlab如何去模拟滤波器。

虽然如此我觉得我设计的滤波器还只是最最简单的，从滤波效果来看，还有许多问题，还需加强对滤波器的学习。

5.2设计的
5.2感悟与总结
通过这次设计，又加深了我对傅里叶变换以及采样定理的理解，深感信号与系统真是博大精深，我们平时所学都只是冰山一角，只是皮毛，更多的内容在等着我们去探索。

一开始想直接用Adobe audition进行滤波，但是说来惭愧，连教程都看不懂，想要将歌曲的伴奏滤都做的很差。

想想也是，哪有那么轻易的事，怎么可能用简答的几步就完成人家苦苦研究的课题。

后来用Matlab编程，感觉自己的编程能力还是很欠缺，遇到了很多问题，一些细节一些小的指令总是被我轻易的忽略。

但我同时也再次感受到了Matlab 的神奇，它有那么多功能，和信号与系统一样，也许我原来知道的只是其中很小的一部分。

总的来说，在以后的学习中，要一手抓理论知识，即掌握老师上课讲的内容，同时也要加强编程能力的锻炼。

对于这次试验所做的内容，其实只是对语音信号的一个粗略的分析，通过结果也可以看出，滤波器设计得并不理想。

由于本人能力和时间有限，所以未进行更深层次的研究，以后一定会加强。

参考资料：
【1】郭宝龙，闫允一，朱娟娟，吴宪祥．工程信号与系统［M］．北京：高等教育出版社．
【2】王祯飞．基于 MATLAB 的语音信号分析和处理［C］．福建师范大学协和学院．
【3】丁玉美，高西全，数字信号处理. 3 版[M]. 西安：西安电子科技大学出版社，2008.
附录
附录一：语音信号的打开与频域分析
%我的语音（原始语音）
fs=44100; %44.1KHz为默认采样频率，可改变[x,fs,bits]=wavread('qingshenshen.wav'); %读入原始语音：情深深雨蒙蒙sound(x,fs,bits);
X=fft(x);
aX=abs(X);
figure; %产生一个新图形界面
n=length(x); %取y数据的长度；
t=(0:n-1)/fs; %求时间单位s的运算公式；subplot(311);plot(t,x);
xlabel('s');
title('时域波形图') %绘制以时间单位s的时域波形subplot(312);plot(X);
title('原始信号频谱'); %画出原始信号频谱图
subplot(313);
yol=(0:length(X)-1)'*fs/length(X); %控制横坐标范围
plot(yol(1:length(yol)/2),aX(1:length(yol)/2)); %画出原始信号频谱的幅值
axis([0 1500 0 3000]);
xlabel('Hz');
title('原始信号幅值');
m2=mean(x);
v2=var(x);
%原唱语音
fs=44100; %44.1KHz为默认采样频率，可改变[x,fs,bits]=wavread('yuanchang.wav'); %读入原始语音：yuanchang.wav sound(x,fs,bits);
X=fft(x);
aX=abs(X);
figure; %产生一个新图形界面
n=length(x); %取y数据的长度；
t=(0:n-1)/fs; %求时间单位s的运算公式；subplot(311);plot(t,x,’g’);
xlabel('s');
title('时域波形图') %绘制以时间单位s的时域波形subplot(312);plot(X);
title('原始信号频谱'); %画出原始信号频谱图
subplot(313);
yol=(0:length(X)-1)'*fs/length(X); %控制横坐标范围
plot(yol(1:length(yol)/2),aX(1:length(yol)/2)); %画出原始信号频谱的幅值
axis([0 1500 0 3000]);
xlabel('Hz');
title('原始信号幅值');
m1=mean(x);
v1=var(x);
附录二：像原始语音中加入噪声的Matlab程序
⑴加入随机噪声
%在原始语音中加入rand产生的随机噪声
[y,fs,bits]=wavread('qingshenshen.wav'); %打开原始语音（qingshenshen.wav)
wavwrite(y,44100,'qingshenshen.wav');
N=length(y); %求出语音信号的长度
Noise=0.01*randn(N,2); %随机函数产生噪声，列数为两列，与原始信号维度相同
Si=y+Noise; %语音信号加入噪声
sound(Si,fs,bits);
subplot(2,1,1);
n=length(y); %取y数据的长度；
t=(0:n-1)/fs; %控制时间轴的坐标s
plot(t,Si); %画出加噪后的时域波形
xlabel('s');
title('加噪语音信号的时域波形');
S=abs(fft(Si)); %傅里叶变换
subplot(4,1,4);
yol=(0:length(S)-1)'*fs/length(S); %控制频域轴坐标
plot(yol(1:length(yol)/2),S(1:length(yol)/2)); %画出加噪后的频谱
xlabel('Hz');
axis([100 1000 0 3000]); %坐标范围控制
title('加噪语音信号的频域波形');
⑵加入高斯白噪声
%加入高斯白噪声
[y,fs,bits]=wavread('qingshenshen.wav');
y=y*3; %放大信号
y2=awgn(y,10,'measured'); %加入高斯白噪声
sound(y2,fs,bits) %试听加高斯白噪声后的音频
fy2=fft(y2); %做傅里叶变换
y2l=(0:length(fy2)-1)'*fs/length(fy2);
subplot(211),plot(y2); %画出加入噪声后的时域图像
title('加高斯白噪声后的歌声信号时域图像');
axis([0,1200000,-1,1]); %坐标范围控制
subplot(212),plot(y2l(1:length(y2l)/2),abs(fy2(1:length(y2l)/2))); %画出加入噪声后的频域图像
axis([0 2000 0 4000]);
title('加噪声后的歌声信号幅度谱');
附录三.滤波器的设计与滤波的Matlab程序
⑴低通滤波器的设计与滤波
%设计低通滤波器
fs=44100;
f=[3000,4000]; %通带、阻带截止频率
n=[1,0];
rip=[0.05,0.01];
[N,fo,mo,w]=firpmord(f,n,rip,fs); %设计滤波器阶数
N=N+1;
hn=firpm(N,fo,mo,w); %制作滤波器
[h,W]=freqz(hn); %求频率响应
h1=abs(h);
h1=20*log10(h1); %转换单位为dB
subplot(211),plot(W*fs/(2*pi),h1); %画出滤波器的幅频响应
xlabel('频率(Hz)');
ylabel('幅度(dB)');
title('低通滤波器幅频特性');
axis([0,8000,-100,10]);
subplot(212),plot(W*fs/(2*pi),angle(h));%画出滤波器的想频响应
xlabel('频率(Hz)');
title('低通滤波器相频特性');
%用设计的低通滤波器滤波
[y,fs,bits]=wavread('qingshenshen.wav');
n=length(y); %求出语音信号的长度
noise=0.01*randn(n,2); %随机函数产生噪声
s=y+noise; %语音信号加入噪声
S=fft(s); %傅里叶变换
z21=fftfilt(hn,s); %让加噪信号通过滤波器hn
sound(z21,fs,bits);
m21=fft(z21); %求滤波后的信号
subplot(2,2,1);
y20=(0:length(S)-1)'*fs/length(S);
plot(y20(1:length(y20)/2),abs(S(1:length(y20)/2)),'g'); %画出未滤波器按信号的频谱
title('滤波前信号的频谱');
subplot(2,2,2);
y2l=(0:length(m21)-1)'*fs/length(m21);
plot(y2l(1:length(y2l)/2),abs(m21(1:length(y2l)/2))); %画出滤波后信号的频谱title('滤波后信号的频谱');
subplot(2,2,3);
plot(s);
title('滤波前信号的波形');
subplot(2,2,4);
plot(z21);title('滤波后的信号波形');
⑵高通滤波器的设计与滤波
%设计高通滤波器
rp=1;
rs=100;
p=1-10.^(-rp/20); %通带阻带波纹
s=10.^(-rs/20);
wp=0.9;
ws=0.7;
fpts=[ws wp];
mag=[0 1];
dev=[p s];
[n23,wn23,beta,ftype]=kaiserord(fpts,mag,dev);
b23=fir1(n23,wn23,'high',Kaiser(n23+1,beta)); %由fir1设计滤波器
[h,w]=freqz(b23,1); %得到频率响应
plot(w*12000*0.5/pi,abs(h)); %画出频率响应图
title('FIR高通滤波器');
axis([3000 6000 0 1.2]);
%用设计的高通滤波器绿除噪声
[y,fs,bits]=wavread('qingshenshen.wav');
n=length(y); %求出语音信号的长度
noise=0.01*randn(n,2); %随机函数产生噪声
s=y+noise; %语音信号加入噪声
S=fft(s); %傅里叶变换
z21=fftfilt(b23,s); %让加噪信号通过滤波器hn
sound(z21,fs,bits);
m21=fft(z21); %求滤波后的信号
subplot(2,2,1);
y20=(0:length(S)-1)'*fs/length(S);
plot(y20(1:length(y20)/2),abs(S(1:length(y20)/2)),'g'); %画出未滤波器按信号的频谱
title('滤波前信号的频谱');
subplot(2,2,2);
y2l=(0:length(m21)-1)'*fs/length(m21);
plot(y2l(1:length(y2l)/2),abs(m21(1:length(y2l)/2))); %画出滤波后信号的频谱title('滤波后信号的频谱');
subplot(2,2,3);
plot(s);
title('滤波前信号的波形');
subplot(2,2,4);
plot(z21);title('滤波后的信号波形');。