语音信号滤波去噪——使用FLATTOPWIN设计的FIR滤波器

三尺青锋

发布于 2022-12-26 14:15:46

1.2K0

发布于 2022-12-26 14:15:46

摘要 本课程设计主要内容是设计利用窗口设计法选择FLATTOPWIN窗设计一个FIR滤波器，对一段含噪语音信号进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程，通过课程设计了解FIR滤波器设计的原理和步骤，掌握用MATLAB语言设计滤波器的方法，了解FLATTOPWIN对FIR滤波器的设计及编程方法。首先利用windows自带的录音机录制一段语音信号，加入一单频噪声，对信号进行频谱分析以确定所加噪声频率，设计滤波器进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。由分析结果可知，滤波后的语音信号与原始信号基本一致，即设计的FIR滤波器能够去除信号中所加单频噪声，达到了设计目的。 关键词 滤波去噪；FIR滤波器；FLATTOPWIN窗；MATLAB

引言

本课程设计主要解决在含噪情况下对语音信号的滤波去噪处理，处理时采用的是利用窗口设计法选择FLATTOPWIN窗设计的FIR滤波器[1]。通过对比滤波前后波形图和滤波前后语音信号的对比，可以看出滤波器对有用信号无失真放大具有重大意义。

课程设计目的

熟悉Matlab语言环境，掌握Matlab语言的编程规则，利用Flattopwin窗函数设计法来设计符合要求的FIR滤波器来实现语音信号的滤波去噪。并绘制滤波前后的时域波形和频谱图。根据图形分析判断滤波器设计的正确性。通过本次课程设计熟悉利用 flattopwin窗函数法设计FIR滤波器的过程[2]。增强自己独立解决问题的能力，提高自己的动手能力。加深对理论知识联系实际问题的理解。为以后的工作奠定坚实的基础。

课程设计要求

录制一段语音，绘制波形并观察其频谱特点，加入一个带外单频噪声，设计一个满足指标的滤波器，对该含噪语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析，根据结果和学过的理论得出合理的结论[3]。

课程设计平台

MATLB是美国MathWorks公司出品的商业教学软件，用于算法开发。数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。 MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。而且新版本的MATLAB提供了完整的联机查询、帮助系统，极大的方便了用户的使用。简单的编程环境提供了比较完备的调试系统，程序不必经过编译就可以直接运行，而且能够及时地报告出现的错误及进行出错原因分析。 MATLAB是一个高级的矩阵/阵列语言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编写好一个较大的复杂的应用程序（M文件）后再一起运行。新版本的MATLAB语言是基于最为流行的C++语言基础上的，因此语法特征与C++语言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于非计算机专业的科技人员使用。而且这种语言可移植性好、可拓展性极强，这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因[4]。 MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能。函数中所使用的算法都是科研和工程计算中的最新研究成果，而且经过了各种优化和容错处理。在通常情况下，可以用它来代替底层编程语言，如C和C++ 。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。MATLAB的这些函数集包括从最简单最基本的函数到诸如矩阵，特征向量、快速傅立叶变换的复杂函数。函数所能解决的问题其大致包括矩阵运算和线性方程组的求解、微分方程及偏微分方程的组的求解、符号运算、傅立叶变换和数据的统计分析、工程中的优化问题、稀疏矩阵运算、复数的各种运算、三角函数和其他初等数学运算、多维数组操作以及建模动态仿真等。

基本理论

FIR滤波器

FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，又称为非递归型滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。工作原理：在进入FIR滤波器前，首先要将信号通过A/D器件进行模数转换，把模拟信号转化为数字信号；为了使信号处理能够不发生失真，信号的采样速度必须满足奈奎斯特定理，一般取信号频率上限的4-5倍做为采样频率；一般可用速度较高的逐次逼进式A/D转换器，不论采用乘累加方法还是分布式算法设计FIR滤波器，滤波器输出的数据都是一串序列，要使它能直观地反应出来，还需经过数模转换，因此由FPGA构成的FIR滤波器的输出须外接D/A模块。FPGA有着规整的内部逻辑阵列和丰富的连线资源，特别适合于数字信号处理任务，相对于串行运算为主导的通用DSP芯片来说，其并行性和可扩展性更好，利用FPGA乘累加的快速算法，可以设计出高速的FIR数字滤波器[5]。

窗口设计法

窗口设计法是一种通过截断和计权的方法使无限长非因果序列成为有限长脉冲响应序列的设计方法。通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的技术指标。在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。窗口设计法步骤如下： （1）根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N。窗函数的类型可根据最小阻带衰减AS独立选择[6]。 （2）根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)。 （3）由性能指标确定窗函数W(n)和长度N。 （4）求得实际滤波器的单位脉冲响应h(n)， h(n)即为所设计FIR滤波器系数向量b(n)。

（2.1）

常见的窗函数性能表如下表2.1所示：表2.1 常见窗函数性能表

FLATTOPWIN窗

w=Flattopwin (L) 返回L-点Flattopwin窗口中列向量。Flattopwin窗的滤波器的过渡带宽为19.6π/M，最小阻带衰减108db。

时间波形和幅度谱如下图2.2、图2.3：

图2.2 时间波形

图2.3 幅度谱

设计步骤

设计流程图

根据设计的要求，首先自己录制一段语音信号，修改语音文件格式，对语音信号加入噪声干扰，再利用Flattopwin窗设计合理的FIR滤波器。最后用滤波器对干扰后的语音信号进行滤波去噪。具体设计流程图如下图3.1所示：

录制语音信号

从电脑上录制一段语音信号，并命名为“cf.wav”，修改语音文件的格式，并放在E盘目录下。在MATLAB软件中调用wavread函数可采集到语音信号。代码如下：

[x,fs,bits]=wavread('e:\cf.wav');

%fs是生成该波形文件的采样频率，bits是波形文件没样本的编码位数

得到原始语音信号时域波形图如图3.2

图3.2 原始语音信号时域波形图

然后对语音号进行快速傅里叶变换，得到信号的频谱特性，并将原始音乐信号的波形图与加干扰后的波形图进行比较。部分代码如下：

y=x+0.1*sin(fn*2*pi*t);    %给原始信号加噪声
X=abs(fft(x));             %对原始信号进行傅里叶变换
Y=abs(fft(y));             %对加噪声后信号进行傅里叶变换

加噪声前后的时域与频谱图如图3.3

图3.3 加噪声前后的时域图与频谱图比较

滤波器设计

在本次的课程设计中我所采用的是利用Flattopwin窗函数来设计FIR滤波器。其中主要用了freqz_m.m和ideal_lp.m两个自编函数。设计的滤波器图如图3.4

图3.4 滤波器图形

信号滤波处理

滤波器设计完成后，在MATLAB平台上用函数fftfilt实现滤波。得到的滤波前后语音信号的时域波形图和频谱图对比图如图3.5、3.6

图3.5 滤波前后语音信号的时域波形图和频谱图

图3.6 滤波前后语音信号的比较

结果分析

在MATLAB中，对原始的语音信号加噪音，利用FLATTOPWIN窗口设计FIR滤波器，通过所设计的滤波器对加噪声后的信号进行处理。再通过sound(x,fs,bits)函数对滤波后的语音信号进行回放，可以听到滤波之后的信号和原始信号一样清晰。具体代码如下：

sound(x,fs,bits);        %播放原始语音信号
sound (y_fil,fs,bits);   % 播放滤波后的音乐信号

所得结果证明了用Flattopwin窗设计的FIR滤波器和音乐信号去噪设计是成功的。

出现的问题及解决方法

在本次课程设计中我遇到的问题如下： 1、对Flattopwin窗不是很了解，对用Flattopwin窗函数设计FIR滤波器的设计步骤很生疏。 2、语音文件的格式有问题，不知道如何修改。 3、不知道如何调用函数，对MATLAB的使用不熟悉。 4、对滤波器一些参数以及通带阻带等概念不是很清楚。 5、在采用Flattopwin窗函数设计的FIR滤波器时得不到理想的滤波器，因而信号的恢复不是特别理想。

针对以上问题，相应的解决方案如下：

1、自己上网查阅资料，以及查看之前《数字信号处理》教材和向图书馆借阅资料，掌握利用Flattopwin窗函数设计FIR滤波器的方法和步骤。 2、询问同学，找同学帮忙。 3、利用MATLAB这款软件自带的help功能，以及查看老师所给提供的课设资料和上网看阅资料。 4、请教值班老师以及在网上查阅资料。 5、通过不断设置参数的值，最终达到最理想的值，设计出理想的滤波器，使信号得到理想恢复。

结束语

本次的课程设计，我的题目是《语音信号滤波去噪——使用Flattopwin设计的FIR滤波》。在本次课程设计之前，我对Flattopwin窗函数完全没有了解，因此在看到这个题目时，我是一头雾水。但是通过自己翻阅资料和询问同学，我掌握了用Flattopwin窗函数设计FIR滤波器的方法步骤，了解了窗函数的基本设计流程。经过这三周的课程设计，我学会了很多东西，不仅仅是技术，更多的是自学能力的提高以及处理问题的方法。我们通信工程专业是个实践性很强的专业，而我们在校大部分的学习时间都是花在理论学习上面，实践的机会很少。因而我对很多所学的理论知识如何跟实践联系的概念很模糊，这次的课程设计给了我这个机会，加深了我对理论联系实际的理解，增强了自己独立分析问题和解决问题的能力，开阔了自己的思维。还有让我看到了自己的不足，自己对本专业的相关知识掌握的还很少，还有很多知识都没掌握，还让我认识到解决问题的方法、途径很多，做事要开阔自己的思维，看待问题要从多个角度看。在此我要感谢学校为我们提供这次课程设计的机会，感谢老师对我的悉心指导，也感谢同学对我的帮助。这次的课程设计让我理论联系实际，不仅巩固了我们的理论知识，还提高了我的动手能力，在这次课程设计中我所学到的知识是我的财富，让我终身受益。

代码

语音信号滤波去噪设计源程序清单

% 程序名称：c.m
% 程序功能：利用FLATTOPWIN设计的FIR滤波对语音信号进行滤波去噪
% 程序作者： 
% 最后修改日期：2017年3月8日
[x,fs,bits]=wavread('e:\cf.wav');  % 输了参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。
sound(x,fs,bits);      % 按指定的采样率和每样本编码位数回放
N=length(x);       % 计算信号X的长度
fn=2200;            % 单频噪声频率
t=0:1/fs:(N-1)/fs;         % 计算时间范围，样本数除以采样频率
x=x(:,1)';               % 将双声道转为单声道
y=x+0.1*sin(fn*2*pi*t);     % 加噪声
sound(y,fs,bits);         % 应该可以明显听出有尖锐的单品啸叫声
X=abs(fft(x));  Y=abs(fft(y));    % 对原始信号和加噪声信号进行fft变换，取幅度值
X=X(1:N/2); Y=Y(1:N/2);       % 截取前半部分
deltaf=fs/N;            % 计算频谱的谱线间隔
f=0:deltaf:fs/2-deltaf;     % 计算频谱频率范围
figure(1)
subplot(2,2,1);plot(t,x);xlabel('时间t');
ylabel('幅度');title('原始语音信号');
axis([0,4,-1.5,1.5]);

subplot(2,2,2);plot(f,X);xlabel('频率f');ylabel('幅度谱');
title('原始语音信号幅度谱');
axis([-10,6000,0,700]);
subplot(2,2,3);plot(t,y);xlabel('时间');ylabel('幅度');
title('加干扰后的语音信号');
axis([0,4,-1.5,1.5]);
subplot(2,2,4);plot(f,Y);xlabel('频率 f');ylabel('幅度谱');
title('加干扰后的语音信号幅度谱');
axis([-10,6000,0,700]);
 
fpd=2100;fsd=2150;fsu=2250;fpu=2300;Rp=1;As=30;      % 阻带滤波器设计指标
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为单位的数字频率 
wsd=fsd/fs*2*pi;wsu=fsu/fs*2*pi;
M=ceil(6.2*pi/dw)+1;      %计算flattopwin窗设计滤波器时需要的阶数
n=0:M-1;         % 定义时间范围
w_fla=(Flattopwin(M));    % 产生M阶的flattopwin窗
hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); 
% 调用自编函数计算理想阻带滤波器的脉冲响应
h_bs=w_fla'.*hd_bs;       % 用窗口法计算实际滤波器脉冲响应
[db,mag,pha,grd,w]=freqz_m(h_bs,1);  % 调用自编函数计算滤波器的频率特征
 
figure(2)
subplot(2,2,1);plot(w,db);title('滤波器幅度响应图');
xlabel('w/pi');ylabel('db');
axis([0,0.8,-60,10]);
line([0,2],[-As,-As],'color','y','linestyle','--','LineWidth',2);
line([0,2],[-Rp,-Rp],'color','g','linestyle','--','LineWidth',2);
line([wsd,wsd],[-30,10],'color','r','linestyle','--','LineWidth',2);
line([wsu,wsu],[-30,10],'color','r','linestyle','--','LineWidth',2);
subplot(2,2,2);plot(w,mag);title('滤波器幅度响应图');
xlabel('w/pi');ylabel('幅度mag');
axis([0,1,-0.5,1.5]);
subplot(2,2,3);plot(w,pha);title('滤波器幅度响应图');
xlabel('w/pi');ylabel('相位pha');
axis([0,1,-4,4]);
subplot(2,2,4);stem(n,h_bs);title('滤波器幅度响应图');
xlabel('w/pi');ylabel('db');
axis([0,1500,0,1]);
 
y_fil=fftfilt(h_bs,y);  %用设计好的滤波器对y进行滤波
Y_fil=fft(y_fil);Y_fil=Y_fil(1:N/2);  %计算频谱取前一半
figure(3)
subplot(3,2,1);plot(t,x);xlabel('时间t');ylabel('幅度');


title('原始语音信号');
subplot(3,2,2);plot(f,X);xlabel('频率f');ylabel('幅度谱');
title('原始语音信号幅度谱');
axis([0,4000,0,600]);
subplot(3,2,3);plot(t,y);xlabel('时间t');ylabel('·幅度');
title('加干扰后语音信号');
subplot(3,2,4);plot(f,Y);xlabel('频率f');ylabel('幅度谱');
title('加干扰后语音信号幅度谱');
axis([0,4000,0,600]);
subplot(3,2,5);plot(t,y_fil);xlabel('时间t');ylabel('幅度');
title('滤波后语音信号');
subplot(3,2,6);plot(f,Y_fil);xlabel('频率f');ylabel('幅度谱');
title('滤波后语音信号幅度谱');
axis([0,4000,0,600]);
figure(4)
subplot(2,1,1);plot(f,20*log10(Y));grid on; 
subplot(2,1,2);plot(f,20*log10(Y_fil));grid on; 
sound(x,fs,bits);   % 播放原始信号
sound (y_fil,fs,bits);  %播放滤波后的信号
 
[C,B,A]=dir2fs(h_bs);

函数FREQZ_M.M定义：

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;

函数IDEAL_LP.M定义：

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);

该文为笔者课程设计论文，如有不妥，还请指出

本文参与腾讯云自媒体分享计划，分享自作者个人站点/博客。

原始发表：2017-05-15，如有侵权请联系 cloudcommunity@tencent.com 删除

matlab

c++

编程算法

本文分享自作者个人站点/博客前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

matlab

c++

编程算法

登录后参与评论

0 条评论

热度