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

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

subplot(2,2,4);stem(n,h_bs);xlabel('n');ylabel('h(n)');title('滤波器脉冲响应图');grid on; y_fil=fftfilt(h_bs,y); %用设计好的滤波器对y进行滤波 Y_fil=abs(fft(y_fil));Y_fil=Y_fil(1:length(Y_fil)/2);%计算频谱取前一半

sound (y_fil,fs,bits); %音乐信号回放 wavwrite(y_fil,fs,'d:\\funk_filted.wav'); %产生去噪后的音乐文件 figure(3); %作图 subplot(4,2,1);plot(t,x);grid on;

title('原始音乐信号');xlabel('时间（t）');ylabel('幅度');axis([0,2.5,-2,2]); subplot(4,2,2);plot(f,X);grid on;

axis([0,3000,0,1500]);title('原始音乐信号幅度谱');xlabel('频率（f）');ylabel('幅度谱'); subplot(4,2,5);plot(t,y_fil);grid on;

axis([0,2.5,-2,2]);title('滤波后的音乐信号');xlabel('时间（t）');ylabel('幅度'); subplot(4,2,6);plot(f,Y_fil);grid on;

axis([0,3000,0,1500]);title('滤波后的音乐信号幅度谱');xlabel('频率（f）');ylabel('幅度谱');

subplot(4,2,3),plot(t,y);grid on;

axis([0,2.5,-2,2]);title('加入噪声干扰后的音乐信号');xlabel('时间（单位：s）');ylabel('幅度');

subplot(4,2,4),plot(f,Y);grid on;

axis([0,3000,0,1500]);title('加入噪声干扰后的音乐信号幅度谱图');xlabel('频率（单位：Hz）');ylabel('幅度谱');

subplot(4,2,7);plot(f,20*log10(abs(Y)/max(abs(X))));grid on;

title('加干扰后的音乐信号幅度谱');xlabel('频率（f）');ylabel('幅度对数坐标'); subplot(4,2,8);plot(f,20*log10(abs(Y_fil)/max(abs(X))));grid on;

axis([0,4000,-100,0]);title('滤波后的音乐信号幅度谱');xlabel('频率（f）');ylabel('幅度对数坐标');

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

附录2：滤波器直接型结构转换频率采样型源程序清单

% 程序名称： tf2fs.m

% 程序功能：将FIR直接型结构的h(n)转换成频率采样结构，得到其各个参数的子程序。

% 程序作者：王杰

% 最后修改日期： 2017-2-28 function [C,B,A]=tf2fs(h) %直接型到频率采样型的转换 %---------------------- %[C,B,A]=tf2fs(h)

%C=包含各并联部分增益的行向量 %B=包含按行排列的分子系数矩阵 %A=包含按行排列的分母系数矩阵 %h=FIR滤波器单位冲激响应向量 %---------------------- N=length(h); H=fft(h,N); r=input('r=');

magH=abs(H);phaH=angle(H)'; %check even or odd N if (N==2*floor(N/2)) L=N/2-1; A1=[1,-r,0;1,r,0];

C1=[real(H(1)),real(H(L+2))]; else

L=(N-1)/2; A1=[1,-r,0]; C1=[real(H(1))]; end

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

k=[1:L]';

%初始数组B和数组A B=zeros(L,2);A=zeros(L,3); %计算分母系数 A(1:L,1)=1;

A(1:L,2)=-2*r*cos(2*pi*k/N); A(1:L,3)=r.^2 A=[A;A1]; %计算分子系数

B(1:L,1)=cos(phaH(2:L+1));

B(1:L,2)=-cos(phaH(2:L+1)-(2*pi*k/N)); %计算增益系数 C=[2*magH(2:L+1),C1]';