目录
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调音符与频率对照表
西电《数字信号处理》第三版