三，滤波器设计：
 
1、相关原理：
 
设计数字滤波器的任务就是寻求一个因果稳定的线性时不变系统，并使系统函数H(z)具有指定的频率特性。
 
数字滤波器从实现的网络结构或者从单位冲激响应分类，可以分成无限长单位冲激响应(IIR)数字滤波器和有限长单位冲激响应(FIR)数字滤波器。
 
数字滤波器频率响应的三个参数：
 
(1) 幅度平方响应：
 
(2) 相位响应
 
其中，相位响应
 
(3) 群时延响应
 
IIR数字滤波器：
 
IIR数字滤波器的系统函数为 的有理分数，即
 
IIR数字滤波器的逼近问题就是求解滤波器的系数 和 ，使得在规定的物理意义上逼近所要求的特性的问题。如果是在s平面上逼近，就得到模拟滤波器，如果是在z平面上逼近，则得到数字滤波器。
 
FIR数字滤波器：
 
设FIR的单位脉冲响应h(n)为实数，长度为N，则其z变换和频率响应分别为
 
按频域采样定理FIR数字滤波器的传输函数H(z)和单位脉冲响应h(n)可由它的N个频域采样值H(k)唯一确定。
 
MATLAB中提供了几个函数，分别用于实现IIR滤波器和FIR滤波器。
 
(1)卷积函数conv
 
卷积函数conv的调用格式为 c=conv(a,b)
 
该格式可以计算两向量a和b的卷积，可以直接用于对有限长信号采用FIR滤波器的滤波。
 
(2)函数filter
 
函数filter的调用格式为 y=filter(b,a,x)
 
该格式采用数字滤波器对数据进行滤波，既可以用于IIR滤波器，也可以用于FIR滤波器。其中向量b和a分别表示系统函数的分子、分母多项式的系数，若a＝1，此时表示FIR滤波器，否则就是IIR滤波器。该函数是利用给出的向量b和a，对x中的数据进行滤波，结果放入向量y。
 
(3)函数fftfilt
 
函数fftfilt的调用格式为 y=fftfilt(b,x)
 
该格式是利用基于FFT的重叠相加法对数据进行滤波，这种频域滤波技术只对FIR滤波器有效。该函数是通过向量b描述的滤波器对x数据进行滤波。
 
关于用butter函数求系统函数分子与分母系数的几种形式。
 
[b,a]=butter(N,wc,'high'):设计N阶高通滤波器，wc为它的3dB边缘频率，以 为单位，故 。
 
[b,a]=butter(N,wc):当wc为具有两个元素的矢量wc=[w1,w2]时，它设计2N阶带通滤波器，3dB通带为 ，w的单位为 。
 
[b,a]=butter(N,wc,'stop'):若wc=[w1,w2]，则它设计2N阶带阻滤波器，3dB通带为 ，w的单位为 。
 
如果在这个函数输入变元的最后，加一个变元“s”，表示设计的是模拟滤波器。这里不作讨论。
 
为了设计任意的选项巴特沃斯滤波器，必须知道阶数N和3dB边缘频率矢量wc。这可以直接利用信号处理工具箱中的buttord函数来计算。如果已知滤波器指标 ， ， 和 ，则调用格式为
 
[N,wc]=buttord(wp,ws,Rp,As)
 
对于不同类型的滤波器，参数wp和ws有一些限制：对于低通滤波器，wpws；对于带通滤波器，wp和ws分别为具有两个元素的矢量，wp=[wp1,wp2]和ws=[ws1,ws2]，并且 ws1
 
2、设计内容：
 
(1)滤波器示例：
 
在这里为了说明如何用MATLAB来实现滤波，特举出一个简单的函数信号滤波实例(对信号x(n)=sin( n/4)+5cos( n/2)进行滤波，信号长度为500点)，从中了解滤波的实现过程。程序如下：
 
Wn=0.2*pi;
 
N=5;
 
[b,a]=butter(N,Wn/pi);
 
n=0:499;
 
x=sin(pi*n/4)+5*cos(pi*n/2);
 
X=fft(x,4096);
 
subplot(221);plot(x);title('滤波前信号的波形');
 
subplot(222);plot(X);title('滤波前信号的频谱');
 
y=filter(b,a,x);
 
Y=fft(y,4096);
 
subplot(223);plot(y);title('滤波后信号的波形');
 
subplot(224);plot(Y);title('滤波后信号的频谱');
 
在这里，是采用了butter命令，设计出一个巴特沃斯低通滤波器，从频谱图中可以很明显的看出来。下面，也就是本课题的主要内容，也都是运用到了butter函数，以便容易的得到系统函数的分子与分母系数，最终以此来实现信号的滤波。
 
(2)N阶高通滤波器的设计(在这里，以5阶为例，其中wc为其3dB边缘频率，以 为单位)，程序设计如下：
 
x=wavread('ding.wav');
 
sound(x);
 
N=5;wc=0.3;
 
[b,a]=butter(N,wc,'high');
 
X=fft(x);
 
subplot(321);plot(x);title('滤波前信号的波形');
 
subplot(322);plot(X);title('滤波前信号的频谱');
 
y=filter(b,a,x);
 
Y=fft(y);
 
subplot(323);plot(y);title('IIR滤波后信号的波形');
 
subplot(324);plot(Y);title('IIR滤波后信号的频谱');
 
z=fftfilt(b,x);
 
Z=fft(z);
 
subplot(325);plot(z);title('FIR滤波后信号的波形');
 
subplot(326);plot(Z);title('FIR滤波后信号的频谱');
 
得到结果如图：
 
(3)N阶低通滤波器的设计(在这里，同样以5阶为例，其中wc为其3dB边缘频率，以 为单位)，程序设计如下：
 
x=wavread('ding.wav');
 
sound(x);
 
N=5;wc=0.3;
 
[b,a]=butter(N,wc);
 
X=fft(x);
 
subplot(321);plot(x);title('滤波前信号的波形');
 
subplot(322);plot(X);title('滤波前信号的频谱');
 
y=filter(b,a,x);
 
Y=fft(y);
 
subplot(323);plot(y);title('IIR滤波后信号的波形');
 
subplot(324);plot(Y);title('IIR滤波后信号的频谱');
 
z=fftfilt(b,x);
 
Z=fft(z);
 
subplot(325);plot(z);title('FIR滤波后信号的波形');
 
subplot(326);plot(Z);title('FIR滤波后信号的频谱');
 
得到结果如图：
 
(4)2N阶带通滤波器的设计(在这里，以10阶为例，其中wc为其3dB边缘频率，以 为单位，wc=[w1,w2],w1 wc w2)，程序设计如下：
 
x=wavread('ding.wav');
 
sound(x);
 
N=5;wc=[0.3,0.6];
 
[b,a]=butter(N,wc);
 
X=fft(x);
 
subplot(321);plot(x);title('滤波前信号的波形');
 
subplot(322);plot(X);title('滤波前信号的频谱');
 
y=filter(b,a,x);
 
Y=fft(y);
 
subplot(323);plot(y);title('IIR滤波后信号的波形');
 
subplot(324);plot(Y);title('IIR滤波后信号的频谱');
 
z=fftfilt(b,x);
 
Z=fft(z);
 
subplot(325);plot(z);title('FIR滤波后信号的波形');
 
subplot(326);plot(Z);title('FIR滤波后信号的频谱');
 
得到结果如图：
 
(5)2N阶带阻滤波器的设计(在这里，以10阶为例，其中wc为其3dB边缘频率，以 为单位，wc=[w1,w2],w1 wc w2)，程序设计如下：
 
x=wavread('ding.wav');
 
sound(x);
 
N=5;wc=[0.2,0.7];
 
[b,a]=butter(N,wc,'stop');
 
X=fft(x);
 
subplot(321);plot(x);title('滤波前信号的波形');
 
subplot(322);plot(X);title('滤波前信号的频谱');
 
y=filter(b,a,x);
 
Y=fft(y);
 
subplot(323);plot(y);title('IIR滤波后信号的波形');
 
subplot(324);plot(Y);title('IIR滤波后信号的频谱');
 
z=fftfilt(b,x);
 
Z=fft(z);
 
subplot(325);plot(z);title('FIR滤波后信号的波形');
 
subplot(326);plot(Z);title('FIR滤波后信号的频谱');
 
得到结果如图：
 
(6)小结：以上几种滤波，我们都可以从信号滤波前后的波形图以及频谱图上看出变化。当然，也可以用sound()函数来播放滤波后的语音，从听觉上直接感受语音信号的变化，但由于人耳听力的限制，有些情况下我们是很难听出异同的。
 
同样，通过函数的调用，也可以将信号的频谱进行“分离观察”，如显出信号的幅值或相位。下面，通过改变系统函数的分子与分母系数比，来观察信号滤波前后的幅值与相位。并且使结果更加明显，使人耳得以很容易的辨听。
 
x=wavread('ding.wav');
 
sound(x);
 
b=100;a=5;
 
y=filter(b,a,x);
 
X=fft(x,4096);
 
subplot(221);plot(x);title('滤波前信号的波形');
 
subplot(222);plot(abs(X));title('滤波前信号的幅值');
 
Y=fft(y,4096);
 
subplot(223);plot(y);title('滤波后信号的波形');
 
subplot(224);plot(abs(Y));title('滤波后信号的幅值');
 
结果如图：
 
>> sound(y);
 
可以听到声音明显变得高亢了。从上面的波形与幅值(即幅频)图，也可看出，滤波后的幅值变成了滤波前的20倍。
 
>> figure,
 
subplot(211);plot(angle(X));title('滤波前信号相位');
 
subplot(212);plot(angle(Y));title('滤波后信号相位');
 
得图：
 
可以看到相位谱没什么变化。
 
(四)、界面设计：
 
直接用M文件编写GUI程序很繁琐，而使用GUIDE设计工具可以大大提高工作效率。GUIDE相当于一个控制面板，从中可以调用各种设计工具以辅助完成界面设计任务，例如控件的创建和布局、控件属性的编辑和菜单设计等。
 
使用GUIDE设计GUI程序的一般步骤如下：
 
1. 将所需控件从控件面板拖拽到GUIDE的设计区域；
 
2. 利用工具条中的工具(或相应的菜单和现场菜单)，快速完成界面布局；
 
3. 设置控件的属性。尤其是tag属性，它是控件在程序内部的唯一标识；
 
4. 如果需要，打开菜单编辑器为界面添加菜单或现场菜单；
 
5. 保存设计。GUIDE默认把GUI程序保存为两个同名文件：一个是.fig文件，用来保存窗体布局和所有控件的界面信息；一个是.m文件，该文件的初始内容是GUIDE自动产生的程序框架，其中包括了各个控件回调函数的定义。该M文件与一般的M文件没有本质区别，但是鉴于它的特殊性，MATALAB把这类文件统称为GUI-M文件。保存完后GUI-M文件自动在编辑调试器中打开以供编辑。
 
6. 为每个回调函数添加代码以实现GUI程序的具体功能。这一步与一般函数文件的编辑调试过程相同。
 
设计过程及内容：
 
在MATLAB版面上，通过键入GUIDE弹出一个菜单栏进入gui制作界面(或者在File到new来进入gui),从而开始应用界面的制作。
 
该界面主要实现了以下几个功能：
 
①打开wav格式的音频文件，并将该音频信号的值读取并赋予某一向量；
 
②播放音频文件，可以选择性的显示该音频信号的波形、频谱、幅值以及相位；
 
③对音频信号进行IIR与FIR的5阶固定滤波处理，可以选择性的显示滤波前后信号的波形、频谱、幅值以及相位，以及播放滤波后的声音。
 
界面如图所示：
 
通过该界面，可以方便用户进行语音信号的处理。
 
界面主程序见附件。
 
(五)、校验：
 
1、本设计圆满的完成了对语音信号的读取与打开，与课题的要求十分相符；
 
2、本设计也较好的完成了对语音信号的频谱分析，通过fft变换，得出了语音信号的频谱图；
 
3、在滤波这一块，课题主要是从巴特沃斯滤波器入手来设计滤波器，也从一方面基本实现了滤波；
 
4、初略的完成了界面的设计，但也存在相当的不足，只是很勉强的达到了打开语音文件、显示已定滤波前后的波形等图。
 
四、 结论：
 
语音信号处理是语音学与数字信号处理技术相结合的交叉学科，课题在这里不讨论语音学，而是将语音当做一种特殊的信号，即一种“复杂向量”来看待。也就是说，课题更多的还是体现了数字信号处理技术。
 
从课题的中心来看，课题是希望将数字信号处理技术应用于某一实际领域，这里就是指对语音的处理。作为存储于计算机中的语音信号，其本身就是离散化了的向量，我们只需将这些离散的量提取出来，就可以对其进行处理了。

http://blog.sina.com.cn/s/blog_574d08530100qu18.html
 

 
 
 
 以下两个滤波器都是切比雪夫I型数字滤波器，不是巴特沃尔滤波器，请使用者注意！
 
 1.低通滤波器
 
 使用说明：将下列代码幅值然后以m文件保存，文件名要与函数名相同，这里函数名：lowp。
 
 function y=lowp(x,f1,f3,rp,rs,Fs)
 %低通滤波
 %使用注意事项：通带或阻带的截止频率的选取范围是不能超过采样率的一半
 %即，f1,f3的值都要小于 Fs/2
 %x:需要带通滤波的序列
 % f 1：通带截止频率
 % f 3：阻带截止频率
 %rp：边带区衰减DB数设置
 %rs：截止区衰减DB数设置
 %FS：序列x的采样频率
 % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值
 % Fs=2000;%采样率
 %
 wp=2*pi*f1/Fs;
 ws=2*pi*f3/Fs;
 % 设计切比雪夫滤波器；
 [n,wn]=cheb1ord(wp/pi,ws/pi,rp,rs);
 [bz1,az1]=cheby1(n,rp,wp/pi);
 %查看设计滤波器的曲线
 [h,w]=freqz(bz1,az1,256,Fs);
 h=20*log10(abs(h));
 figure;plot(w,h);title('所设计滤波器的通带曲线');grid on;
 %
 y=filter(bz1,az1,x);%对序列x滤波后得到的序列y
 end
 
 --------------------------------------
 
 低通滤波器使用例子的代码
 
 fs=2000;
 t=(1:fs)/fs;
 ff1=100;
 ff2=400;
 x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t);
 figure;
 subplot(211);plot(t,x);
 subplot(212);hua_fft(x,fs,1);
 %低通测试
 % y=filter(bz1,az1,x);
 y=lowp(x,300,350,0.1,20,fs);
 figure;
 subplot(211);plot(t,y);
 subplot(212);hua_fft(y,fs,1);%hua_fft()函数是画频谱图的函数，代码在下面给出，要保存为m文件调用
 
 %这段例子还调用了我自己写的专门画频谱图的函数，也给出，不然得不出我的结果
 
 %画信号的幅频谱和功率谱
 %频谱使用matlab例子表示
 function hua_fft(y,fs,style,varargin)
 %当style=1,画幅值谱；当style=2,画功率谱;当style=其他的，那么花幅值谱和功率谱
 %当style=1时，还可以多输入2个可选参数
 %可选输入参数是用来控制需要查看的频率段的
 %第一个是需要查看的频率段起点
 %第二个是需要查看的频率段的终点
 %其他style不具备可选输入参数，如果输入发生位置错误
 nfft= 2^nextpow2(length(y));%找出大于y的个数的最大的2的指数值（自动进算最佳FFT步长nfft）
 %nfft=1024;%人为设置FFT的步长nfft
   y=y-mean(y);%去除直流分量
 y_ft=fft(y,nfft);%对y信号进行DFT，得到频率的幅值分布
 y_p=y_ft.*conj(y_ft)/nfft;%conj()函数是求y函数的共轭复数，实数的共轭复数是他本身。
 y_f=fs*(0:nfft/2-1)/nfft;�T变换后对应的频率的序列
 % y_p=y_ft.*conj(y_ft)/nfft;%conj()函数是求y函数的共轭复数，实数的共轭复数是他本身。
 if style==1
     if nargin==3
         plot(y_f,2*abs(y_ft(1:nfft/2))/length(y));%matlab的帮助里画FFT的方法
         %ylabel('幅值');xlabel('频率');title('信号幅值谱');
         %plot(y_f,abs(y_ft(1:nfft/2)));%论坛上画FFT的方法
     else
         f1=varargin{1};
         fn=varargin{2};
         ni=round(f1 * nfft/fs+1);
         na=round(fn * nfft/fs+1);
         plot(y_f(ni:na),abs(y_ft(ni:na)*2/nfft));
     end
 
 elseif style==2
             plot(y_f,y_p(1:nfft/2));
             %ylabel('功率谱密度');xlabel('频率');title('信号功率谱');
     else
         subplot(211);plot(y_f,2*abs(y_ft(1:nfft/2))/length(y));
         ylabel('幅值');xlabel('频率');title('信号幅值谱');
         subplot(212);plot(y_f,y_p(1:nfft/2));
         ylabel('功率谱密度');xlabel('频率');title('信号功率谱');
 end
 end
 
 下面三幅图分别是滤波前的时频图，滤波器的滤波特性曲线图和滤波后的时频图，通过图可以看出成功留下了100Hz的低频成分而把不要的高频成分去除了。
 
 
 
 
 
 
 
 2.高通滤波器
 
 function y=highp(x,f1,f3,rp,rs,Fs)
 %高通滤波
 %使用注意事项：通带或阻带的截止频率的选取范围是不能超过采样率的一半
 %即，f1,f3的值都要小于 Fs/2
 %x:需要带通滤波的序列
 % f 1：通带截止频率
 % f 2：阻带截止频率
 %rp：边带区衰减DB数设置
 %rs：截止区衰减DB数设置
 %FS：序列x的采样频率
 % rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值
 % Fs=2000;%采样率
 %
 wp=2*pi*f1/Fs;
 ws=2*pi*f3/Fs;
 % 设计切比雪夫滤波器；
 [n,wn]=cheb1ord(wp/pi,ws/pi,rp,rs);
 [bz1,az1]=cheby1(n,rp,wp/pi,'high');
 
 %查看设计滤波器的曲线
 [h,w]=freqz(bz1,az1,256,Fs);
 h=20*log10(abs(h));
 figure;plot(w,h);title('所设计滤波器的通带曲线');grid on;
 y=filter(bz1,az1,x);
 end
 
 下面是高通滤波器的例子
 
 fs=2000;
 t=(1:fs)/fs;
 ff1=100;
 ff2=400;
 x=sin(2*pi*ff1*t)+sin(2*pi*ff2*t);
 figure;
 subplot(211);plot(t,x);
 subplot(212);hua_fft(x,fs,1);
 
 %------高通测试
 z=highp(x,350,300,0.1,20,fs);
 figure;
 subplot(211);plot(t,z);
 subplot(212);hua_fft(z,fs,1);
 
 下面三幅图分别是滤波前的时频图，滤波器的滤波特性曲线图和滤波后的时频图，通过图可以看出成功留下了400Hz的高频成分而把不要的低频成分100Hz去除了。
 
 
 
  
 
 
 
 
 

freqz函数
 此函数可以求取离散系统的频响特性
 画出幅频特性曲线和相频特性曲线
 
unwrap函数介绍
 此函数可以解卷绕 ，可以防止相位跳变
 unwrap(pha,tol,dim)；
 pha 可为数列或矩阵，tol 为判断的阈值，dim 参数可以设置对矩阵的行，还是列进行unwrap操作，
 如果对列每一列操作，那么第三个参数可以省略或填1
 而要对每一行操作，那么第三个参数填2
 
filter函数
 以最简单的 y = filter(b,a,X) 为例
 y = filter(b,a,X) 滤除向量X中的数据，其中b是分子系数向量，a是分母系数向量。如果a(1)不等于1的话，则就利用a(1)标准化滤波器系数，可以利用多项式除法使分母变为1；如果 a(1) 等于0，滤波器返回错误值。
 
hankel矩阵
 
toeplitz函数的功能是生成托普bai利兹（toeplitz）矩阵du。
 托普利兹矩阵的特点是：除第一行、第一列外，其他每个元素都与它左上角的元素相同。

【注】：主要说滤波器的保存和调用。
 一、打开App, 找到filter design;
 二、根据要求设计你的滤波器；
 以巴特沃夫带通滤波器为例：
 （1）如果解决了你的问题给我点个赞吧。
 
 （2）点击DesignFilter, 设计好后，点击File->Export；
 如下图所示设置：（也可以自己探索，下图带通以测试可用）
 
 测试代码：
 
clear all;
close all;
clc;
% 采样频率：1000Hz
% 采样时间：就是画图显示的时间，适用于离线操作。
fs = 1000;
T = 1/fs; 
L = 1500; 
t = (0:L-1)*T;
xt = 0.7*sin(2*pi*2*t) + sin(2*pi*120*t) + 5*sin(2*pi*200*t) + 3*sin(2*pi*495*t);
figure(1); % 图1
plot(fs*t(1:1500),xt(1:1500))
title('原始信号');
xlabel('t (ms)');
ylabel('X(t)');

Y = fft(xt);
Y_Mognitude = fftshift(abs(Y))/L; % 求幅值
f = fs*(-L/2:L/2-1)/L;
figure(2); % 图2
plot(f,Y_Mognitude);
xlabel('f/Hz');
ylabel('原始幅频特性');

Y_Frequency = fftshift(angle(Y));
% 设计滤波器
load('testSaveF.mat', 'bp');
bp_xt = filter(bp,xt);
bp_Y_Mognitude = fftshift(abs(fft(bp_xt)))/L;
bp_Y_Frequency = fftshift(angle(fft(bp_xt)));

figure(5); % 图5
plot(f,bp_xt);
title('滤波后时域信号');
xlabel('t/ms');
ylabel('bp_x(t)');

figure(6); % 图6
plot(f,bp_Y_Mognitude);
title('滤波后幅频特性');
xlabel('Hz');
ylabel('幅值');