数字信号处理之低通滤波器设计
目录
1 clc; 2 clear; 3 close all; 4 5 [x1,fs,bits]=wavread(\'qq.wav\'); %读取语音信号的数据,赋给变量x1 6 L=length(x1); 7 sound(x1,fs); %播放语音信号 8 9 N=L; 10 k=[0:N-1]; 11 y1=fft(x1,N); %对信号做L个点FFT变换 12 figure(1); 13 subplot(2,1,1); 14 plot(k(1:20000),abs(y1(1:20000))); %做原始语音信号的FFT频谱图 15 title(\'原始语音信号FFT频谱\'); 16 17 f=fs*(0:N-1)/N; 18 subplot(2,1,2); 19 plot(f(1:20000),abs(y1(1:20000))); 20 title(\'原始语音信号频谱\'); 21 xlabel(\'f Hz\');ylabel(\'fuzhi\');
说明:程序第17行的解释见“四、DFT的周期性”公式四。
说明:上图是FFT变换的频谱图(数字频率k),下图是转换成模拟频率下的频谱图。
1 %给原始的语音信号加上一个高频余弦噪声,频率为5kHz。画出加噪后的语音信号时域和频谱图,与原始信号对比,可以很明显的看出区别。 2 t=[0:1/fs:(L-1)/fs]; %将所加噪声信号的点数调整到与原始信号相同,构造采样时间点(模拟时间) 3 Au=0.03; %噪声幅度 4 d=[Au*cos(2*pi*8000*t)]\'; %噪声为8kHz(2*pi*8000/(2*pi))的余弦信号(模拟时间) 5 x2=x1+d; 6 7 figure(2); 8 subplot(2,1,1); 9 plot(x1); %做原始语音信号的时域图形 10 title(\'原始语音信号\');xlabel(\'time n\');ylabel(\'fuzhi\'); 11 12 sound(x2,fs); %播放加噪声后的语音信号 13 subplot(2,1,2); 14 plot(x2); %做原始语音信号的时域图形 15 title(\'加噪后的信号\');xlabel(\'time n\');ylabel(\'fuzhi\'); 16 17 y2=fft(x2,N); %对信号做L个点FFT变换 18 figure(3); 19 subplot(2,1,1); 20 plot(k(1:20000),abs(y2(1:20000))); %做原始语音信号的FFT频谱图 21 title(\'加躁语音信号FFT频谱\'); 22 23 subplot(2,1,2); 24 plot(f(1:20000),abs(y2(1:20000))); 25 title(\'加躁语音信号频谱\'); 26 xlabel(\'f Hz\');ylabel(\'fuzhi\');
说明:程序前5行是构造高频杂波信号,高频噪声信号也是按照fs为采样频率的,而且采样点数也是L(与原始信号等长)。
说明:通过频谱图就可以看到所加噪声信号是8Khz。
1 wp=2*pi*4000; %通带边界角频率 2 ws=2*pi*5000; %阻带边界角频率 3 Rp=1; %通带最大衰减 4 Rs=15; %阻带最小衰减 5 [NN,Wn]=buttord(wp,ws,Rp,Rs,\'s\'); %选择滤波器的最小阶数 6 [Z,P,K]=buttap(NN); %创建butterworth模拟滤波器 7 [Bap,Aap]=zp2tf(Z,P,K); 8 [b,a]=lp2lp(Bap,Aap,Wn); 9 [bz,az]=bilinear(b,a,fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换 10 11 [H,W]=freqz(bz,az); %绘制频率响应曲线 12 figure(4); 13 plot(W*fs/(2*pi),abs(H)); 14 grid; 15 xlabel(\'频率/Hz\'); 16 ylabel(\'频率响应幅度\'); 17 title(\'Butterworth\');
说明:
(1)前4行是设置波特沃斯滤波器的指标
(2)第5行是根据前4行的参数计算滤波器的阶数NN和3db截止频率Wn。
(3)第6行根据阶数NN计算巴特沃斯归一化模拟低通原型滤波器系统函数的零、极点和增益因子。
(4)第9行将模拟滤波器转换到数字滤波器,而bz、az分别是分子和分母的系数。
到此,有了bz、az,滤波器就构造出来了。
(5)第11行计算滤波器的频率响应,freqz函数的使用见freqz()–matlab函数。
(6)第13行的f(模拟频率)=W*fs/(2*pi)解释见“2.2.2采样后的离散傅里叶频谱”。
说明:可以看到巴特沃斯滤波器的特性是与之前设置的参数匹配的。
x3=filter(bz,az,x2); figure(5); subplot(2,1,1); plot(t,x2); %画出滤波前的时域图 title(\'滤波前的时域波形\'); subplot(2,1,2); plot(t,x3); %画出滤波后的时域图 title(\'滤波后的时域波形\'); sound(x3,fs); %播放滤波后的信号 y2=fft(x2,N); %对信号做L个点FFT变换 figure(6); subplot(2,1,1); plot(f(1:20000),abs(y2(1:20000))); title(\'加躁语音信号频谱\'); xlabel(\'f Hz\');ylabel(\'fuzhi\'); y2=fft(x3,N); %对信号做L个点FFT变换 subplot(2,1,2); plot(f(1:20000),abs(y2(1:20000))); title(\'去噪语音信号频谱\'); xlabel(\'f Hz\');ylabel(\'fuzhi\');
说明:程序的第一行“x3=filter(bz,az,x2);”就是调用上边构造的巴特沃斯低通滤波器对含躁信号x2滤波,其实就是数值运算。
说明:通过对比可以看到,8Khz的高频噪声确实是被滤掉了。
附:C调音符与频率对照表
西电《数字信号处理》第三版