[1]张志涌．精通MATLAB 6.5版[M]．北京：北京航空航天大学出版社，2003． [2]约翰·Ｇ·普罗克斯．数字信号处理[M]．西安：西安交通大学出版社，2009． [3]张小虹．信号系统与数字信号处理[M]．第1版．西安：西安电子科技出版社，2002． [4]谢德芳．数字信号处理[M]．北京：科学出版社，2005．
[5]郝建军．数字通信[M]．第2版．北京：北京邮电大学出版社，2010．
[6]张威．MATLAB 基础与编程入门[M]．西安：西安电子科技大学出版社，2010．
附录一： 语音信号滤波去噪——使用TUKEYWIN滤波器
第 11 页
%程序功能：在Matlab中，用窗口设计法设计FIR滤波器 %程序作者： %最后修改日期：
%程序一：
[x,fs,bits]=wavread('speech_dftwav');%读取音乐信号 sound(x,fs,bits);%播放音乐信号 N=length(x);%计算音乐信号的长度 fn=2136; % 单频噪声频率
t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率 x=x(:,1)'; % 将双声道转为单声道 y=x+sin(fn*2*pi*t);%添加噪声
sound(y,fs,bits);%播放干扰后的音乐信号，应该可以明显听出有尖锐的单频啸叫声
%程序二：
X=abs(fft(x)); Y=abs(fft(y)); % 对原始信号和加噪信号进行fft变换,取幅度谱 X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分 figure(1)
subplot(2,2,1);plot(t,x);xlabel('时间(t)');ylabel('幅度'); title('原始语音信号');axis([0,2.5,-0.7,0.7]);
subplot(2,2,2);plot(f,X);xlabel('频率(f)');ylabel('幅度谱');title('原始语音信号幅度谱');axis([0,4000,0,300]);
subplot(2,2,3);plot(t,y);xlabel('时间(t)');ylabel('幅度');title('加噪声后的语音信号');axis([0,2.5,-0.7,0.7]);
subplot(2,2,4);plot(f,Y);xlabel('频率(f)');ylabel('幅度谱');title('加噪声后的语音信号幅度谱');axis([0,4000,0,300]);
%程序三：
fpd=2000;fsd=2050;fsu=2150;fpu=2200;% 带阻滤波器设计指标 Rp=1;As=22; % 通带波纹和阻带衰减
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(2.4*pi/dw)+1; % 计算Tukeywin窗设计该滤波器时需要的阶数
第 12 页
n=0:M-1; % 定义时间范围
w_tuk=tukeywin(M); % 产生M阶的Tukeywin窗 hd_bs=ideal_lp(wcd,M)+ideal_lp(pi,M)-ideal_lp(wcu,M); % 调用自编函数计算理想带阻滤波器的脉冲响应
h_bs=w_tuk'.*hd_bs; % 用窗口法计算实际滤波器脉冲响应
[db,mag,pha,grd,w]=freqz_m(h_bs,1); % 调用自编函数计算滤波器的频率特性 figure(2)
subplot(2,2,1);plot(w/pi,db);xlabel('w/pi');ylabel('db');title('FIR滤波器的幅度响应图');axis([0.3 0.7 -70 10]);
line([0,0.7],[-1,-1],'Color','r','LineWidth',2,'LineStyle','--'); line([0,0.7],[-22,-22],'Color','r','LineWidth',2,'LineStyle','--');
line([wsd/pi,wsd/pi],[-70,10],'Color','r','LineWidth',2,'LineStyle','--'); line([wsu/pi,wsu/pi],[-70,10],'Color','r','LineWidth',2,'LineStyle','--'); grid on;
subplot(2,2,2);plot(w/pi,mag);xlabel('w/pi');ylabel('幅度mag');title('FIR滤波器的幅度响应图');axis([0 1 -0.2 1.2]);grid on;
subplot(2,2,3);plot(w/pi,pha);xlabel('w/pi');ylabel('相位pha');title('滤波器相位响应图');axis([0 1 -4 4]);grid on;
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:N/2); % 计算频谱取前一半 sound (y_fil,fs,bits);
wavwrite(y_fil,fs,bits,'1111恢复.wav'); %保存恢复的语音信号 wavwrite(y,fs,bits,'1111加噪.wav'); %保存加噪的语音信号 figure(3)
subplot(3,2,1);plot(t,x);xlabel('时间(t)');ylabel('幅度'); title('原始语音信号');axis([0,2.5,-0.7,0.7]);
subplot(3,2,2);plot(f,X);xlabel('频率(f)');ylabel('幅度谱');title('原始语音信号幅度谱');axis([0,4000,0,300]);
subplot(3,2,3);plot(t,y);xlabel('时间(t)');ylabel('幅度');title('加噪声后的语音信号');axis([0,2.5,-0.7,0.7]);
subplot(3,2,4);plot(f,Y);xlabel('频率(f)');ylabel('幅度谱');title('加噪声后的语音信号幅度谱');axis([0,4000,0,300]);
subplot(3,2,5);plot(t,y_fil);xlabel('时间(t)');ylabel('幅度'); title('滤波后语音信号');axis([0,2.5,-0.7,0.7]);
subplot(3,2,6);plot(f,Y_fil);xlabel('频率(f)');ylabel('幅度谱');title('滤波后语音信号幅度谱');axis([0,4000,0,300]); 附录二：
函数FREQZ_M.M定义：
function [db,mag,pha,grd,w] = freqz_m(b,a);
第 13 页
% 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; 附录三：
函数IDEAL_LP.M定义： function hd = ideal_lp(wc,M); % 理想低通滤波器计算
第 14 页
% -------------------------------- % [hd] = ideal_lp(wc,M)
% hd = 0 to M-1之间的理想脉冲响应 % wc = 截止频率(弧度) % M = 理想滤波器的长度 %
alpha = (M-1)/2; n = [0:1:(M-1)]; m = n - alpha + eps; hd = sin(wc*m) ./ (pi*m);
第 15 页