音频处理Matlab程序

				[x,fs,nbit]=wavread('Test.wav');%读取Test音频文件
				sound(x,fs)
				N=length(x);
				n=1:N;
				c=0.1*sin(2*pi*1000*n./fs);
				x1=x'+c;%将正弦波加载于原音频文件
				sound(x1,fs);
				subplot(2,2,1);stem(n,x);
				xlabel('采样点');
				ylabel('模拟输入 (V)');
				title('原始声音信号'),xlabel('n')
				X=fft(x);magX=abs(X(1:(N+1)/2+1));
				k=0:1:(N+1)/2;w=2*pi/N*k;
				subplot(2,2,2);stem(w/pi,magX);
				title('原始声音信号幅频图');xlabel('以pi为单位的频率')
				subplot(2,2,3);stem(n,x1);
				xlabel('采样点');
				ylabel('模拟输入 (V)');
				title('加入噪声的声音信号');xlabel('n')
				X1=fft(x1);magX1=abs(X1(1:(N+1)/2+1));
				subplot(2,2,4);stem(w/pi,magX1);
				title('加入噪声信号幅频图');xlabel('以pi为单位的频率')
				wp=0.02*pi;%通带截止频率
				ws=0.04*pi;%阻带截止频率
				Rp=1;%通带纹波
				As=20;%阻带纹波
				T=1;
				OmegaP=(2/T)*tan(wp/2);%原型通带频率
				OmegaS=(2/T)*tan(ws/2);%原型阻带频率
				ep=sqrt(10^(Rp/10)-1);%通带纹波参数
				Ripple=sqrt(1/(1+ep*ep));%通带纹波
				Attn=1/(10^(As/20));%阻带衰减
				[cs,ds]= p_butt(OmegaP,OmegaS,Rp,As);
				[b,a]= bilinear(cs,ds,T);%双线性变换
				[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);
				grd=grpdelay(b,a,w);
				y=filter(b,a,x1);%b,a为巴特沃兹低通滤波器分子、分母系数
				Y1=fft(y);magY1=abs(Y1(1:(N+1)/2+1));k=0:1:(N+1)/2;w=2*pi/N*k;
				stem(w/pi,magY1);title('带噪声信号的滤波后的幅频图')
				sound(y,fs);