赞
踩
前段时间在对用MATLAB处理试验数据时,需要对多路信号进行滤波后做同步,在这个过程中使用MATLAB中自带的 fir1 函数以及Filter Designer工具箱(FDATool)进行FIR滤波器设计时,信号经滤波后会向后时延一定点数,造成同步的不一致。
翻阅本科《数字信号处理》课程相关资料后,发现FIR滤波器具有线性相位的特点,主要表现为:信号通过滤波器前后波形基本不变,但信号整体会产生一个时延 ,也就是滤波器的群时延。
同时,用MATLAB的 filter 函数对输入信号进行滤波时,如果将原信号直接作为输入信号,那么由于FIR滤波器这个时延的特征,输出的信号会整体向后移动长度,移动出去的这小段并没有被保留到输出信号向量中,而是直接溢出了,取而代之的是在信号向量的前面补充了个0,这就导致了输出信号和输入信号无法对齐,且输出信号末端溢出了一部分,这主要是由于filter这个函数的性质造成的。
经过研究我找到了解决方法,因此本文将结合MATLAB中fir1、filter函数以及FIR滤波器的相关知识,来讲解一下如何使FIR滤波前后输入信号和输出信号对齐,且不产生损失。
FIR滤波器的线性相位条件:
设FIR滤波器的长度为N,若FIR滤波器的单位冲激响应 h(n) 满足:
即关于 (N-1)/2 偶对称,则满足第一类线性相位:
若FIR滤波器的单位冲激响应 h(n) 满足:
即关于 (N-1)/2 奇对称,则满足第二类线性相位:
FIR滤波器分类:
FIR滤波器根据 (N-1)/2 点处的对称性质可分为第一类和第二类线性相位,这其中还存在一个问题:N是奇数还是偶数?实际上,N取奇数或偶数并不影响滤波器的线性相位特征,但会使滤波器产生不同的零点,影响滤波器的性能,关于零点的计算这里不具体介绍了,感兴趣的可以去翻数字信号处理的书,只给结论——基于以上性质,FIR滤波器可分为4类:
因此,本文主要讨论是具有第一类线性相位的FIR滤波器,也就是类型Ⅰ和类型Ⅱ,其中,先放一个结论:只建议使用类型Ⅰ,不建议使用其他类型的FIR滤波器,因为这会带来很多问题。
我们用MATLAB的 fir1 函数进行滤波器设计,当然,你也可以选用 Filter Designer 工具箱进行设计,fir1 中的 1 指的并不是指第一类线性相位,也不是指Ⅰ型FIR滤波器,而是利用窗函数设计法设计FIR滤波器的函数,这一点可以参考fir1函数的help文档,写得很清楚。另外,MATLAB还有 fir2 函数,是利用频率采样法设计任意幅频特性的FIR滤波器。本文只讨论日常使用较多的 fir1 函数。
事实上,经过实际使用测试,发现 fir1 函数设计出的滤波器全部都是第一类线性相位的(h(n) 关于 (N-1)/2 偶对称),具体是Ⅰ型还是Ⅱ型,要看滤波器长度N,N为奇数为Ⅰ型,N为偶数则为Ⅱ型。
根据前面介绍的内容,如果你需要设计高通、带阻滤波器,则滤波器长度N必须取为奇数(也就是说这种情况只能选择Ⅰ型),否则,MATLAB会报警告,如图,并自动帮你加1点,以保证选用的是Ⅰ型滤波器,因为Ⅱ型FIR滤波器在高频处有零点,不能用于高通、带阻滤波器的设计。
利用 FIR1 生成一个截止频率为FH的64阶FIR滤波器的代码:
- FIR_n = 64; % FIR滤波器的阶数,(差分方程中,x[n-1],x[n-2]……
- % 这些项的数量为FIR滤波器的阶数FIR_n,不包含x[n]和y[n]
- % 即对应的冲激响应h[n]中,δ[n-1],δ[n-2]这些项的数量为FIR_n,不包含δ[n],
- % 所以,对于FIR滤波器,冲激响应函数序列长度应该是:N = FIR_n+1点)
- % 若要设计Ⅰ型FIR滤波器,则滤波器冲激响应的长度N必须为奇数,那么,
- % 滤波器阶数FIR_n就必须为偶数(= N-1),否则MATLAB会自动帮你补为偶数
-
- Wn = fH * 2/fs; % 将滤波器截至频率转为数字角频率,作为传递进fir1函数的参数
-
- b = fir1(FIR_n,Wn,'low'); % 利用fir1函数生成一个低通FIR滤波器
- % b 可以理解为该滤波器的冲激响应序列 由δ[n] δ[n-1] δ[n-2]……组成
- % 对应滤波器差分方程中x[n] x[n-1]…等各阶项前面的系数
由于FIR滤波器的输出只取决于输入,因此,滤波器长度N = 滤波器阶数+1,具体可阅读上方代码的注释部分。利用 fir1 函数得到的序列b可以认为是滤波器的冲激响应,其长度就是N。
实际上,b也是滤波器有理传递函数分子部分的系数,b序列将在后续作为参数传递给 filter 函数,才能对信号进行滤波。
时延点数
在前一节,我们已经推出FIR滤波器的频域相位函数
可见,FIR滤波器都会使输入信号向后时延 (N-1)/2 个点,表现如下:
(本文中所演示的信号均在所设计FIR滤波器的通带内,可以顺利通过滤波器)
根据此,我们可以对滤波后的产生时延予以消除,这将在下面的内容中结合程序介绍。
另外值得注意的是,当N为偶数时,由于 (N-1)/2 不是整数,这种情况一般可认为时延了 N/2 个点,这一情况在最后讨论。
使用MATLAB的 filter 函数进行滤波,以及输出信号溢出的问题:
前面的内容中,在设计好FIR滤波器后,得到序列b,查阅 filter 函数的help文档可知,b 是其有理传递函数分子部分的系数,filter 函数所对应的滤波器结构可以在help文档中的“More About(详细信息)”部分找到,由于FIR滤波器的输出只取决于输入,结合文档中给出的滤波器结构图,FIR滤波器的有理传递函数分母部分的系数应为1。
因此,接着上文中的程序,之后滤波操作如下:
x_filter = filter(b,1,x); % 利用设计好的滤波器 b 对 x 序列进行滤波
然而,filter 函数有一个性质:输出的 x_filter 序列长度与输入的 x 序列长度是保持一致的,由于FIR滤波器对输入信号产生的时延, 输入信号会被向后移动 (N-1)/2 点,因此这 (N-1)/2 点会被移动到 x 序列长度之外,即溢出了,演示如下:
红框所示的部分就是因为 filter 函数输入输出序列长度保持一致以及FIR滤波器的时延所造成的信号溢出问题,这在信号处理当中往往是不可忽视的,因此,需要对 filter 函数的输入信号序列 x 和其输出序列 x_fliter 进行调整,来避免这一溢出问题,滤波后得到完整信号。
由于信号序列 x 经过 filter 滤波后会因FIR滤波器时延特性而向后溢出一部分,且溢出点数为 (N-1)/2,因此,我们可以在信号序列 x 进入滤波器之前,先在其末尾加入 (N-1)/2 个0点,然后将其送入 filter 函数进行滤波,得到的输出恰好将末尾这(N-1)/2 个0点移了出去,保留了完整信号,同时,由于信号向后时延,会在信号首端多出(N-1)/2 个0点,为了满足输出与输入序列长度的一致性,需要将输出序列前面的(N-1)/2 个0点也删除,最终得到和输入信号序列等长度且对齐的输出信号序列。演示如下,其中红框部分是有用信号:
代码实现如下:
- N = length(b); % 对于FIR滤波器,长度 N = FIR_n+1
- x_filter = filter(b,1,[x zeros(1,(N-1)/2)]); % 对输入序列 x 末端补 (N-1)/2 个 0
- x_filter = x_filter((N-1)/2+1:end); % 去除输出序列中前面的 0
其中的 (N-1)/2 可以也写成 fix(N/2),以应使得长度N为偶数的FIR滤波器也适用。
至此,我们分析并解决了MATLAB中FIR滤波器对输入信号的时延及其溢出问题,使输入信号和输出信号的长度得以对应上,且保持对齐。
然而,还有一个不可回避的问题:当滤波器长度N为偶数时(即Ⅱ型FIR滤波器),时延点数 (N-1)/2 不再是一个整数,那么信号到底向后时延了多少点呢?
这里我认为,取 N/2 或 (N-2)/2 个点均可,因为经过研究发现:Ⅱ型FIR滤波器的效果相当于对输入信号进行了“重采样”,具体表现为输入信号序列上的某一点对应的是输出信号序列上的左右两点之间,各采样点之间无法对应上,而是相差了半个点,如下图所示:
而Ⅰ型FIR滤波器便不存在这个问题。所以,对于Ⅱ型FIR滤波器,实际物理上好像真的把信号时延了 (N-1)/2 个单位长度,只是对于离散的信号,我们只能取整数点来补偿时延,由于 (N-1)/2 不是整数,因此可以取 N/2 或 (N-2)/2 个点作为时延点数进行补偿,具体如何选择要看你所面对的实际问题,如果需要多路严格同步,那么每路都要取一样的点数 ,保持一致。
此外,由于这个“重采样”效果,导致无论你取 N/2 或 (N-2)/2 个点作为时延点数,都无法和原信号的采样点对应上,倘若你后续要做相干解调等对波形相干性要求比较高的操作,那么非常不建议你使用Ⅱ型FIR滤波器。
综上,任何情况下都不建议使用Ⅱ型FIR滤波器,因为问题太多。相比之下,Ⅰ型FIR可以满足各种滤波器的设计需要,且不会出现上述问题,所以还是尽量使用Ⅰ型FIR滤波器吧。
第一次写博文,介绍了MATLAB中FIR滤波器使用的一些问题,笔者刚读研究生,科研能力欠佳,缺乏实际经验,学的也不够深入,如果有错误和问题欢迎第一时间在评论区指出,我会对文章进行修正。
我本人是研究一些信号处理算法的,通信和探测都有涉猎,此后我会不定期分享一些科研过程中遇到的问题,感谢各位老板阅读!
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。