语音信号滤波去噪——使用flattopwin窗设计的FIR滤波器要点

黄珺 《吹管滤波去噪—使用FLATTOPWIN窗设计的FIR滤波器》 第17页，共21页

附录一： 音乐信号滤波去噪——使用FLARTOPWIN滤波器

%程序功能：在Matlab中，用窗口设计法设计FIR滤波器 %程序作者：赵伯政 %最后修改日期：2012-3-4

%程序一：

>> [x,fs,bits]=wavread(yinyue.wav'');%读取音乐信号 >> sound(x,fs,bits);%播放音乐信号 >> N=length(x);%计算音乐信号的长度 >> fn=2200; % 单频噪声频率

>> t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率 >> x=x(:,1)'; y=x+sin(fn*2*pi*t);

>> plot(t,x); xlabel('t ins');ylabel('x');%画原始音乐信号的时域波形图

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(2)

>> subplot(2,2,1);plot(t,x);xlabel('时间(t)');ylabel('幅度'); title('原始音乐信号'); >> subplot(2,2,2);plot(f,X);xlabel('频率(f)');ylabel('幅度谱');title('原始音乐信号幅度谱');axis([0,8000,0,600]);

>> subplot(2,2,3);plot(t,y);xlabel('时间(t)');ylabel('幅度');title('加干扰后的音乐信号'); >> subplot(2,2,4);plot(f,Y);xlabel('频率(f)');ylabel('幅度谱');title('加干扰后的音乐信号幅度谱');axis([0,8000,0,600]);

%程序三：

黄珺 《吹管滤波去噪—使用FLATTOPWIN窗设计的FIR滤波器》 第18页，共21页

>> fpd=2100;fsd=2150;fsu=2250;fpu=2300;Rp=1;As=20;% 带阻滤波器设计指标 >> 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.1*pi/dw)+1; % 计算窗设计该滤波器时需要的阶数 >> n=0:M-1; % 定义时间范围

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

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

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

subplot(2,2,1);plot(w,db);title('滤波器幅度响应图');xlabel('w/pi');ylabel('db'); axis([0,0.5,-30,10]);

>> line([0,0.5],[-As,-As],'color','r','linestyle','--','LineWidth',2); >> line([0,0.5],[-Rp,-Rp],'color','r','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,3,-4,4]);

>> subplot(2,2,4);stem(n,h_bs);title('滤波器脉冲响应图');xlabel('n');ylabel('h(n)'); axis([0,3000,0,1]);

%程序四：

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

黄珺 《吹管滤波去噪—使用FLATTOPWIN窗设计的FIR滤波器》 第19页，共21页

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

>> subplot(3,2,1);plot(t,x);xlabel('时间(t)');ylabel('幅度'); title('原始音乐信号'); >> subplot(3,2,2);plot(f,X);xlabel('频率(f)');ylabel('幅度谱');title('原始音乐信号幅度谱');axis([0,8000,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,8000,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,8000,0,600]);

黄珺 《吹管滤波去噪—使用FLATTOPWIN窗设计的FIR滤波器》 第20页，共21页

附录二：

函数FREQZ_M.M定义：

function [db,mag,pha,grd,w] = freqz_m(b,a); % freqz 子程序的改进版本 % ------------------------------------ % [db,mag,pha,grd,w] = freqz_m(b,a);

% db = [0 到pi弧度]区间内的相对振幅(db) % mag = [0 到pi弧度]区间内的绝对振幅 % pha = [0 到pi弧度]区间内的相位响应 % grd = [0 到pi弧度]区间内的群迟延

% w = [0 到pi弧度]区间内的501个频率样本向量 % b = Ha(z)的分子多项式系数(对FIR b=h) % a = Ha(z)的分母多项式系数(对 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;