Scilab -找到Butterworth系数和过滤模拟数据

Question

对于一个带有加速度计的项目，我正在寻找一种过滤一些高频的方法。我不是一个信号专家，但我读到，通常为我的需要，巴特沃斯滤波器是使用。

为此，我使用Scilab，并且我正在与对我的结果的解释进行斗争:我编写了以下代码(底部)并比较了Scilab输出函数，以便在Arduino中实现这个滤波方程：

我的问题是：

• cnum函数和过滤器函数有什么区别？
• 分析结果系数分别为：

系数: Num / Den： 1.

1. 0.063662 0.0020264 0.0000323

其中，我期望一个结果为标准形，其期望值为

b1 0.331785003；b2 0.99535501；b3 0.99535501；b4 0.331785003 a1 -0.965826145；a2 -0.582614466；a3 -0.106171201

从外部软件(https://www.meme.net.au/butterworth.html)计算N=3，低通，Fc=5，Fs=15。

你能看我的代码，并给我一些建议来纠正它，并有正确的系数ai和bi？

/////////////////// cleanup
clear;
//clc;
close;


////////////////////// variable declaration
fcut = 5; //cut off frequency hz (delta 1)
fsampl = 15 ; //sampling frequency hz (delta 2)
delta1_in_dB = -3; // attenuation value at fcut
delta2_in_dB = -21; // final attenuation
delta1 = 10^(delta1_in_dB/20)
delta2 = 10^(delta2_in_dB/20)
epsilon = 0; //ripple value [0 1]
rp = [epsilon epsilon] // ripple vector for analpf 

// conversion of attenuation
delta1 = 10^(delta1_in_dB/20);
delta2 = 10^(delta2_in_dB/20);

// N computation
N = log10((1/(delta2^2))-1)/(2*log10(fspan/fcut));
N = ceil(N);
disp("Order",N);

/////////////////// compute different functions to compaire Butterworth
[poleZP,gainZP]=zpbutt(N,fcut*2*%pi);

[hsAna,poleAna,zeroAna,gainAna]=analpf(N,'butt',rp,fcut*2*%pi);

//disp("Pole : Zpbutt ",poleZP , "Analpf ",poleAna)
//disp("Gain : Zpbutt ",gainZP , "Analpf ",gainAna)
//disp("function :",hsAna)
//       conclusion : zpbutt et analpf donnent la même sortie

/////////////////// paramters in linear system
// Generate the equivalent linear system of the filter
num   = gainAna * real(poly(zeroAna,'s'));
den   = real(poly(poleAna,'s'));
elatf = syslin('c',num,den);

Cnum=coeff(num);
Cden=coeff(den)/Cnum;
Cnum=1;

disp('coefficients : Num / Den : ',Cnum,Cden)





/////////////////// plot an exemple to compare csim and filter
rand('normal');
Input = rand(1,1000); // Produce a random gaussian noise
t     = 1:1000;
t     = t*0.01; // Convert sample index into time steps

y_csim = csim(Input,t,elatf); // Filter the signal with csim
y_res = filter(Cnum, Cden, Input) // Filter the signal with filter


// plot curves
subplot(3,1,1);
plot(t,Input);
xtitle('The gaussian noise','t','y');
subplot(3,1,2);
plot(t,y_csim,'b');
xtitle('The ''csim'' filtered gaussian noise','t','y');
subplot(3,1,3);
plot(t,y_res,'r');
xtitle('The ''filter'' filtered gaussian noise','t','y');

谢谢您的支持！

Answer 1

在将该网站的输出转换为规范形式时，您似乎错过了一个符号反转。

该网站的产出如下：

3.014⋅yi = (1⋅xi + 3⋅xi-1 + 3⋅xi-2 + 1⋅xi-3) + -2.911⋅yi-1 + -1.756⋅yi-2 + -0.320⋅yi-3

当我把它转换成规范形式时，我得到：

yi   (0.331785 + 0.995355*z^-1 + 0.995355*z^-2 + 0.331785*z^-3)
-- = -----------------------------------------------------------
xi      (1 + 0.96582*z^-1 + 0.582614*z^-2 + 0.10617*z^-3)

这意味着：

[b0..b3] = [0.331785   0.99535501 0.99535501 0.331785  ]
[1 a1..a3] = [1.         0.96582614 0.58261447 0.1061712 ]

在scilab中获得相同答案的一个简单方法是使用它们的iir滤波器设计函数，它可以一步完成模拟设计和双线性变换，如下所示：

--> hz=iir(N,'lp','butt',5./15.,[0 0])
 hz  = 

                                                    
   0.3318051 +0.9954154z +0.9954154z² +0.3318051z³  
   -----------------------------------------------  
       0.1060171 +0.5826442z +0.9657797z² +z³       

你可以说系数是向后看的，但是一定要注意到z的幂的符号；整个方程被乘以z/z，相对于你想要的标准形式。处理这个问题的简单方法就是获取系数和把它们从左向右翻转。

--> coeff(hz.num)(:,$:-1:1)
 ans  =

   0.3318051   0.9954154   0.9954154   0.3318051

--> coeff(hz.den)(:,$:-1:1)
 ans  =

   1.   0.9657797   0.5826442   0.1060171

Answer 2

下面是我的更新代码，其中包含来自@Mark的建议解决方案(再次感谢!)我仍然有一些问号，主要是面向Scilab的，如果有人能教我(但这不是我的项目的阻塞点)。

flts和filter函数有什么区别？我的输出结果与那些函数不同。

代码：

////////////////////// cleanup
clear;
//clc; clear console
close;

////////////////////// variable declaration
fcut = 10000; //cut off frequency hz (delta 1)
fsampl = 100000 ; //sampling frequency hz (delta 2)
delta1_in_dB = -3; // attenuation value at fcut
delta2_in_dB = -21; // final attenuation

// conversion of attenuation
delta1 = 10^(delta1_in_dB/20);
delta2 = 10^(delta2_in_dB/20);

// order N computation
N = log10((1/(delta2^2))-1)/(2*log10(fsampl/fcut));
N = ceil(N);
disp("Order",N);

///////////////////// computer Butterworth transfer function
hz=iir(N,'lp','butt',fcut/fsampl,[ ]);

///////////////////////////// extract coefficient from transfer function
num   = coeff(hz(2));
den   = -1*coeff(hz(3));
gain = 1/num(1)

// display value and graph
disp('gain',gain)
disp('''b'' coefficients : num xi',num)
disp('''a'' coefficients : den yi',den(N:-1:1))

[hzm,fr]=frmag(hz,256);
f = scf(0)
plot(fr,hzm)


/////////////////// plot an exemple
rand('normal');
Input = rand(1,1000); // Produce a random gaussian noise
t     = 1:1000;
t     = t*0.01; // Convert sample index into time steps

//Change from transfer function to linear system
    sl= tf2ss(hz)

//Filter the signal
    fs=flts(Input,sl);

// plot curves of raw and filtered data
f = scf(1)
subplot(2,1,1);
plot(t,Input,'b');
xtitle('The gaussian noise','t','y');
subplot(2,1,2);
plot(t,fs,'r');
xtitle('The filtered gaussian noise','t','y');




///////////////////// export on CSV
filename = fullfile(TMPDIR, "data_test_filter.csv");
separator = ";"
M = [t' Input' fs']; //data matrix
csvWrite(M, filename, separator);

filename = fullfile(TMPDIR, "coeff.csv");
C = [gain num den(N:-1:1)] ;//coeff matrix
csvWrite(C, filename, separator);


disp("DONE !");