高通滤波算法

frankneu

加速度计积分得到速度和位移，速度和位移曲线需要进行高通滤波方可消除低频误差，特向大家请教如何实现软件高通滤波，谢谢！！

songzy41

数字滤波器（包括高通）的设计方法有很多，楼主最好先看一下有关数字信号处理的书藉。

frankneu

谢谢赐教！！
我现在需要用加速度传感器测量一个频率在0.2Hz左右的运动信号，由于我所进行的是在线信号处理，先用数据采集卡采进加速度信号后进行数字高通滤波算法处理，

其中用Matlab设计的滤波器的传递函数如下：

0.99539 (s-1)^4
-------------------------------------------------
(s^2 - 1.997s + 0.9973)(s^2 - 1.993s + 0.9935)
算 法语句如下：
w2 = 1.997 * w1 - 0.9973 * w0 + x2 - 2 * x1 + x0
x0 = x1
x1 = x2
v2 = 1.993 * v1 - 0.9935 * v0 + w2 - 2 * w1 + w0
w0 = w1
w1 = w2
y2 = 0.99539 * v2
v0 = v1
v1 = v2
就是不知道这么计算是不是正确的数字滤波方法

songzy41

楼主给出的传递函数用S来表示，不知这传递函数是在S平面（模拟滤波器）还是在Z平面（数字滤波器）。楼主采集到的是数字信号，要处理的话应用数字滤波器。

frankneu

传递函数中s应该是Z,我用Matlab中的数字滤波器设计模块设计的数字高通滤波器的Z传递函数,算法语句是我参考程佩青<数字信号处理>中所介绍的数字滤波器软件实现方法所编写的,由于第一次接触,所以不知道所谓的软件数字滤波就是通过那些语句实现的??

songzy41

按下式传递函数，S为Z，
0.99539 (s-1)^4
-------------------------------------------------
(s^2 - 1.997s + 0.9973)(s^2 - 1.993s + 0.9935)
则要变成Z^(-1)的表达式：
0.99539 (1-2Z^(-1)+Z^(-2))                 (1-2Z^(-1)+Z^(-2))
-----------------------------------------  ×  -----------------------------------------
(1 - 1.997 Z^(-1) + 0.9973Z^(-2))    (1 - 1.993Z^(-1) + 0.9935Z^(-2))
这是2个二阶滤波器的串接，我先检查一下它的频响，程序是：
b0=0.99539;
B=[1 -2 1;1 -2 1];
A=[1 -1.997 0.9973; 1 -1.993 0.9935];
[b, a]=cas2dir(b0, B, A);
freqz(b,a);
从频响来看完全不是高通滤波器，所以请楼主列出设计滤波器系数的程序。

frankneu

我使用Matlab\Simulink\filtering\filter-designs\digital filterdesigil工具箱设计数字高通滤波器，其中参数如下：
desgn method:IIR Butterworth;Fs=10,Fstop=0.001,Fpass=0.01,Astop=90,Apass=1将该滤波器导入到Workspace得到的参数如下：
Den=1        -3.9908        5.9723        -3.9724        0.99081
Num=0.99539        -3.9816        5.9724        -3.9816        0.99539
然后就变换成上面的Z传递函数，其中s用该用Z代替。
我一直不是很明白的就是这种设计和软件实现数字滤波器的做法到底是否正确，敬请指教，谢谢！

songzy41

按照楼主提供的系数，求出了滤波器的频响曲线，程序为：
Den=[1        -3.9908        5.9723        -3.9724        0.99081];
Num=[0.99539        -3.9816        5.9724        -3.9816        0.99539];
freqz(Num,Den);

frankneu

得到滤波器后我的做法如下：
pk(tf(Num,Den))     %求零极点模型

Zero/pole/gain:
0.99539 (s-1.069) (s-0.9351) (s^2  - 1.996s + 1)
------------------------------------------------
  (s-1.1) (s-0.9055) (s^2  - 1.985s + 0.9943)
然后利用滤波器级联结构（程佩青，《数字信号处理》p324）得到上面的滤波器算法语句，不知数字滤波器是不是这样实现的？？

songzy41

如果已知num和den，要求滤波器的零极点，可用tf2zp函数；如果已知num和den而要把它展开成级联形式，可用dir2cas函数。dir2casr函数并不是MATLAB中带有的，但在很多基于MATLAB的信号处理书中都有。

[ 本帖最后由 songzy41 于 2006-12-11 20:19 编辑 ]

pigmin

多谢赐教

frankneu

多些赐教，在下不才，还望继续执教：已知Den，Num，如何利用[b, a]=cas2dir(b0, B, A);怎么才能展开成级联形式呀，是否能提供以下cas2dir函数的源程序呀，谢谢！
用tf2zp函数得到的结果如下，
ans =
   1.1004         
   0.9924 + 0.0971i
   0.9924 - 0.0971i
   0.9055  
也不知道哪个是零点，哪个是极点呀，Matlab学艺不精，见笑见笑了！

songzy41

有Den和Num后，利用tf2zp:
[Z,P,K] = tf2zp(Num,Den)得
Z =
   1.0694         
   0.9978 + 0.0668i
   0.9978 - 0.0668i
   0.9351         
P =
   1.1004         
   0.9924 + 0.0971i
   0.9924 - 0.0971i
   0.9055         
K =
    0.9954

songzy41

我在笫10楼写错了，应用dir2cas函数，源代码为：
function [b0, B, A]=dir2cas(b,a);
%
b0=b(1); b=b/b0;
a0=a(1); a=a/a0;
b0=b0/a0;
%
M=length(b); N=length(a);
if N>M
   b=[b zeros(1,N-M)];
elseif M>N
   a=[a zeros(1,M-N)]; N=M;
end
%
K=floor(N/2); B=zeros(K,3); A=zeros(K,3);
if K*2==N;
   b=[b 0];
   a=[a 0];
end
%
broots=cplxpair(roots(b));
aroots=cplxpair(roots(a));
for i=1:2:2*K
   Brow=broots(i:1:i+1,:);
   Brow=real(poly(Brow));
   B(fix(i+1)/2,:)=Brow;
   Arow=aroots(i:1:i+1,:);
   Arow=real(poly(Arow));
   A(fix(i+1)/2,:)=Arow;
end

frankneu

非常感谢，今天又有学习的东西了。