音乐信号滤波去噪

音乐信号滤波去噪

使用三角窗设计的FIR滤波器

学生姓名：林应盛指导教师：胡双红

摘要本次课程设计是用采集一段音乐信号并对该音乐信号进行滤波去噪处理。课程设计平台为MATLAB7.0。设计步骤为：首先采集一段音乐信号并观察其频谱，然后设计一个三角窗FIR滤波器，最后对该信号进行滤波。信号在进行滤波处理后，观察并记录滤波前后波形和频谱的变化，能够听到滤波后的音乐信号和滤波前相比明显的变得清晰，基本达到了设计目的。

关键词课程设计；滤波去噪；FIR滤波器；三角窗；MATLAB7.0

1 引言

在数字信号处理中，滤波占有极其重要的作用，数字滤波器是谱分析、雷达信号处理、通信信号处理应用中的基本处理算法。目前常用的滤波器设计方法普遍采用Matlab 仿真，DSP实现。音乐信号的处理与滤波的设计主要是用Matlab作为工具平台，设计中涉及到音乐的录制、播放、存储和读取，音乐信号的抽样、频谱分析，滤波器的设计及音乐信号的滤波，通过数字信号处理课程的理论知识的综合运用。从实践上初步实现对数字信号的处理。

1.1课程设计的目的

设计一个FIR滤波器，可以有多种方法，窗函数法是设计FIR数字滤波器的最简单的方法。它在设计FIR数字滤波器中有很重要的作用，正确地选择窗函数可以提高设计数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶次。常用的窗函数有以下几种：矩形窗(Rectangular window)、三角窗(Triangular window)、汉宁

窗(Hanning window)、汉明窗(Hamming window)、布拉克曼窗(Blackman window)、切比雪夫窗(Chebyshev window)、巴特里特窗(Bartlett window)及凯塞窗(Kaiser window)。

在本次课程设计的目的是如何设计一个三角窗FIR滤波器，从而达到对音乐信号滤波的效果。

1.2课程设计的要求

（1）滤波器指标必须符合工程实际。

（2）设计完后应检查其频率响应曲线是否满足指标。

（3）处理结果和分析结论应该一致，而且应符合理论。

（4）独立完成课程设计并按要求编写课程设计报告书。

1.3设计平台MATLAB

MATLAB是由美国Math Works公司20世纪80年代中期推出的数学软件。MATLAB是“Matric Laboratory”的缩写，意及“矩阵实验室”，优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。Matlab已经发展成为多学科、多种工作平台的功能强大的大型软件。在欧美的高校和研究机构中，MATLAB是一种非常流行的计算机语言，许多重要的学术刊物上发表的论文均是用MATLAB来分析计算以及绘制出各种图形[1]。

MATLAB是一完整的并可扩展的计算机环境，是一种进行科学和工程计算的交互式程序语言。它的基本数据单元是不需要指定维数的矩阵，它可直接用于表达数学的算式和技术概念，而普通的高级语言只能对一个个具体的数据单元进行操作。因此，解决同样的数值计算问题，使用MATLAB要比使用Basic、Fortran和C语言等提高效率许多倍。许多人赞誉它为万能的数学“演算纸”。MATLAB采用开放式的环境，你可以读到它的算法，并能改变当前的函数或增添你自己编写的函数。

MATLAB包含的内容非常丰富，功能强大，可以概括为以下几个方面：

1)可以在多种操作系统下运行，如DOS、Windows 95/98/2000/2000/NT、Compaq Alpha、LinuxSun Solaris等。

2)有超过500种的数学、统计、科学及工程方面的函数，使用简单快捷，并且有很强的用户自定义函数的能力。

3)有强大的图形绘制和可视化功能，可以进行视觉数据处理和分析，进行图形、图像的显示及编辑，能够绘制二维、三维图形，使用户可以制作高质量的图形，从而写出图文并茂的文章。

4)有从外部文件及外部硬件设备读入数据的能力。

5)有丰富的工具箱〔toolbox〕。各个领域的专家学者将众多学科领域中常用的算法编写为一个个子程序，即m文件，这些m文件包含在一个个工具箱中。其工具箱可以分为两大类，即功能性工具箱和科学性工具箱。功能性工具箱主要用来扩充MATLAB的符号计算、图形可视化、建模仿真、文字处理等功能以及与硬件实时交互的功能。学科性工具箱是按学科领域来分类的，如信号处理、控制、通信、神经网络图像处理、系统辨识、鲁棒控制、模糊逻辑、小波等工具箱[2]。

MATLAB中的信号处理工具箱内容丰富，使用简便。在数字信号处理中常用的算法，如FFT，卷积，相关，滤波器设计，参数模型等，几乎都只用一条语句即可调用。数字信号处理常用的函数有波形的产生、滤波器的分析和设计、傅里叶变换、Z变换等。

2 设计原理

2.1 FIR滤波器

数字滤波器根据其冲激响应函数的时域特性，可分为2种，即无限长冲激响（IIR）滤波器和有限长冲激响应（FIR）滤波器。FIR和IIR的滤波原理都是进行卷积，就是对数入信号进行某种计算。FIR用处就在于对数字信号进行必要的处理，得到所需的输出信号。

FIR系统有自己突出的优点：系统总是稳定的；易实现线性相位；允许设计多通带（或多阻带）滤波器，后两项都是IIR系统不易实现的。FIR数字滤波器的设计方法有多种，如窗函数设计法、频率采样法和Chebyshev逼近法等。随着Matlab软件尤其是

Matlab 的信号处理工作箱的不断完善，不仅数字滤波器的计算机辅助设计有了可能，而且还可以使设计达到最优化。

FIR 数字滤波器设计的基本步骤如下： （1）确定技术指标

在设计一个滤波器之前，必须首先根据工程实际的需要确定滤波器的技术指标。在很多实际应用中，数字滤波器常被用来实现选频操作。因此，指标的形式一般在频域中给出幅度和相位响应。幅度指标主要以2种方式给出。第一种是绝对指标。他提供对幅度响应函数的要求，一般应用于FIR 滤波器的设计。第二种指标是相对指标。他以分贝值的形式给出要求。本文中滤波器的设计就以线性相位FIR 滤波器的设计为例。 （2）逼近

确定了技术指标后，就可以建立一个目标的数字滤波器模型（通常采用理想的数字滤波器模型）。之后，利用数字滤波器的设计方法（窗函数法、频率采样法等），设计出一个实际滤波器模型来逼近给定的目标。 （3）性能分析和计算机仿真

上两步的结果是得到以差分或系统函数或冲激响应描述的滤波器。根据这个描述就可以分析其频率特性和相位特性，以验证设计结果是否满足指标要求；或者利用计算机仿真实现设计的滤波器，再分析滤波结果来判断。

2.2窗口设计法

FIR 滤波器的设计方法有许多种，如窗函数设计法、频率采样设计法和最优化设计法等。窗函数设计法的基本原理是用一定宽度窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列，主要设计步骤为：

(1) 通过傅里叶逆变换获得理想滤波器的单位脉冲响应hd(n)。

())

()]

(sin[21

a n a n w dw e n h c ax

ax

jwn d --=

=

?

-ππ

（2-1）

(2) 由性能指标确定窗函数W(n)和窗口长度N 。

(3) 求得实际滤波器的单位脉冲响应h(n)， h(n)即为所设计FIR 滤波器系数向量b(n)。 )()()(n W n h n h d ?= （2-2）

(4) 检验滤波器性能。

2.3三角窗

三角窗是最简单的频谱函数w(e jwt )为非负的一种窗函数。三角窗函数的时域形式可以表示为： 当n 为奇数时

W(k)=???

????≤≤+++-+≤≤+n k n n k n n k n k

21,1)1(22

11,12 （2-3）

当n 为偶数时

W(k)=???

????≤≤+-≤≤-n k n n k n n k n k 2,)1(22

1,1

2 （2-4）

三角窗函数的主瓣宽度为N

π

8，比矩形窗函数的主瓣宽度增加了一倍，但是它的旁宽度却小得多。

Triang 函数：生成三角窗

调用方式 w=triang(n);输入参数是n 是窗函数的长度；输出参数w 是窗函数的值组成的n 阶向量。

三角窗也是两个矩形窗的卷积。三角窗函数的首尾两个数值通常是不为零的。当n 是偶数时，三角窗的傅里叶变换总是非负数。

三角窗与矩形窗比较,主瓣宽约等于矩形窗的两倍,但旁瓣小,而且无负旁瓣 。 三角窗：w=triang(M)，产生一个M 点的汉明窗函数。

表2-1 常见窗函数性能表

3设计步骤

3.1流程

采集音乐信号后对原始音乐信号进行加噪声，然后设计出合适的滤波器，最后用设计好的滤波器对音乐信号进行滤波去噪，观察并记录结果。设计流程如下图所示

图3-1设计流程图

3.2音乐信号的采集

录制一段格式.wav的音乐信号，时间为1-2S。然后在Matlab软件平台下，利用函数wavread对音乐信号进行采样，记住采样频率和采样点数。采集完成后在信号中加入一个单频噪声，试听效果。

[x,fs,bits]=wavread('D:\我的文档\桌面\gg.wav'); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。

sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放

3.3对音乐信号加噪声

计算出音乐信号的长度，然后设定单频噪声频率，并计算出时间范围，样本数除以采样频率，再回放音乐应该可以明显听出有尖锐的单频啸叫声。

N=length(x); % 计算信号x的长度

fn=1800; % 单频噪声频率，此参数可改

t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率

x=x'; y=x+0.05*sin(fn*2*pi*t);

sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声

3.4 分析加噪声前后音乐信号的波形及频谱变化

给音乐信号加入单频噪声后，对原始信号和加噪声信号进行fft变换，去幅度谱。X=abs(fft(x)); Y=abs(fft(y)); % 对原始信号和加噪信号进行fft变换,取幅度谱

X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分

deltaf=fs/2/N; % 计算频谱的谱线间隔

f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围

图3-4 音乐信号加噪声前后时域图和频谱图

3.5 滤波器设计

对音乐信号加噪声后，接下来的工作是设计一个三角窗滤波器，

>> fpd=1700;fsd=1750;fsu=1850;fpu=1900;Rp=1;As=24; % 带阻滤波器设计指标

>> fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min((fsd-fpd),(fpu-fsu)); % 计算上下边带中心频率，和频率间隔

>> wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率

>> M=ceil(6.1*pi/dw)+1 % 计算三角窗设计该滤波器时需要的阶数

M =

1347

>> n=[0:1:M-1]; % 定义时间范围

>> w_tri=(triang(M))'; 产生M阶的三角窗

>> hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应

>> h= hd_bs.* w_tri; % 用窗口法计算实际滤波器脉冲响应

>> [db,mag,pha,grd,w]=freqz_m(h,[1]); % 调用自编函数计算滤波器的频率特性

图3-5滤波器幅度和相位响应图

3.6用滤波器对音乐信号进行处理

设计好滤波器后，我们将要对音乐信号进行滤波，仔细对比滤波前和滤波后的音乐信号。

y_fil=filter(h,1,y); % 用设计好的滤波器对y进行滤波

Y_fil= fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半

图3-6滤波前和滤波后音乐信号的波形及频谱的对比

3.7结果分析

sound(x,fs,bits);

回放原始音乐信号

sound(y,fs,bits);

回放滤波后的音乐信号

通过观察上图，滤波后的音乐信号发生了衰减，说明滤波器起到了滤波作用，同时通过幅频对比，可以看出滤波器滤掉了一部分频率范围内的信号。分别听原始音乐和滤波后的音乐信号，发现滤波后的音乐信号单频噪声消失，同时原始信号强度稍有减弱，基本达到了滤波的效果。

4出现的问题及解决方法

4.1出现的问题

（1）音乐无法播放，经过同学帮助，原来是要把音乐信号改为单声道信号，把音乐信号改为左声道后能放出音乐。

（2）原始信号和加噪声后的波形和频谱图图显示不出来。

（3）滤波前后单频噪声无法消失。

4.2 解决方法

（1）把音乐信号改为单声道信号后能放出音乐。

（2）改变坐标范围，能显示波形和频谱图。

（3）设计滤波器的时候多写了一行错误的代码，去除后再运行，发现单频噪声消失。

5结束语

这次课程设计基本上是对上学期《数字信号处理》的复习，通过这次课程设计后心得颇多，感受也不少。最令我感触最深的是，我们以前只是按着书本学习，也就是重理论轻实践。通过这次课程设计我是知道自己动手能力差了，理论固然重要，但理论最终还是要应用到实践的。我想今后的学习一定得重视实践，只有理论与实践结合才能学到东西。同时我也还有个比较大的感受就是，一个人的力量是有限的，当遇到问题的时候，仅仅只凭自己单个的力量往往无法将问题解决，得靠他人的帮助，也就是得有团队精神。还有就是学海无涯，尽管自己知道一些，但别人比自己动更多，要虚心请教。

在这我要真诚的感谢胡双红老师，还有那些帮助过我的同学，谢谢你们，你们让我学到的不仅仅是一点知识，还有生活的道理。在未来求学的道路上我将更努力去面对遇到的困难。我下定决心今后一定要把知识学扎实了，这一才能在面对问题是更从容。

参考文献

[1] 程佩青．数字信号处理教程［M］．北京：清华大学出版社，2002

[2] 薛年喜主编．MATLAB 在数字信号处理中的应用[M]．北京：清华大学出版社，2002

[3] 维纳·K·恩格尔，约翰·G·普罗克斯.数字信号处理[M]. 西安：西安交通大学出版社，2002.6

[4] 孙祥主编.MATLAB7.0基础教程［M］．北京：清华大学出版社，2005

[5] 谢德芳.数字信号处理[M].北京：科学出版社，2005

附件1：MATLAB源程序

[x,fs,bits]=wavread('D:\我的文档\桌面\gg.wav');

>> sound(x,fs,bits);

>> N=length(x);

>> fn=1800;

>> t=0:1/fs:(N-1)/fs;

>> x=x'; y=x+0.05*sin(fn*2*pi*t);

>> sound(y,fs,bits);

>> X=abs(fft(x)); Y=abs(fft(y));

>> X=X(1:N/2); Y=Y(1:N/2);

>> deltaf=fs/N;

>> f=0:deltaf:fs/2-deltaf;

>> subplot(2,2,1);plot(t,x);title('原始音乐信号')

axis([0 6 -1 1]);xlabel('时间（单位：s）');ylabel('幅度');grid

>> subplot(2,2,2); plot(f,X);title('原始音乐信号信号幅度谱');

axis([0 2000 0 4000]);xlabel('频率（单位：s）');ylabel('幅度谱');grid

>> subplot(2,2,3);plot(t,y);title('加入单频干扰后的音乐信号')

axis([0 20 -2 2]);xlabel('时间（单位：s）');ylabel('幅度');grid

>> subplot(2,2,4); plot(f,Y);title('加入干扰后的音乐信号信号幅度谱');

axis([0 2000 0 4000]);xlabel('频率（单位：s）');ylabel('幅度谱');grid

>> fpd=1700;fsd=1750;fsu=1850;fpu=1900;Rp=1;As=24; % 带阻滤波器设计指标

>> fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min((fsd-fpd),(fpu-fsu)); % 计算上下边带中心频率，和频率间隔

>> wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; % 将Hz为单位的模拟频率换算为rad为单位的数字频率

>> M=ceil(6.1*pi/dw)+1 % 计算三角窗设计该滤波器时需要的阶数

M =

1347

>> n=[0:1:M-1]; % 定义时间范围

>> w_tri=(triang(M))'; 产生M阶的三角窗

>> hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应

>> h= hd_bs.* w_tri; % 用窗口法计算实际滤波器脉冲响应

>> [db,mag,pha,grd,w]=freqz_m(h,[1]); % 调用自编函数计算滤波器的频率特性

>> subplot(2,2,1);plot(w/pi,db);title('滤波器幅度响应图')

axis([0 0.5 -30 1]);xlabel('w/pi');ylabel('dB');grid

>> subplot(2,2,2);plot(w/pi,mag);title('滤波器幅度响应图')

axis([0 0.5 0 1.15]);xlabel('w/pi');ylabel('幅度mag');grid

>> subplot(2,2,3);plot(w/pi,pha);title('滤波器相位响应图')

axis([0 1 -3 3]);xlabel('w/pi');ylabel('相位pha');grid

>> subplot(2,2,4);stem(n,hd_bs);title('滤波器脉冲响应图');

axis([0 1500 -0.1 1]);xlabel('n');ylabel('h(n)');grid

y_fil=filter(h,1,y); % 用设计好的滤波器对y进行滤波

Y_fil= fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半

subplot(3,2,1);plot(t,x);axis([0 20 -1 1])

xlabel('时间（单位：s）');ylabel('幅度');title('原始音乐信号x');

subplot(3,2,2);plot(f,X);axis([0 5000 0 2500])

xlabel('频率（单位：Hz）');ylabel('幅度谱');title('原始音乐信号幅度谱图X');

subplot(3,2,3);plot(t,y); xlabel('时间（单位：s）');ylabel('幅度');title('加入单频噪声音乐信号y');

axis([0 20 -1 1])

subplot(3,2,4);plot(f,Y); axis([0 5000 0 2500])

xlabel('频率（单位：Hz）');ylabel('幅度谱');title('加入单频噪声信号幅度谱图Y');

subplot(3,2,5);plot(t,y_fil); axis([0 20 -1 1])

xlabel('时间（单位：s）');ylabel('幅度');title('滤波后音乐信号y_fil');

subplot(3,2,6);plot(f,Y_fil); axis([0 5000 0 2500])

xlabel('频率（单位：Hz）');ylabel('幅度谱');title('滤波后音乐信号Y_fil')

附件2：自编函数源程序

function [db,mag,pha,grd,w] = freqz_m(b,a);

% Modified version of freqz subroutine

% ------------------------------------

% [db,mag,pha,grd,w] = freqz_m(b,a);

% db = Relative magnitude in dB computed over 0 to pi radians % mag = absolute magnitude computed over 0 to pi radians

% pha = Phase response in radians over 0 to pi radians

% grd = Group delay over 0 to pi radians

% w = 501 frequency samples between 0 to pi radians

% b = numerator polynomial of H(z) (for FIR: b=h)

% a = denominator polynomial of H(z) (for FIR: a=[1])

%

[H,w] = freqz(b,a,1000,'whole');

H = (H(1:1:501))'; w = (w(1:1:501))';

mag = abs(H);

db = 20*log10((mag+eps)/max(mag));

pha = angle(H);

% pha = unwrap(angle(H));

grd = grpdelay(b,a,w);

% grd = diff(pha);

% grd = [grd(1) grd];

% grd = [0 grd(1:1:500); grd; grd(2:1:501) 0];

% grd = median(grd)*500/pi;

function hd=ideal_lp(wc,M);

%Ideal Lowpass filter computation

%------------------------------------

%[hd]=ideal_lp(wc,M)

% hd=ideal impulse response between 0 to M-1

% wc=cutoff frequency in radians % M=length of the ideal filter alpha=(M-1)/2;

n=[0:1:(M-1)];

m=n-alpha+eps;

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