Pulsar 中Integerator基本功能包含微分、积分与低通滤波,其中低通滤波为一阶低通滤波器,本身滤波性能一般,主要是通过高通与低通的组合,实现带有低通滤波的微分及积分(可高通滤波)功能,高通滤波器本身可等效为微分器,低通滤波器可等效为积分器的功能。
微分Differential (Z domain)的z域表达式
dtdy=Tszz−1
X(z)KHPF+z−1(1−KHPF)Y(z)=Y(z)
⇒X(z)Y(z)=z+KHPF−1zKHPF
在Pulsar中是设定系数KHPF来决定截止频率,因此需要确认出KHPF与τ的关系
H(z)=1+Tsτzz−11
1+Tsτzz−11=z+KHPF−1zKHPF
⇒zKHPF+Tsτ(z−1)KHPF=z+KHPF−1⇒z(1+Tsτ)KHPF−TsτKHPF=z+KHPF−1
⇒⎩⎪⎪⎪⎨⎪⎪⎪⎧(1+Tsτ)KHPF=1−TsτKHPF=KHPF−1
⇒τ=1−KHPFKHPFTs=2πfc1
⇒fc=2π(1−KHPF)TsKHPF
当Ts=1/2000s,KHPF=0.2 时,计算得到fc=79.57754Hz。
Simulink中仿真结果如下:
在Pulsar中设置相同参数,2000Hz下,生成1-200Hz扫频信号,经滤波后结果如下,与仿真结果一致:
X(z)Y(z)=(1−z+KHPF−1zKHPF)=z−(1−KHPF)(z−1)(1−KHPF)
第二个滤波器为了可实现绝对积分的功能,信号走向及滤波方式进行了微调,如下图,直观上当K2=1时,为绝对的积分功能,当K2不为1时,进行等效变换到模型下方的形式,仍然为低通滤波器。
XY=z−Droopz⋅Droop
以截止频率相同的高通低通相连的两个滤波器进行串联后,相位角相加,幅值响应相乘(换成对数或dB后相加),可见在低频阶段f≪51fcut−off以高通滤波器的低频衰减为主(幅度增益20dB/10Oct,相位超前90°,这实际上是典型的微分特征),而在高频部分,f≫5fcut−off以高通滤波器的低频衰减为主(幅度增益−20dB/10Oct,相位滞后90°,这实际上是典型的积分特征)。
当低通滤波器与高通滤波器串联后,
传递函数如下:
X(z)Y(z)=z−(1−KHPF)(z−1)(1−KHPF)⋅z−Droopz⋅Droop
通过第一个低通滤波器实现,参见“First LowPass Filter”,不再赘述。
当第一级低通滤波器截止频率为0Hz时,即KHPF=0所有信号全部通过,同时KDroop=1
系统传递函数为:
X(z)Y(z)=z−(1−KHPF)(z−1)(1−KHPF)⋅z−Droopz⋅Droop∣∣∣∣∣∣KHPF=0KDroop=1
⇒X(z)Y(z)=z−(1−0)(z−1)(1−0)⋅z−1z⋅1=z−1z
对比标准积分的z变换表达式:
∫ydt=Tsz−1z
不考虑单位及Signal Conditioner Scale的前提下,为满足X(z)Y(z)⋅λ=Tsz−1z时,λ=Ts,当Fs=2000Hz时,λ=Ts=Fs1=0.0005,与Integrator.exe计算一致。
考虑到Scale的影响,该系数 IntegratorScale 实际应为
λ=TskUnitScalekOutputScalekInputScale
纯积分的状态下,无法消除静态偏差,即LowPassFilter输出一直为0,静态偏差持续输出至第二个低通滤波器,由于KDroop=1,导致输出会持续累加(积分)输出。【常规情况下不建议单独直接使用】
为了解决当输入信号存在一个静态的偏移而不希望输出持续积分输出的情况下,通过设置一个截止频率相对较低高通的滤波器后再进行积分是比较稳妥的选择。
积分状态下,通常要求KHPF尽可能小但不为零(Integrator.exe要求不能大于0.1,推荐尽可能小),Droop尽可能大(Integrator.exe要求不能小于0.9),在公式中可认为KHPF→0,Droop→1,
X(z)Y(z)=z−(1−KHPF)(z−1)(1−KHPF)⋅z−Droopz⋅Droop∣∣∣∣∣∣KHPF→0KDroop→1≈z−1z
因此此时的Integrator Scale与纯积分时一致:
λ=TskUnitScalekOutputScalekInputScale
以迭代频率2000Hz,截止频率fcut−off=0.12Hz对应设定KHPF=0.00037685,Droop=0.99962为例,输入扫频速率4Oct/min的1-200Hz扫频信号,采集输入与输出,并将其与离线采集的正弦输入直接进行微分进行比较:
subplot(2)中Pulsar与Simulink仿真结果完全一致,与sDAP直接对原始信号进行积分存在细微差异。
subplot(4)中可得到高通部分幅值基本没衰减,相位为0(该频段超过了10倍截止频率)
subplot(4)中可得到第二个低通频率幅值基本在1Hz附近衰减幅度超过了-20dB/10Oct,导致幅值略小于真实的积分值,相位滞后90°,与预期一致,在下图的sDAP分析中也可验证得到相同的结论。
X(z)Y(z)=Tszz−1⋅z−1Tsz⋅z−(1−KHPF)(z−1)(1−KHPF)⋅z−Droopz⋅Droop
X(z)Y(z)=Tszz−1⋅Ts⋅z−(1−KHPF)z(1−KHPF)⋅z−Droopz⋅Droop
X(z)Y(z)=Tszz−1⋅Ts⋅z−(1−KHPF)z⋅KHPF⋅KHPF1−KHPF⋅z−Droopz⋅(1−Droop)⋅1−DroopDroop
X(z)Y(z)=Tszz−1⋅Ts⋅z−(1−KHPF)zKHPF⋅KHPF1−KHPF⋅z−Droopz⋅(1−Droop)⋅1−DroopDroop
当f<51fcutoff时,低通幅值基本不衰减,即z−(1−KHPF)zKHPF≈1,z−Droopz⋅(1−Droop)≈1
X(z)Y(z)=Tszz−1⋅(KHPF1−KHPF⋅1−DroopDroop⋅Ts)
将微分的z变换与Integrator的z变换结果进行比较,假定相差λ的系数:
dtdY=Tszz−1=XYλ=XY/(KHPF1−KHPF⋅1−DroopDroop⋅Ts)
λ=1/(KHPF1−KHPF⋅1−DroopDroop⋅Ts)
λ=(1−KHPFKHPF⋅Droop1−Droop⋅Fs)
当Droop=0.8,UnitInputScale=1mm,UnitOutputScale=1m/s时,IntegratorScaleλ=125,与Integrator Setup.exe计算结果一致。
Pulsar进行微分后,在频率较高后,信号的绝对微分与Integrator的输出存在较大的差异,与预期一致。
截止频率79.57Hz时,Integrator微分与实际信号的真实微分吻合的阶段基本控制在1/10fcut−off,及在8Hz以内。主要是随着频率的升高,频率越接近甚至超过截止频率时,高通滤波器已经无法维持与理想微分器20dB/10Oct的幅度增益,导致存在衰减,但这对于Pulsar中实际应用影响较小,一般Pulsar中对位移进行微分主要有两方面左右
1)位移微分出速度,一般速度段在5Hz以内,该区间结果比较吻合,不受影响;
2)位移Command微分得到速度Command进行前馈,这部分通常高频不需要太多增益,过高的高频增益可能会带来系统的不稳定。
对比Pulsar输出与Simulink仿真结果,两者进行XY绘图时,可见时域完全一致。
Integrator的微分和积分两种情况正是上面的高低通滤波器的串联组合,但是作为积分使用时,将截止频率调低,频率响应尽快进入到积分状态,比如设置fcut−off=0.12Hz时,基本可消除直流分量/静态偏移带来的积分误差,而在1Hz后,完全等效积分功能。而作为微分使用时将截止频率调高,当迭代频率2000Hz时,对应KHPF>0.6时,fcut−off=477Hz,基本可认为100Hz以内完全等效微分的功能,同时在更高频有非常好的低通滤波效果,避免高频噪声微分值过大的负面作用。
using BBK.Numeric;
using NWaves.Filters.Base;
using NWaves.Signals;
using MathNet.Filtering.IIR;
using NWaves.Utils;
var data=new List<double>();
int Fs = 1024;
double khp = 0.4, kdroop = 0.8;
double invertKhp = 1 - khp;
IntegratorSet integratorSet = new(khp, kdroop, Fs, 1, false);
Console.WriteLine($"Freq HP:{integratorSet.FreqHp:f2} Hz,Freq Droop:{integratorSet.FreqDroop:f2} Hz");
var iir = new IirFilter(integratorSet.B, integratorSet.A);
var tf = new TransferFunction(integratorSet.B, integratorSet.A);
var tfHp = new TransferFunction(new double[] { invertKhp, -invertKhp }, new double[] { 1, -invertKhp });
var tfDroop = new TransferFunction(new double[] { kdroop/(kdroop/(1- kdroop)), 0 }, new double[] { 1, -kdroop });
var tfRsp = tf.FrequencyResponse(Fs);
var tfHpRsp = tfHp.FrequencyResponse(Fs);
var tfDroopRsp = tfDroop.FrequencyResponse(Fs);
var magnitudeResponse = tfHp.FrequencyResponse(Fs).Magnitude.Select(x => x).ToArray();
switch (ChannelName)
{
case "1":
data = tfRsp.Magnitude.Select(x => x).ToList();
ChannelName = "iirGain";
break;
case "2":
data = tfRsp.Phase.Select(x => x / Math.PI * 180).ToList();
ChannelName = "iirPhase";
break;
case "3":
data = tfHpRsp.Magnitude.Select(x => x).ToList();
ChannelName = "hpGain";
break;
case "4":
data = tfHpRsp.Phase.Select(x => x / Math.PI * 180).ToList();
ChannelName = "hpPhase";
break;
case "5":
data = tfDroopRsp.Magnitude.Select(x => x).ToList();
ChannelName = "droopGain";
break;
case "6":
data = tfDroopRsp.Phase.Select(x => x / Math.PI * 180).ToList();
ChannelName = "droopPhase";
break;
default:
DiscreteSignal sig1 = new((int)SefInfo.SampleRate, data.ToArray().ToFloats());
data = iir.ApplyTo(sig1).Samples.ToDoubles().ToList();
break;
}