最近初入信号处理，对带通滤波器的资料进行总结。

转载自：


如何在matlab用带通滤波器进行滤波？


matlab中fdatool使用说明


Matlab fir滤波(conv)


以上链接详细记录了带通滤波器的设计及matlab自带fdatool的使用步骤。
自己编的一个带通滤波器的程序及解释
fs=500;                                                 % 采样率
n=length(data);                                         % 滤波数据长度
t = 0:1/fs:(n-1)/fs;                                    % 横坐标刻度值
lpass=100;
hpass=200;                                              % 带通滤波的带宽
filterorder =2;                                         % 滤波器阶数
filtercutoff = [2*lpass/fs 2*hpass/fs];                 % 频率归一化 对应频率/nyquist频率，nyquist频率是采样率的一半 
[f_b, f_a] = butter(filterorder,filtercutoff);          % 调用butter()函数返回巴特沃斯滤波器传递函数的零极点系数和增益
sig= filtfilt(f_b,f_a,data);                            % 零相位数字滤波 通过x在正向和反向两个方向上处理输入数据来执行零相位数字滤波 。在向前过滤数据后，反转过滤后的序列，然后将其运回过滤器。
noise=data-sig;                                         % 噪声为原始信号减去滤波后信号的剩余信号
SNR=10*log((sig.^2)/(noise.^2));                        % 10lg(Ps/Pn)，其中Ps和Pn分别代表信号和噪声的有效功率，也可以换算成电压幅值的比率关系：20Lg(Vs/Vn)，Vs和Vn分别代表信号和噪声电压的“有效值”。用dB表示。

有需要可以自行修改，不保证正确。

    看了很多介绍设计FIR滤波器的，但鲜有告诉你如何应用的。本文以工程师的角度，从介绍、特点、设计使用三个方面出发，并结合代码介绍如何设计并应用FIR滤波器。同时本文也是个人的学习笔记，学习链接也放在了下面，如果不足，请多指导。
介绍：What is FIR filter?   
特点：Why is FIR filter?
如何设计并使用：How to apply a designed FIR filter?
一，介绍：What is FIR filter? 
线性时不变系统（LTI）冲激响应按照其是有限长还是无限长可分为FIR(Finite Impulse Response)有限长冲激响应系统以及无限长冲激响应IIR(Infinite Impulse Response)系统。
关于有限长和无现长的理解，如下图的，该图的冲激响应有无限多个，所以就是无限长冲击响应系统；如果冲激响应是有限个，就是有限长冲击响应系统。
二，特点：Why is FIR filter?
1，传递函数：
2，差分方程
其中N代表滤波器阶数，N越大，该滤波器的幅频响应就会越理想，过渡带就会越陡峭，但缺点是带来了更多的计算量。需要综合考量选择。
FIR是全零点系统，也即在Z复平面上Z传递函数的极点全在Z=0处。
FIR滤波器具有多种实现形式，比如直接型、二阶级联型、Lattice结构，都只是上述基本传递函数的不同数学表达形式，没有本质区别，只是在具体算法实现上各具特点。这里将二阶级联形式描述如下。
二阶级联的意思是将上述传递函数分解为二阶多项式块连乘的形式，其数学表达如下：
援引部分来自手把手教系列之FIR滤波器设计 @逸珺 ，如有侵权，立即删除。
三，如何设计并使用：How to apply a designed FIR filter?
设计方法
FIR滤波器主要设计方法有窗函数法、“最优法”（切比雪夫逼近法、最小均方差）等。其中窗函数法使用最为广泛。“最优法”也比较常用。
最优法”主要思路就是找到一组脉冲响应，让它的频域响应
与期望的滤波器的频域响应
尽可能的一致，主要通过两种方法来实现，一个是最小二乘法，另一个是切比雪夫法。
关于最小二乘法（最小均方差法）和切比雪夫逼近法，以及窗函数方法设计原理和流程，已有大牛介绍的比较好，详见
J Pan：如何快速设计一个FIR滤波器（二）​zhuanlan.zhihu.com
从工程师的角度看这篇文章，虽理论性特别强（有关于连续信号的介绍，个人建议不要太深究。），但是缺乏更贴切的实践和使用介绍。我当时读了几遍之后对实际应用还是有一些疑问。接下来我举个例子，重点从仿真和实际使用来介绍一下。所以这里就不得不提到matlab了。
如何利用MATLAB设计FIR滤波器
如何快速设计一个FIR滤波器（一） 也介绍到，可以通过一种简单设计FIR的方法——零极点法 设计FIR滤波器。
这个方法非常简单，稍加培训，用笔和纸就能完成；当然缺点也很显而易见：零极点设计出的滤波器，只能给出大概的频率响应，对于一些要求较高的系统，显得无能为力。今天我们介绍一种更加严谨的方法。
matlab可以很方便的设计各种滤波器。具体就是命令行输入‘filterDesigner’弹出设计框。如下图，图上方的几个小方框对应着幅频响应、相频响应等。左下方可选择滤波器类型和具体参数等。
举个例子，实现采样频率2kHz，带宽为100Hz～300KHz带通滤波器。
设计一个128阶的FIR带通滤波器，Fstop1为100Hz, Fpass1为110Hz，Fpass2为290Hz, Fstop2为300Hz，Wstop1 为30dB, Wstop2 为30dB。
分析：从下图可以看出，FIR滤波器的相位是线性的。
然后可以拷贝其系数，根据差分方程，进行滤波。
matlab code:
Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = Fs*1;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz , 5.8 , 500 , 120 Hz sinusoid
y = 1*sin(2*pi*50*t) + sin(2*pi*120*t) +  sin(2*pi*5.8*t) +  sin(2*pi*500*t); 
y_target = sin(2*pi*120*t)

N             = 128;     % Order
Fc1           = 100;    % First Cutoff Frequency
Fc2           = 300;    % Second Cutoff Frequency
flag          = 'scale';  % Sampling Flag
SidelobeAtten = 100;      % Window Parameter
% Create the window vector for the design algorithm.
win = chebwin(N+1, SidelobeAtten);

% Calculate the coefficients using the FIR1 function.
b  = fir1(N, [Fc1 Fc2]/(Fs/2), 'bandpass', win, flag);
Hd = dfilt.dffir(b);

figure
freqz(b)

filteredSignal = filter(Hd.Numerator,1,y);
% filteredSignal = filter(b,1,y);

figure
subplot(3,1,1)
plot(t,y)
title('Original Signal')
ys = ylim;
 
subplot(3,1,2)
plot(t,filteredSignal)
title('Target Bandpass Signal')
xlabel('Time (s)'); ylim(ys)

subplot(3,1,3)
plot(t,y_target)
title('Filtered BandPass Signal')
xlabel('Time (s)'); ylim(ys)
C code:  手把手教系列之FIR滤波器设计
其实这部分，对工程师来说很关键啊。
References:
手把手教系列之FIR滤波器设计 
如何快速设计一个FIR滤波器（一） 
如何快速设计一个FIR滤波器（二）

 





 

%本文将针对一个含有 5Hz 、 15Hz 和 30Hz 的混和正弦波信号， 设计一个 FIR 带通滤波器，

%给出利用 MATLAB 实现的三种方法： 程序设计法、 FDATool 设计法和 SPTool 设计法。 参

%数要求：采样频率 fs=100Hz ，通带下限截止频率 fc1=10 Hz ，通带上限截止频率 fc2=20Hz ，过渡带宽 6 Hz，通阻带波动 0.01 ，采用凯塞窗设计。

fc1 =10 ;

fc2 =20 ;

fs=100 ;

[n,Wn,beta,ftype]=kaiserord([7 13 17 23],[0 1 0],[0.01 0.01 0.01],100);

%得出滤波器的阶数 n=38 ， beta=3.4

w1=2*fc1/fs; w2=2*fc2/fs;% 将模拟滤波器的技术指标转换为数字滤波器的技术指

window=kaiser(n+1,beta);% 使用 kaiser 窗函数

b=fir1(n,[w1 w2],window); %使用标准频率响应的加窗设计函数 fir1

freqz(b,1,512);% 数字滤波器频率响应

t = (0:100)/fs;

s = sin(2*pi*t*5)+sin(2*pi*t*15)+sin(2*pi*t*30);% 混和正弦波信号

sf = filter(b,1,s); %对信号 s 进行滤波

 

采用[n,Wn,beta,ftype] = kaiserord(f,a,dev)函数来估计滤波器阶数等，得到凯塞窗滤波器。

　　这里的函数kaiserord(f,a,dev)或者kaiserord(f,a,dev,fs)：

　　f为对应的频率，fs为采样频率；当f用数字频率表示时，fs则不需要写。

　　a=[1 0]为由f指定的各个频带上的幅值向量，一般只有0和1表示；a和f长度关系为（2*a的长度）- 2=（f的长度）

　　devs=[0.05 10^(-2.5)]用于指定各个频带输出滤波器的频率响应与其期望幅值之间的最大输出误差或偏差，长度与a相等。

 

其中fir1(n,Wn)归一化频率的含义

Wn是一个归一化频率，在滤波器设计中是用fs/2进行归一。

在fir1滤波器设计时采用的是归一化频率。实际采样频率为fs，实际的截止频率为fc，设归一化截止频率为fcm，fcm=fc/(fs/2)。

filter函数   参考 https://blog.csdn.net/u012111020/article/details/73744234

b=fir1(n,Wn)返回截止频率为Wn的N阶FIR低通滤波系数行列向量b。

 

 

Matlab 是一个强大的工具，可以用来做各种各样的仿真设计、数字信号处理和科学计算。

由于工作的需要，需要进行数字信号处理，首先就得做仿真，然后将仿真得到的结果再c++上面去实现。

1.先来看看fir滤波器是个啥玩意？

在信号处理领域中，对于信号处理的实时性、快速性的要求越来越高。而在许多信息处理过程中，如对信号的过滤、检测、预测等，都要广泛地用到滤波器。

其中数字滤波器具有稳定性高、精度高、设计灵活、实现方便等许多突出的优点，避免了模拟滤波器所无法克服的电压漂移、温度漂移和噪声等问题，因而随着数字技术的发展，用数字技术实现滤波器的功能越来越受到人们的注意和广泛的应用。其中有限冲激响应(FIR)滤波器能在设计任意幅频特性的同时保证严格的线性相位特性，在语音、数据传输中应用非常广泛。

滤波器分两大类，一种是模拟(ANALOG)滤波器，另外一种是数字(DIGITAL)滤波器。模拟(ANALOG)滤波器是由模拟电路构成，而数字(DIGITAL)滤波器是由数字处理集成电路模块(DSP)和相应的软件构成。

数字(DIGITAL)滤波器是可编程的，所以相对于模拟(ANALOG)滤波器有很多优点。其中最大的优点是通过改变程序或改变程序变量就可设计出不同特点的滤波器，而且数字滤波器可以精确的处理低频率信号。

前面我们提到，什么是滤波器？滤波器就是把噪音去掉，把感兴趣的信号，或者说我们感兴趣的频率信号，从大量信号中提取出来。这如果要用数学语言来表达，那么就是用一个期望的频率特征函数H(f)去乘以输入信号频率X(f)。我们知道，输入信号是有时间性的，它是随着时间的改变而改变。就是说信号是发生在时间空间(时空，TIME DOMAIN)里的，那么，“期望的频率特征函数H(f)去乘以输入信号频率X(f)”这个数学表达在时间空间里是怎样的一个表达式呢？根据傅立叶变换定律，“期望的频率特征函数H(f)去乘以输入信号频率X(f)”在时间空间里就是“这个期望的频率特征函数H(f)在时间空间里的表达式h(t)去和输入信号x(t)做一个卷积”。



具体什么是卷积？用一句经典的话概括：卷积就是各个时刻的输入信号各自乘以相对应的衰减或增幅，然后叠加在一起作为输出信号输出，这里的衰减或增幅就对应与系统的单位冲激响应。——加权叠加。



物理意义：
卷积的重要的物理意义是：一个函数（如：单位响应）在另一个函数（如：输入信号）上的加权叠加。
详细的卷积原理可参考知乎网上的介绍，里面举了非常详细的例子介绍。
卷积讲解可参考：http://blog.csdn.net/bitcarmanlee/article/details/54729807
卷积的应用

用一个模板和一幅图像进行卷积，对于图像上的一个点，让模板的原点和该点重合，然后模板上的点和图像上对应的点相乘，然后各点的积相加，就得到了该点的卷积值。对图像上的每个点都这样处理。由于大多数模板都是对称的，所以模板不旋转。卷积是一种积分运算，用来求两个曲线重叠区域面积。可以看作加权求和，可以用来消除噪声、特征增强。
FIR滤波器是非递归型滤波器的简称，又叫有限长单位冲激响应滤波器。带有常系数的FIR滤波器是一种LTI(线性时不变)数字滤波器。冲激响应是有限的意味着在滤波器中没有发反馈。长度为N的FIR输出对应于输入时间序列x(n)饿关系由一种有限卷积和的形式给出，具体形式如下：

直接形式FIR滤波器图解：

输入信号是有时间性的，随着时间的改变而改变，FIR滤波器最终的输出是各个时刻的输入乘以相应的权重（系数），然后进行叠加，输出。
FIR原理

fir原理
FIR数字滤波器“移动平均数”为例子：

“移动平均数”就是按我们事先设定的信号个数将输入信号加以平均。譬如，如果我们按每4个信号就做一次平均，那么这个4点的“移动平均数”滤波器就如下图所示：



![在这里插入图片描述](https://img-blog.csdnimg.cn/20200409170522920.png?x-oss-process=im
下图是经过11点和51点“移动平均数”滤波器过滤的信号图：

“移动平均数”滤波器的频率响应如下图所示：


如上图所示，随着点数的增加，滚降(ROLLOFF)变陡了，但对旁瓣(sidelobe，衰减部分)的高低影响不大。但是如果我们考虑对滤波器的每个系数采用不同的权重(加权)，而不是像“移动平均数”滤波器那样，用相同的权重(1/4，对4点“移动平均数”滤波器来说)，那么可以期待旁瓣的大小会大大的降低。
对系数采用不同权重的滤波器，我们可以用下面的数学公式来表达：

这就是FIR数字滤波器的一般表达式。
下面我们以设计一个低通滤波器(LPF)为例，来说明FIR数字滤波器窗函数法的设计要点。
假设采样频率为Fs，滤波器的截止(CUT-OFF)频率为Fco，滤波器的长度为Nfir，那么用图形表示出来就如下图所示：

假设Nfir=128，Nco=13 注：Fco=Fs*(Nco/Nfir)，h(t)的计算为：

那么这个低通滤波器的有限冲激响应就如下图所示：

这样我们就设计出了一个FIR低通滤波器。为了检测这个滤波器的性能，我们用信号发生器产生包含从直流到频率为采样频率的一组信号，如下图所示：


我们把这组信号与前面设计的FIR低通滤波器做卷积运算，并将结果(输出)进行快速傅立叶变换(FFT)，得到的频率响应如下图所示：

除了以上方法获得加权系数（抽头系数）外，还可以通过MATLAB获取。
总之，FIR滤波器的目的是滤除不需要的成分，留下需要的成分，如何留下就是通过加权叠加的方式实现。
滤波器除了低通外，还有高通，带通及带阻。
参考文献：http://blog.sina.com.cn/s/blog_74504f8f0100p5ub.html
FIR滤波器设计方法：
直接窗函数设计方法
等波纹设计方法
2.matlab设计64阶FIR带通滤波器

使用matlab自带的滤波器设计工具fdatool，specify order设置比滤波器阶数少一，选择带通滤波器，FIR矩形窗函数。点击designer，即可生成对应的滤波器，选择对应的可以看到相频图，零极点图。





滤波器生成完成后，根据这个滤波器生成相应的系数转化成c语言的头文件，在c平台上即可进行相应的滤波器实现。





c语言代码实现。显示结果，可以明显看到滤波器实现成功，有效的降低了其他带宽，只允许56.8MHZ-96.8MHZ频率通过。