plot(f,abs(y1(1:512))); title('原始语音信号频谱') xlabel('Hz'); ylabel('fuzhi'); 程序2：
fs=22050; %语音信号采样频率为 22050
x1=wavread('Windows Critical Stop.wav'); %读取语音信号的数据，赋给变量 x1 t=0:1/22050:(size(x1)-1)/22050;
y1=fft(x1,1024); %对信号做 1024点 FFT变换 f=fs*(0:511)/1024;
x2=randn(1,length(x1)); %产生一与 x长度一致的随机信号 sound(x2,22050); figure(1)
plot(x2) %做原始语音信号的时域图形 title('高斯随机噪声'); xlabel('time n'); ylabel('fuzhi n'); randn('state',0); m=randn(size(x1)); x2=0.1*m+x1;
sound(x2,22050);%播放加噪声后的语音信号 y2=fft(x2,1024); figure(2) plot(t,x2)
title('加噪后的语音信号'); xlabel('time n'); ylabel('fuzhi n'); figure(3)
subplot(2,1,1);
plot(f,abs(y2(1:512))); title('原始语音信号频谱'); xlabel('Hz'); ylabel('fuzhi'); subplot(2,1,2);
plot(f,abs(y2(1:512)));
title('加噪后的语音信号频谱'); xlabel('Hz');
ylabel('fuzhi');
根据以上代码，你可以修改下面有错误的代码 程序3:双线性变换法设计 Butterworth滤波器 fs=22050;
x1=wavread('h:\\课程设计 2\\shuzi.wav'); t=0:1/22050:(size(x1)-1)/22050; Au=0.03;
d=[Au*cos(2*pi*5000*t)]'; x2=x1+d; wp=0.25*pi; ws=0.3*pi; Rp=1; Rs=15; Fs=22050; Ts=1/Fs;
wp1=2/Ts*tan(wp/2); %将模拟指标转换成数字指标 ws1=2/Ts*tan(ws/2);
[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s'); %选择滤波器的最小阶数 [Z,P,K]=buttap(N); %创建 butterworth模拟滤波器 [Bap,Aap]=zp2tf(Z,P,K); [b,a]=lp2lp(Bap,Aap,Wn);
[bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器 的转换
[H,W]=freqz(bz,az); %绘制频率响应曲线 figure(1)
plot(W*Fs/(2*pi),abs(H)) grid
xlabel('频率／Hz') ylabel('频率响应幅度') title('Butterworth') f1=filter(bz,az,x2); figure(2)
subplot(2,1,1)
plot(t,x2) %画出滤波前的时域图 title('滤波前的时域波形'); subplot(2,1,2)
plot(t,f1); %画出滤波后的时域图 title('滤波后的时域波形');
sound(f1,22050); %播放滤波后的信号 F0=fft(f1,1024);
f=fs*(0:511)/1024; figure(3)
y2=fft(x2,1024); subplot(2,1,1);
plot(f,abs(y2(1:512))); %画出滤波前的频谱图 title('滤波前的频谱') xlabel('Hz'); ylabel('fuzhi'); subplot(2,1,2)
F1=plot(f,abs(F0(1:512))); %画出滤波后的频谱图 title('滤波后的频谱') xlabel('Hz'); ylabel('fuzhi');
程序4:窗函数法设计滤波器： fs=22050;
x1=wavread('h:\\课程设计 2\\shuzi.wav');
t=0:1/22050:(size(x1)-1)/22050; Au=0.03;
d=[Au*cos(2*pi*5000*t)]'; x2=x1+d; wp=0.25*pi; ws=0.3*pi; wdelta=ws-wp;
N=ceil(6.6*pi/wdelta); wn=(0.2+0.3)*pi/2;
b=fir1(N,wn/pi,hamming(N+1)); figure(1)
%取整
%选择窗函数，并归一化截止频率
freqz(b,1,512)
f2=filter(bz,az,x2) figure(2)
subplot(2,1,1) plot(t,x2)
title('滤波前的时域波形'); subplot(2,1,2) plot(t,f2);
title('滤波后的时域波形');
sound(f2,22050); %播放滤波后的语音信号 F0=fft(f2,1024); f=fs*(0:511)/1024; figure(3)
y2=fft(x2,1024); subplot(2,1,1);
plot(f,abs(y2(1:512))); title('滤波前的频谱') xlabel('Hz'); ylabel('fuzhi'); subplot(2,1,2)
F2=plot(f,abs(F0(1:512))); title('滤波后的频谱') xlabel('Hz'); ylabel('fuzhi');
1：使用 wavrecord函数 2：使用 awgn等函数
3：使用 fdatool可以设计滤波器 5：wavplay 6: guide matlab 滤波
matlab滤波的方法包括平滑滤波，fir，irr，小波滤波，小波包滤波，自适应
lms，rls滤波，最佳 fir滤波，fadtool，或者使用 filter design工具箱。
filter适合于平滑滤波，对于无特殊要求的滤波。y = filter(b, a, x)，其中b,a为滤波 器系数。计算系统在输入x作用下的零状态响应y[k]。可以结合filter函数和fir滤波或者
irr滤波。b,a可以使fir或者irr产生的滤波器系数。
fir滤波线性相位滤波器适用于比较严格的场合，irr对于一般的滤波，对相位无严格 要求。小波滤波可以实现多频率的，多尺度的信号分解，把信号分为高频和低频，在振动 信号处理中，高频信号一般为噪声。小波包是小波的升级和改进，可以对高频实现更高的 分辨率分解。
1 基于 filter平滑滤波（效果没有 fir，iir好），这个滤波器的指标应该是最小 的均方差。
f1=3;f2=40;fs=100;t=0:1/fs:1;
x=sin(2*pi*t*f1)+.25*sin(2*pi*t*f2);
b=ones(1,10)/10; y=filter(b,1,x); plot(t,x,'r',t,y,'g')。
从图中看出 filter方法可实现平滑滤波，但有相对迟滞。
2 filtfilt 零相位平滑滤波（不推荐，这里 b的长度取值对滤波后信号的幅 值有影响）
f1=3;f2=40;fs=100;t=0:1/fs:1;
x=sin(2*pi*t*f1)+.25*sin(2*pi*t*f2); b=ones(1,10)/10 yy=filtfilt(b,1,x); plot(t,x);hold on; plot(t,x,'r',t,yy,'g')
3 IIR滤波器,结合零相位滤波 filtfilt
fp1=3;fs1=10;f2=40;Fs=100;t=0:1/fs:1;Rp=1;As=30;%此处采样率并不是越 高越好，一般设为信号最高频率 5-15倍即可。 x=sin(2*pi*fp1*t)+.25*sin(2*pi*f2*t); [n,wn]=buttord(2*fp1/Fs,2*fs1/Fs,Rp,As); [b,a]=butter(n,wn);%默认为低通滤波 yf=filter(b,a,x);%用平滑滤波 y=filtfilt(b,a,x) ;
plot(t,x,'r',t,y,'g',t,yf,'b')
4 传统 FIR滤波器
方法一采用窗的精度计算所需滤波器阶数。 fp1=3;fs1=4 ;f2=40;Fs=100;t=0:1/Fs:4; x=sin(2*pi*t*fp1)+.25*sin(2*pi*t*f2) ; wp = fp1*pi/50;ws =fs1*pi/50;
deltaw= ws - wp;% 过渡带宽 Δω的计算
N = ceil(6.6*pi/ deltaw) + 1; % 按海明窗计算所需的滤波器阶数N0 wdham = (hamming(N+1))';% 海明窗计算 Wn=(3+4)/100;% 计算截止频率 b=fir1(N,Wn,wdham); y=filter(b,1,x) ;
yt=filtfilt(b,1,x) ;此函数要求 x长度至少为滤波器阶数的 3倍。 方法2 直接设出滤波器阶数，然后凑试