数字信号处理的课程设计---数字滤波器的设计：歌曲中的人声抑制

第1页共21页实验报告数字滤波器的设计课程名称数字信号处理题目名称_歌曲中的人声抑制学生学院专业班级学号学生姓名指导教师2012年7月1日第2页共21页广东工业大学课程设计任务书题目名称歌曲中的人声抑制学生学院专业班级姓名学号一、课程设计的内容运用数字信号处理的理论知识进行频谱分析和滤波器设计。先通过理论分析和推导设计出符合技术指标要求的数字滤波器，并推断出相应的结论；再利用MATLAB或C/C++作为编程工具进行数字滤波的计算机实现。二、课程设计的要求与数据1.选择一段带有演唱人声的高保真歌曲文件，利用某种格式转换软件将其转换为单声道的.wav文件，并记录下采样频率等参数。用MATLAB或C/C++编程显示出该.wav文件的时域波形图和频谱图。2.针对FIR和IIR数字滤波器，各选择一种，根据以下给定的技术指标设计数字滤波器得出滤波器系数（给出理论计算过程）。并显示出各滤波器的频率响应。滤波器技术指标类型：带阻滤波器3dB截止频率：fc1=500Hz，fc2=2700Hz通带内幅频响应波纹起伏：2dB阻带截止频率：fs1=650Hz，fs2=1900Hz阻带内幅频响应衰减：-20dB3.用自己设计的滤波器对以上.wav文件进行滤波，显示出滤波后信号的时域波形图和频谱图，并从主观听觉和客观数据上对滤波前后的信号进行对比，分析信号的变化。三、课程设计应完成的工作运用数字信号处理的理论知识进行频谱分析和滤波器设计。先通过理论分析和推导设计出符合技术指标要求的数字滤波器，并推断出相应的结论；再利用MATLAB或C/C++作为编程工具进行数字滤波的计算机实现。第3页共21页四、应收集的资料及主要参考文献1、赵健，李勇，数字信号处理，清华大学出版社，北京，2010发出任务书日期：2012年6月12日指导教师签名：计划完成日期：2012年7月1日基层教学单位责任人签章：主管院长签章：第4页共21页目录1.实验目的与要求——————————————————————————————52.实验方案———————————————————————————————53.实验结果及数据处理————————————————————————————204.结论—————————————————————————————————235.问题与讨论——————————————————————————————23参考文献—————————————————————————————————24附录———————————————————————————————————24第5页共21页一、实验目的与要求1．实验目的高保真歌曲可涵盖20-20000Hz频率范围，而其中的演唱人声的人耳敏感内容则只占一部分中间频段。我们可以利用带阻滤波器来对歌曲中的演唱人声加以抑制，从而获得一个伴唱版的歌曲。本实验综合运用数字信号处理的理论知识进行频谱分析和滤波器设计。先通过理论分析和推导设计出符合技术指标要求的数字滤波器，并推断出相应的结论；再利用MATLAB或C/C++作为编程工具进行数字滤波的计算机实现。1)掌握数字信号处理的基本概念、基本理论和基本方法；2)掌握用MATLAB或C/C++设计FIR和IIR数字滤波器的方法；3)实现对一个音频信号进行数字滤波的完整过程。2.实验要求1.选择一段带有演唱人声的高保真歌曲文件，利用某种格式转换软件将其转换为单声道的.wav文件，并记录下采样频率等参数。用MATLAB或C/C++编程显示出该.wav文件的时域波形图和频谱图。2.针对FIR和IIR数字滤波器，各选择一种，根据以下给定的技术指标设计数字滤波器得出滤波器系数（给出理论计算过程）。并显示出各滤波器的频率响应。滤波器技术指标类型：带阻滤波器3dB截止频率：fc1=500Hz，fc2=2700Hz通带内幅频响应波纹起伏：2dB阻带截止频率：fs1=650Hz，fs2=1900Hz阻带内幅频响应衰减：-20dB3.用自己设计的滤波器对以上.wav文件进行滤波，显示出滤波后信号的时域波形图和频谱图，并从主观听觉和客观数据上对滤波前后的信号进行对比，分析信号的变化。二、实验方案1.转换格式。2.Matlab的功能十分强大，已经为我们提供了很多设计滤波器的函数。下面列举一些下面需要用到的函数。（一）。语音的录入与打开：①[y,fs,bits]=wavread('Blip',[N1N2]);用于读取语音，采样值放第6页共21页在向量y中，fs表示采样频率(Hz)，bits表示采样位数。[N1N2]表示读取从N1点到N2点的值（若只有一个N的点则表示读取前N点的采样值）。②sound(x,fs,bits);用于对声音的回放。向量y则就代表了一个信号（也即一个复杂的“函数表达式”）也就是说可以像处理一个信号表达式一样处理这个声音信号。③函数FFT用于序列快速傅立叶变换。函数的一种调用格式为y=fft(x)其中，x是序列，y是序列的FFT，x可以为一向量或矩阵，若x为一向量，y是x的FFT。且和x相同长度。若x为一矩阵，则y是对矩阵的每一列向量进行FFT。④函数abs(x)用于计算复向量x的幅值，函数angle(x)用于计算复向量的相角，介于和之间，以弧度表示。（二）、滤波器设计：1、相关原理：设计数字滤波器的任务就是寻求一个因果稳定的线性时不变系统，并使系统函数H(z)具有指定的频率特性。数字滤波器从实现的网络结构或者从单位冲激响应分类，可以分成无限长单位冲激响应（IIR）数字滤波器和有限长单位冲激响应（FIR）数字滤波器。数字滤波器频率响应的三个参数：（1）幅度平方响应：第7页共21页（2）相位响应（3）群时延响应IIR数字滤波器：IIR数字滤波器的系统函数为的有理分数，即IIR数字滤波器的逼近问题就是求解滤波器的系数和，使得在规定的物理意义上逼近所要求的特性的问题。如果是在s平面上逼近，就得到模拟滤波器，如果是在z平面上逼近，则得到数字滤波器。FIR数字滤波器：设FIR的单位脉冲响应h(n)为实数，长度为N，则其z变换和频率响应分别为按频域采样定理FIR数字滤波器的传输函数H(z)和单位脉冲响应h(n)可由它的N个频域采样值H(k)唯一确定。MATLAB中提供了几个函数，分别用于实现IIR滤波器和FIR滤波器。（1）卷积函数conv卷积函数conv的调用格式为c=conv(a,b)该格式可以计算两向量a和b的卷积，可以直接用于对有限长信号采用FIR滤波器的滤波。（2）函数filter函数filter的调用格式为y=filter(b,a,x)第8页共21页该格式采用数字滤波器对数据进行滤波，既可以用于IIR滤波器，也可以用于FIR滤波器。其中向量b和a分别表示系统函数的分子、分母多项式的系数，若a＝1，此时表示FIR滤波器，否则就是IIR滤波器。该函数是利用给出的向量b和a，对x中的数据进行滤波，结果放入向量y。（3）函数fftfilt函数fftfilt的调用格式为y=fftfilt(b,x)该格式是利用基于FFT的重叠相加法对数据进行滤波，这种频域滤波技术只对FIR滤波器有效。该函数是通过向量b描述的滤波器对x数据进行滤波。（4）关于用butter函数求系统函数分子与分母系数的几种形式。[b,a]=butter(N,wc,'high'):设计N阶高通滤波器，wc为它的3dB边缘频率，以为单位，故。[b,a]=butter(N,wc):当wc为具有两个元素的矢量wc=[w1,w2]时，它设计2N阶带通滤波器，3dB通带为，w的单位为。[b,a]=butter(N,wc,'stop'):若wc=[w1,w2]，则它设计2N阶带阻滤波器，3dB通带为，w的单位为。如果在这个函数输入变元的最后，加一个变元“s”，表示设计的是模拟滤波器。这里不作讨论。为了设计任意的选项巴特沃斯滤波器，必须知道阶数N和3dB边缘频率矢量wc。这可以直接利用信号处理工具箱中的buttord函数来计算。如果已知滤波器指标，，和，则调用格式为第9页共21页[N,wc]=buttord(wp,ws,Rp,As)对于不同类型的滤波器，参数wp和ws有一些限制：对于低通滤波器，wpws；对于高通滤波器，wpws；对于带通滤波器，wp和ws分别为具有两个元素的矢量，wp=[wp1,wp2]和ws=[ws1,ws2]，并且ws1wp1wp2ws2；对于带阻滤波器wp1ws1ws2wp2。（5）窗函数滤波器，即FIR有限冲激相应滤波器，这样取得：b=fir1（n，Wn,'stop'）,n是一个合适的阶数，适当选取，要求不高，Wn上面取得了，stop表示带阻。y=filter2（b，x）x为滤波器输入，y即输出。(6)Bartlett函数：生成巴特利特窗n=51;window=bartlett(n);[h,w]=freqz(window,1);subplot(1,2,1)stem(window);subplot(1,2,2);plot(w/pi,20*log(abs(h)/abs(h(1))));3.记录下采样频率等参数。用MATLAB或C/C++编程显示出该.wav文件的时域波形图和频第10页共21页谱图。[x,fs,bits]=wavread('slowsoul.wav',[10245120]);sound(x,fs,bits);X=fft(x,44100);magX=abs(X);angX=angle(X);subplot(221);plot(x);title('原始信号波形');subplot(222);plot(X);title('原始信号频谱');subplot(223);plot(magX);title('原始信号幅值');subplot(224);plot(angX);title('原始信号相位');原始信号波形记录：采样频率为：44100HZ4.进行IIR滤波器的设计。MATLAB强大的功能已经具备IIR设计的函数，所以我们可以直接调第11页共21页用得出滤波器的各种参数及结果。首先确定IIR滤波器的频率响应：[x,fs,nbits]=wavread('slowsoul.wav');fp1=500;%通带频率fp2=2700;%通带频率fc1=650;%阻带频率fc2=1900;%阻带频率wp1=2*fp1/fs;%采样频率归一化的通带角频率ws1=2*fc1/fs;%采样频率归一化的阻带角频率wp2=2*fp2/fs;%采样频率归一化的通带角频率ws2=2*fc2/fs;%采样频率归一化的阻带角频率wp=[wp1,wp2];ws=[ws1,ws2];rp=2;%通带波动系数rs=20;%阻带波动系数[n,wc]=buttord(wp,ws,rp,rs);%求巴特沃兹滤波器的系数及3db截止频率[b,a]=butter(n,wc,'stop');%计算滤波器系数[h,f]=freqz(b,a,2*fs,fs);mag=abs(h);ph=angle(h);ph=ph*180/pi;subplot(2,1,1);plot(f,mag);gridon;xlabel('频率(Hz)');ylabel('幅频特性');subplot(2,1,2);plot(f,ph);gridon;xlabel('频率(Hz)');ylabel('相频特性');结果：第12页共21页接着进行滤波，观察滤波后波形。代码：[x,fs,nbits]=wavread('slowsoul.wav');fp1=500;%通带频率fp2=2700;%通带频率fc1=650;%阻带频率fc2=1900;%阻带频率wp1=2*fp1/fs;%采样频率归一化的通带角频率ws1=2*fc1/fs;%采样频率归一化的阻带角频率wp2=2*fp2/fs;%采样频率归一化的通带角频率ws2=2*fc2/fs;%采样频率归一化的阻带角频率wp=[wp1,wp2];ws=[ws1,ws2];rp=2;%通带波动系数rs=20;%阻带波动系数