切比雪夫I滤波器

%程序名称：qbxf.m%程序功能：采用基于切比雪夫设计法，设计IIR滤波器对含噪语音进行滤波去噪处理。%程序作者：%最后修改日期：2014-3-3[x,fs,bits]=wavread('D:\aaa\zhouyi1.wav');%输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。%sound(x,fs,bits);%按指定的采样率和每样本编码位数回放N=length(x);%计算信号x的长度fn=2400;%单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs;%计算时间范围，样本数除以采样频率x=x';y=x+0.1*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);%截取前半部分deltaf=fs/N;%计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf;%计算频谱频率范围figure(1)subplot(2,2,1);plot(t,x);gridon;ylabel('幅度');xlabel('时间(单位：s)'),title('原始语言信号');subplot(2,2,2);plot(f,X);gridon;axis([030000500]);%xmin,xmax分别设定了图中x轴上的最小值和最大值ylabel('幅度谱');xlabel('频率(单位：Hz)');title('语言信号幅度频谱图');subplot(2,2,3);plot(t,y);gridon;axis([02-11]);%xmin,xmax分别设定了图中x轴上的最小值和最大值ylabel('幅度');xlabel('时间(单位：s)');title('加入单频干扰后的语言信号');subplot(2,2,4);plot(f,Y);gridon;axis([030000500]);%xmin,xmax分别设定了图中x轴上的最小值和最大值ylabel('幅度谱');xlabel('频率(单位：Hz)');title('加入干扰后的语言信号幅度谱图');fp=fn-400;fc=fn-300;%定义通带和阻带截止频率Rp=1;As=16;%定义通带波纹和阻带衰减wp=2*pi*fp/fs;ws=2*pi*fc/fs;%计算对应的数字频率T=1;%定义采样间隔%双线性变换法设计%第一步预畸变Omegap=(2/T)*tan(wp/2);Omegas=(2/T)*tan(ws/2);%截止频率预畸变%第二步根据指标设计模拟滤波器[c,d]=afd_chb1(Omegap,Omegas,Rp,As);%计算滤波器阶数和截止频率%第三步模数滤波器变换[b,a]=bilinear(c,d,1/T);%双线性变换得到数字滤波器系统函数分子分母系数%设计结束[db,mag,pha,grd,w]=freqz_m(b,a);%验证滤波器是否达到指定性能%计算脉冲响应delta=[1,zeros(1,99)];ha=filter(b,a,delta);figure(2)Subplot(2,2,1);plot(w/pi,db);axis([0.40.6-202]);%xmin,xmax分别设定了图中x轴上的最小值和最大值%绘制相对幅度响应曲线图，限制横纵坐标范围，添加网格X_l=[0.4,0.4,wp/pi,ws/pi;0.6,0.6,wp/pi,ws/pi];Y_l=[-As,-Rp,-20,-20;-As,-Rp,2,2];%在wp，ws，Rp，As处画线以更直观判断设计是否达标，每列参数是每个线条的端点坐标line(X_l,Y_l,'Color','r','LineWidth',2,'LineStyle','--')%添加线条，红色，线宽为2title('相对振幅');xlabel('w/pi');ylabel('dB');Subplot(2,2,2);plot(w/pi,mag);axis([01-0.21.2]);%xmin,xmax分别设定了图中x轴上的最小值和最大值title('绝对振幅');xlabel('w/pi');ylabel('mag');Subplot(2,2,3);plot(w/pi,pha);gridon;axis([01-44]);%xmin,xmax分别设定了图中x轴上的最小值和最大值title('相位响应');xlabel('w/pi');ylabel('pha');Subplot(2,2,4);plot(ha);gridon;axis([060-0.20.4]);%xmin,xmax分别设定了图中x轴上的最小值和最大值title('脉冲响应');xlabel('时间');ylabel('幅度');y_fil=filter(b,a,y);%IIR滤波器对信号signal进行滤波处理，其中b,a为上面设计好的滤波器参数。Y_fil=abs(fft(y_fil));Y_fil=Y_fil(1:N/2);%计算频谱取前一半sound(y_fil,fs,bits);%按指定的采样率和每样本编码位数回放figure(3)subplot(3,2,1);plot(t,x);axis([01.75-11]);%xmin,xmax分别设定了图中x轴上的最小值和最大值ylabel('幅度');xlabel('时间(单位：s)'),title('原始语言信号');subplot(3,2,2);plot(f,X);axis([030000500]);%xmin,xmax分别设定了图中x轴上的最小值和最大值ylabel('幅度谱');xlabel('频率(单位：Hz)');title('语言信号幅度频谱图');subplot(3,2,3);plot(t,y);axis([01.75-11]);%xmin,xmax分别设定了图中x轴上的最小值和最大值ylabel('幅度');xlabel('时间(单位：s)');title('加入单频干扰后的语言信号');subplot(3,2,4);plot(f,Y);axis([030000500]);%xmin,xmax分别设定了图中x轴上的最小值和最大值ylabel('幅度谱');xlabel('频率(单位：Hz)');title('加入干扰后的语言信号幅度谱图');subplot(3,2,5);plot(t,y_fil);axis([01.75-11]);%xmin,xmax分别设定了图中x轴上的最小值和最大值ylabel('幅度');xlabel('时间(单位：s)'),title('滤波后语言信号');subplot(3,2,6);plot(f,Y_fil);axis([030000500]);%xmin,xmax分别设定了图中x轴上的最小值和最大值ylabel('幅度谱');xlabel('频率(单位：Hz)');title('滤波后语言信号幅度频谱图');wavwrite(y,fs,bits,'zhouyi2.wav');%输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。wavwrite(y_fil,fs,bits,'zhouyi3.wav');%输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。[C,B,A]=dir2cas(b,a);CBA