对于带有加速度计的项目,我正在寻找一种过滤一些高频的方法。我不是信号专家,但我通常会根据需要阅读,使用巴特沃斯滤波器。
我为此目的使用 Scilab,我正在努力解释我的结果:我编写了以下代码(在底部)并比较了 Scilab 输出函数以在 Arduino中实现这个过滤方程:
我的问题是:
- cnum 函数和 filter 函数有什么区别?
- analpf(butt) 或 zpbutt 的结果系数为
“系数:Num / Den:”
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
从这里的外部软件(链接)计算,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');
提前感谢您的支持!!