音乐信号滤波去噪-基于Gaussian的频率采样型FIR滤波器设计(2)

王杰《音乐信号滤波去噪--使用GAUSS窗设计的频率采样型FIR滤波器》第6页共23页

2.3 GAUSSIAN窗

高斯窗是一种指数窗。主瓣较宽，故而频率分辨力低；无负的旁瓣，第一旁瓣衰减达一55dB。常被用来截短一些非周期信号，如指数衰减信号等。对于随时间按指数衰减的函数，可采用指数窗来提高信噪比。 Gausswin的时域表达式可以表示为：

w[k?1]?e1k?N22?(?)2N/2(2-4)

其中k与?的取值范围必须满足0?k?N,??2[2]

2.4滤波器结构

本次课程设计的滤波器采用的是频率采样型结构。 (1) 频率采样型结构的导入

若FIR DF的冲激响应为有限长(N点)序列h(n)，则有：

图2-1 关于h(n)的推导

所以，对h(n)可以利用DFT得到H(k),再利用内插公式：

1N?1H(k)H(z)?(1?z)?(2-5) ?k?1Nk?01?WNz?N来表示系统函数。 (2) 频率采样型滤波器结构

由式2-5得到FIR滤波器的另外一种结构：频率采样型结构。它是由两部分级联而成。

1N?1H(z)?Hc(z)?Hk(z) (2-6)

Nk?0王杰《音乐信号滤波去噪--使用GAUSS窗设计的频率采样型FIR滤波器》第7页共23页

?NH(z)?(1?z)， c其中：级联中的第一部分为梳状滤波器

第二部分由N个谐振器组成的谐振贵。

Hk(z)?H(k)(2-7) ?k?11?WNz(3) 频率采样型结构流图如图2-2所示

图2-2 频率采样型结构流程图

(4) 频率采样型结构特点

它的系数H(k)直接就是在wk?2?k处的频率响应。因此，控制滤波器的频率响N应是很直接的。但所有的相乘系数及H(k)都是复数，应将它们先化成二阶的实数，这样乘起来比较复杂，增加了乘法次数及储存量；且所有的谐振器的极点都是在单

?k位圆上，由WN决定。考虑到系数量化的影响，当系数量化时，极点会移动，有些

极点就不能被梳状滤波器的零点锁抵消。(零点由延时单元绝对，不受量化的影响)系统就会变得不稳定[1]。

王杰《音乐信号滤波去噪--使用GAUSS窗设计的频率采样型FIR滤波器》第8页共23页

3设计步骤

3.1设计流程图

滤波器设计流程图如图3-1所示。

图3-1 滤波器设计流程图

3.2下载并截取音乐信号

从网上下载一段音乐，从中截取一段格式为.wav 的语音信号，时间为2—3秒，并加入噪声。在MALAB 平台上，观察原始语音信号与加入噪声后的时域和频谱图。

源程序如下所示：

[x,fs,bits]=wavread('G:\\数字信号处理课程设计_rE\\funk.wav');%读取音乐信号的数据 sound(x,fs,bits);

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

王杰《音乐信号滤波去噪--使用GAUSS窗设计的频率采样型FIR滤波器》第9页共23页

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

fn=2436; %单频噪声功率 y=x+0.1*sin(fn*2*pi*t); %加单频噪声

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; %计算频谱的频率范围

原始信号与含噪信号的时域图和频域图如图3-2所示，由时频域分析图可以看出，音乐信号能量主要集中在1000Hz以内。

图3-2 原始信号与含噪信号的时频域分析图

3.3滤波器设计

截取好原始信号，接下来的工作是设计一个Gaussian窗滤波器。 利用公式[1]：

RP=-min(db([1:1:ceil(wp1/dww)+1,ceil(wp2/dww)+1:1:end])) (3-1)

王杰《音乐信号滤波去噪--使用GAUSS窗设计的频率采样型FIR滤波器》第10页共23页

AS=-max(db(wsd/dww+1:wsu/dww+1)) (3-2)

分别得到通带波纹及通带最小阻带衰减，并输入相应滤波器的频率参数，接着用freqz_m求得滤波器的响应频率特性。

源程序如下所示:

fp1=1836;fs1=2406;fs2=2476;fp2=3036; %带阻滤波器设置指标 Rp=6;As=28;

fcd=(fp1+fs1)/2; %计算上下边带中心频率 fcu=(fp2+fp2)/2;

df=min((fs1-fp1),(fp2-fs2));

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

M=ceil(5.8*pi/dw)+1; %计算Gauss窗设计该滤波器时需要的阶数 n=0:M-1; % 定义时间范围 gausswf=gausswin(M); %利用gausswin函数产生一个Gaussian窗函数 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); %计算理想阻带滤波器的脉冲响应

h_bs=gausswf'.*hd_bs; %窗口法计算实际滤波脉冲响应 [db,mag,pha,grd,w]=freqz_m(h_bs,[1]); %计算滤波器的频率特性 wp1=fp1/fs*2*pi;wp2=fp2/fs*2*pi;dww=2*pi/1000;

RP=-min(db([1:1:ceil(wp1/dww)+1,ceil(wp2/dww)+1:1:end])) %计算通带最大衰减 AS=-max(db(wsd/dww+1:wsu/dww+1)) %计算出阻带最小衰减

滤波器幅度响应图、相位响应图、脉冲响应图如图3-3所示。由图可知，滤波器的实际通带最大衰减高于设定的通带最大衰减，而实际的阻带最小衰减低于设定的阻带最小衰减，性能达标。