跳转到内容
VitaLogos
新建笔记

信号与系统 06：滤波、采集与振动应用推导

学习教程更新于 编辑本文

完整入口 · 上一章：DFT 与频谱标定

本章把已推导的数学模型接回传感器、ADC、数字滤波与振动测量。所有器件参数都要另外核对具体硬件；模型推导不等于某个芯片的实测响应。

1. 滤波为什么同时改变幅值和相位

跳转到“1. 滤波为什么同时改变幅值和相位”

连续/离散 LTI 系统由 hh 描述；频域由 HH 描述。低通、高通、带通和带阻是对频率响应形状的分类，而不是另一套与卷积无关的运算。

类型主要保留例子
低通低频变化RC 电容输出、移动平均
高通高频变化差分 x[n]−x[n−1]x[n]-x[n-1]
带通某个中间频带包络解调前选择共振频带
带阻某个频带之外抑制工频附近分量

差分器：

H(z)=1−z−1⇒H(ejΩ)=2je−jΩ/2sin⁡(Ω/2).H(z)=1-z^{-1} \Rightarrow H(e^{j\Omega})=2j e^{-j\Omega/2}\sin(\Omega/2).

它使直流为零，幅值随频率变化；它不是在全频带都等于理想连续求导 jωj\omega。在低频 Ω≪1\Omega\ll1，1−e−jΩ≈jΩ1-e^{-j\Omega}\approx j\Omega，除以 TsT_s 后才近似连续导数。

2. RC 如何得到可运行的近似递推

跳转到“2. RC 如何得到可运行的近似递推”

连续 RC：

dydt=x−yτ.\frac{dy}{dt}=\frac{x-y}{\tau}.

用前向差分近似导数：

y[n+1]−y[n]Ts≈x[n]−y[n]τ.\frac{y[n+1]-y[n]}{T_s}\approx\frac{x[n]-y[n]}{\tau}.

整理：

y[n+1]=(1−Tsτ)y[n]+Tsτx[n].\boxed{y[n+1]=\left(1-\frac{T_s}{\tau}\right)y[n]+\frac{T_s}{\tau}x[n]}.

这叫前向欧拉近似。其零输入误差每步乘 1−Ts/τ1-T_s/\tau，稳定条件：

∣1−Tsτ∣<1⇒0<Tsτ<2.\left|1-\frac{T_s}{\tau}\right|<1 \Rightarrow 0<\frac{T_s}{\tau}<2.

稳定不等于精确。1<Ts/τ<21<T_s/\tau<2 时会交替过冲，虽然最终趋向稳态，但不符合真实一阶 RC 的单调充电形状。

τ=1\tau=1 s、输入 1 V、y[0]=0y[0]=0、Ts=0.1T_s=0.1 s：

y[1]=0.1,y[2]=0.19,y[3]=0.271.y[1]=0.1,\quad y[2]=0.19,\quad y[3]=0.271.

y[10]=1−0.910≈0.651322y[10]=1-0.9^{10}\approx0.651322 V，而精确 1 s 值 1−e−1≈0.6321211-e^{-1}\approx0.632121 V。

3. 零阶保持条件下的精确离散化

跳转到“3. 零阶保持条件下的精确离散化”

假设区间 nTs≤t<(n+1)TsnT_s\le t<(n+1)T_s 内输入保持为 x[n]x[n]。第二章直接积分求解 RC 的方法给出：

y((n+1)Ts)=x[n]+(y(nTs)−x[n])e−Ts/τ.y((n+1)T_s)=x[n]+(y(nT_s)-x[n])e^{-T_s/\tau}.

定义：

a=e−Ts/τ,y[n]=y(nTs),a=e^{-T_s/\tau},\qquad y[n]=y(nT_s),

得到：

y[n+1]=ay[n]+(1−a)x[n].\boxed{y[n+1]=a y[n]+(1-a)x[n]}.

这是该理想输入保持模型的精确采样关系。将下标改成当前时刻 nn：

y[n]=ay[n−1]+(1−a)x[n−1].y[n]=a y[n-1]+(1-a)x[n-1].

零状态 Z 变换：

Y=az−1Y+(1−a)z−1X,Y=az^{-1}Y+(1-a)z^{-1}X, HZOH(z)=(1−a)z−11−az−1.\boxed{H_{\mathrm{ZOH}}(z)=\frac{(1-a)z^{-1}}{1-az^{-1}}}.

注意输入 x[n]x[n] 施加在接下来的一段区间，输出 y[n+1]y[n+1] 是区间末端电压，因此有一拍因果对齐。常用数字平滑器 y[n]=ay[n−1]+(1−a)x[n]y[n]=ay[n-1]+(1-a)x[n] 用当前点计算，其传递函数少了这个 z−1z^{-1}；不能在不说明时间对齐时把两者完全等同。

double state = 0.0;
double update_rc(double held_input, double a) {
    state = a * state + (1.0 - a) * held_input;
    return state; // 返回下一个采样边界处的电压
}

当 Ts=0.1T_s=0.1 s、τ=1\tau=1 s，a=e−0.1a=e^{-0.1}。递推 10 次后 a10=e−1a^{10}=e^{-1}，精确得到 1−e−11-e^{-1}，可与欧拉结果比较。代码假设参数已由正的 Ts,τT_s,\tau 算好；全局 state 仅代表一个通道，复位和多通道实例需要各自管理状态。

4. LL 点移动平均与群延迟

跳转到“4. LLL 点移动平均与群延迟”
y[n]=1L∑m=0L−1x[n−m].y[n]=\frac1L\sum_{m=0}^{L-1}x[n-m].

零状态传递函数：

H(z)=1L∑m=0L−1z−m=1−z−LL(1−z−1).H(z)=\frac1L\sum_{m=0}^{L-1}z^{-m} =\frac{1-z^{-L}}{L(1-z^{-1})}.

在 z=1z=1 用原来的有限和取极限，直流增益为 1。单位圆上：

H(ejΩ)=1L1−e−jLΩ1−e−jΩ=e−j(L−1)Ω/2sin⁡(LΩ/2)Lsin⁡(Ω/2).\begin{aligned} H(e^{j\Omega}) &=\frac1L\frac{1-e^{-jL\Omega}}{1-e^{-j\Omega}}\\ &=e^{-j(L-1)\Omega/2} \frac{\sin(L\Omega/2)}{L\sin(\Omega/2)}. \end{aligned}

实数幅值因子在零点之间可能改变符号，相位会跳 π\pi；在无零点区间展开相位，群延迟为：

τg=−darg⁡HdΩ=L−12 sample.\tau_g=-\frac{d\arg H}{d\Omega}=\frac{L-1}2\ \text{sample}.

换算秒是 (L−1)/(2fs)(L-1)/(2f_s)。L=5L=5、fs=1000f_s=1000 Hz 时延迟 2 ms。启动缓冲区未填满时存在暂态，滤波延迟不能只用“第一个输出几点出现”来定义。

FIR 可设计成线性相位但不必然全部线性相位；IIR 也不能笼统说“没有确定相位”。实际使用应看 HH 的幅频和相频。

5. 抽取为何必须先抗混叠

跳转到“5. 抽取为何必须先抗混叠”

抽取率为整数 RR，操作为 y[m]=v[mR]y[m]=v[mR]，输出采样率：

fs,out=fs,in/R.f_{s,\mathrm{out}}=f_{s,\mathrm{in}}/R.

它并非仅仅让文件变小。频谱会折叠到更窄频带。用整数筛选恒等式：

1R∑ℓ=0R−1ej2πℓn/R={1,n 是 R 的倍数0,否则.\frac1R\sum_{\ell=0}^{R-1}e^{j2\pi\ell n/R} =\begin{cases}1,&n\text{ 是 }R\text{ 的倍数}\\0,&\text{否则}.\end{cases}

把它插入 DTFT 求和，可得：

Y(ejΩ)=1R∑ℓ=0R−1V(ej(Ω+2πℓ)/R).\boxed{ Y(e^{j\Omega}) =\frac1R\sum_{\ell=0}^{R-1} V\left(e^{j(\Omega+2\pi\ell)/R}\right) }.

这是多个频谱分支相加。为保留基带信号，抽取前低通应抑制高于新 Nyquist 频率的分量。例：输入 2000 Hz、R=4R=4，输出 500 Hz，新 Nyquist 为 250 Hz；300 Hz 若未充分抑制，抽取后会显示为 200 Hz。

插值整数倍 RR 的第一步可在样本间插入 R−1R-1 个零，但这会产生频谱镜像；随后需要重建低通及恰当增益。单纯插零或抽点都不是完整的高质量重采样。

6. 从几何和理解 CIC/Sinc

跳转到“6. 从几何和理解 CIC/Sinc”

为避免与 RC 电阻混淆，本节 RR 只表示抽取率；梳状差分延迟 DD 按输出采样计，级数 PP。令 K=RDK=RD。

一阶未归一化等效高采样率滤波器：

1−z−K1−z−1=1+z−1+⋯+z−(K−1).\frac{1-z^{-K}}{1-z^{-1}} =1+z^{-1}+\cdots+z^{-(K-1)}.

右侧是 KK 点求和，所以直流增益为 KK。级联 PP 次、按输入采样率定义的归一化等效 FIR：

Heq(z)=1KP(1−z−K1−z−1)P.\boxed{ H_{\mathrm{eq}}(z)=\frac1{K^P} \left(\frac{1-z^{-K}}{1-z^{-1}}\right)^P }.

单位圆上：

∣Heq∣=∣sin⁡(KΩ/2)Ksin⁡(Ω/2)∣P.|H_{\mathrm{eq}}| =\left|\frac{\sin(K\Omega/2)}{K\sin(\Omega/2)}\right|^P.

这与移动平均级联的响应相同；随后再抽取 RR 倍。CIC 可以用高采样率积分器、抽取、低采样率梳状器实现等效操作，常规结构见 MathWorks CIC 文档。

非零零点由 sin⁡(KΩ/2)=0\sin(K\Omega/2)=0 得到：

Ω=2πqK,q≢0(modK).\Omega=\frac{2\pi q}{K},\qquad q\not\equiv0\pmod K.

换算输入频率：

fq=qfs,inK=qfs,outD.f_q=\frac{q f_{s,\mathrm{in}}}{K} =\frac{q f_{s,\mathrm{out}}}{D}.

原点用极限取直流增益 1，不是零点。当 D=1D=1 时，零点位于输出采样率的非零整数倍，但零点不等于整个阻带都足够衰减。

等效 FIR 每级延迟 (K−1)/2(K-1)/2 个输入样本：

τg=P(K−1)2fs,in s.\tau_g=\frac{P(K-1)}{2f_{s,\mathrm{in}}}\ \mathrm s.

以上群延迟针对给出的等效 FIR 及其时间对齐。积分器的延迟放置、抽取相位、流水寄存器和接口缓冲会增加或改变实际输出时刻；比较器件时序或不同实现时，应另列这些延迟，不能只凭幅频曲线相同就认定输出逐点对齐。

未归一化直流增益 KPK^P，粗略最大位宽增长需考虑 ⌈Plog⁡2K⌉\lceil P\log_2K\rceil；实际定点实现还要规定回绕、截位、舍入和饱和方式。

“sinc 形状”通常指这种正弦比。在 Ω\Omega 足够小时，用 sin⁡(Ω/2)≈Ω/2\sin(\Omega/2)\approx\Omega/2 才接近连续 sinc，不能把两者在全频带严格等同。“Sinc5”可指五级形状，但不独立决定具体 ADC 的数据率或稳定时间。

7. 振动的位移、速度、加速度

跳转到“7. 振动的位移、速度、加速度”

令位移：

d(t)=Adcos⁡(ωt+ϕ),Ad 单位 m.d(t)=A_d\cos(\omega t+\phi),\qquad A_d\ \text{单位 m}.

逐次求导：

v(t)=d′(t)=−Adωsin⁡(ωt+ϕ)=Adωcos⁡(ωt+ϕ+π/2),\begin{aligned} v(t)&=d'(t)=-A_d\omega\sin(\omega t+\phi)\\ &=A_d\omega\cos(\omega t+\phi+\pi/2), \end{aligned} a(t)=v′(t)=−Adω2cos⁡(ωt+ϕ)=Adω2cos⁡(ωt+ϕ+π).\begin{aligned} a(t)&=v'(t)=-A_d\omega^2\cos(\omega t+\phi)\\ &=A_d\omega^2\cos(\omega t+\phi+\pi). \end{aligned}

所以峰值关系：

Av=ωAd,Aa=ω2Ad.A_v=\omega A_d,\qquad A_a=\omega^2A_d.

速度相对位移超前 90∘90^\circ，加速度相对位移差 180∘180^\circ；必须始终沿用相同余弦与参考方向。

例：Ad=10 μm=10−5A_d=10\,\mu\mathrm m=10^{-5} m、f=100f=100 Hz：

ω=2π100≈628.319 rad/s,\omega=2\pi100\approx628.319\ \mathrm{rad/s}, Av≈0.00628319 m/s=6.28319 mm/s,A_v\approx0.00628319\ \mathrm{m/s} =6.28319\ \mathrm{mm/s}, Aa≈3.94784 m/s2.A_a\approx3.94784\ \mathrm{m/s^2}.

这些是峰值；零均值正弦的 RMS 再除 2\sqrt2。频域积分为 V=A/(jω)V=A/(j\omega)、D=A/(jω)2D=A/(j\omega)^2，只在相应频率非零且条件成立时使用。DC 处除法不存在，低频噪声经过积分会被强烈放大，漂移需要单独处理。

8. IQ 解调怎样读出幅度与相位

跳转到“8. IQ 解调怎样读出幅度与相位”

令输入窄带正弦：

x(t)=Acos⁡(ωct+ϕ).x(t)=A\cos(\omega_ct+\phi).

选参考 2cos⁡ωct2\cos\omega_ct 与 −2sin⁡ωct-2\sin\omega_ct。乘法恒等式给出：

2x(t)cos⁡ωct=Acos⁡ϕ+Acos⁡(2ωct+ϕ),−2x(t)sin⁡ωct=Asin⁡ϕ−Asin⁡(2ωct+ϕ).\begin{aligned} 2x(t)\cos\omega_ct &=A\cos\phi+A\cos(2\omega_ct+\phi),\\ -2x(t)\sin\omega_ct &=A\sin\phi-A\sin(2\omega_ct+\phi). \end{aligned}

低通滤掉 2ωc2\omega_c 附近项，保留基带：

I=Acos⁡ϕ,Q=Asin⁡ϕ.I=A\cos\phi,\qquad Q=A\sin\phi.

因此：

A=I2+Q2,ϕ=atan2⁡(Q,I),\boxed{A=\sqrt{I^2+Q^2},\qquad\phi=\operatorname{atan2}(Q,I)}, I+jQ=Aejϕ.I+jQ=Ae^{j\phi}.

本节参考包含系数 2，所以幅度直接为 AA；若不用 2，基带幅度为 A/2A/2。若 Q 用正正弦，符号也会变，不能跨约定直接比较相位。

载波频率、采样时钟和低通带宽需合适；频率不匹配时，相位会随时间旋转。调制的基本频移关系可对照 MIT 连续时间调制讲次。

9. 包络与解析信号

跳转到“9. 包络与解析信号”

对窄带幅度调制 x(t)=A(t)cos⁡(ωct+ϕ)x(t)=A(t)\cos(\omega_ct+\phi)，当 A(t)A(t) 相对载波变化慢且频带适当分离，可用正交分量形成包络。

Hilbert 变换的频率响应为：

HH(jω)=−j sgn⁡(ω).H_{\mathrm H}(j\omega)=-j\,\operatorname{sgn}(\omega).

它使正频率相位移 −90∘-90^\circ、负频率相位移 +90∘+90^\circ；cos⁡\cos 变成 sin⁡\sin。解析信号：

xa(t)=x(t)+jx^(t),x_{\mathrm a}(t)=x(t)+j\widehat x(t), aenv(t)=∣xa(t)∣=x(t)2+x^(t)2.a_{\mathrm env}(t)=|x_{\mathrm a}(t)| =\sqrt{x(t)^2+\widehat x(t)^2}.

对单一正弦 x^=Asin⁡(ωct+ϕ)\widehat x=A\sin(\omega_ct+\phi)，因此包络恒为 ∣A∣|A|。多分量宽带信号的模仍可计算，但不一定有简单的物理“调制幅度”意义。解析信号与包络的计算定义见 SciPy Hilbert 文档。

机械冲击诊断常先带通选共振频带，再取包络并分析包络谱；包络峰值不能独立证明某类故障，需转速、结构与其他证据。详细方法见包络专题。

10. 倒谱为什么能突出重复间隔

跳转到“10. 倒谱为什么能突出重复间隔”

实倒谱的一种定义是：

c[q]=IDFT⁡{ln⁡(∣X[k]∣+ε)},ε>0.c[q]=\operatorname{IDFT}\{\ln(|X[k]|+\varepsilon)\},\qquad\varepsilon>0.

qq 是“倒频率”下标，对应采样间隔可用 q/fsq/f_s 表示。ε\varepsilon 防止对零取对数，并影响弱谱解释；实际还须规定是否采用幅度、功率或复对数。

实倒谱定义可对照 MathWorks rceps 文档。若时域 y=h∗xy=h*x，则频域 ∣Y∣=∣H∣∣X∣|Y|=|H||X|，忽略零点和正则化影响时：

ln⁡∣Y∣=ln⁡∣H∣+ln⁡∣X∣.\ln|Y|=\ln|H|+\ln|X|.

频域的乘法被对数变成加法，逆变换后不同结构有机会在倒谱中分离。

例如回声模型 y(t)=x(t)+αx(t−τd)y(t)=x(t)+\alpha x(t-\tau_d)：

Y=X(1+αe−jωτd).Y=X(1+\alpha e^{-j\omega\tau_d}).

∣α∣<1|\alpha|<1 时，复对数级数：

ln⁡(1+αe−jωτd)=∑m=1∞(−1)m+1mαme−jmωτd.\ln(1+\alpha e^{-j\omega\tau_d}) =\sum_{m=1}^{\infty}\frac{(-1)^{m+1}}m \alpha^m e^{-jm\omega\tau_d}.

这些指数对应延迟 mτdm\tau_d，说明倒谱为何能在重复延迟附近出现结构。实倒谱取对数幅值后还有相应对称项。周期谱线间距 Δf\Delta f 常关联倒谱位置约 1/Δf1/\Delta f，但并非所有倒谱峰都能直接解释为真实回声。工程细节见倒谱专题。

11. 阶次、转速与角域重采样

跳转到“11. 阶次、转速与角域重采样”

转速 rr 单位 rpm，转频：

frot=r60 Hz.f_{\mathrm{rot}}=\frac r{60}\ \mathrm{Hz}.

阶次 oo 是信号频率与转频的比：

o=ffrot,f=or60.\boxed{o=\frac f{f_{\mathrm{rot}}},\qquad f=o\frac r{60}}.

3000 rpm 对应 50 Hz，一阶为 50 Hz、二阶为 100 Hz。转速变化时固定阶次的 Hz 也变化，所以长时间固定时间轴 FFT 可能将峰展宽。

定义轴转角：

θ(t)=θ(0)+2π∫0tfrot(ξ) dξ.\theta(t)=\theta(0)+2\pi\int_0^t f_{\mathrm{rot}}(\xi)\,d\xi.

若单个阶次分量为 Acos⁡(oθ(t)+ϕ)A\cos(o\theta(t)+\phi)，在等间隔转角 θm=mΔθ\theta_m=m\Delta\theta 处重采样：

xθ[m]=Acos⁡(omΔθ+ϕ).x_\theta[m]=A\cos(om\Delta\theta+\phi).

旋转速度的变化已被角域坐标吸收，该模式在角域保持固定阶次。重采样需可靠转速/键相信号、插值与抗混叠；不能仅用一个平均转速除整段频率就称为完整阶次跟踪。进一步见阶次与赫兹。

12. 从传感器到测量结果的单位链

跳转到“12. 从传感器到测量结果的单位链”

若传感器灵敏度为 SS V/(m/s²)，模拟前端增益为 GG，理想 ADC 区间宽度为 Δ\Delta V/code，去偏移码值为 c[n]−c0c[n]-c_0，简化换算：

vADC[n]≈Δ(c[n]−c0),v_{\mathrm{ADC}}[n]\approx\Delta(c[n]-c_0), a[n]≈Δ(c[n]−c0)GS.\boxed{a[n]\approx\frac{\Delta(c[n]-c_0)}{GS}}.

单位检查：V 除以 V/(m/s²)，得到 m/s²。真实转换需具体 ADC 码制、校准偏置和增益；这里仅说明单位链。

对于无混叠、工作频带内的单频成分，模拟频响、数字滤波频响和时间偏差可以分段计算：

Ameas=Atrue∣Hsensor(jω)∣∣Hanalog(jω)∣∣Hdigital(ejΩ)∣,A_{\mathrm{meas}} =A_{\mathrm{true}}|H_{\mathrm{sensor}}(j\omega)| |H_{\mathrm{analog}}(j\omega)| |H_{\mathrm{digital}}(e^{j\Omega})|, ϕmeas=ϕtrue+arg⁡Hsensor+arg⁡Hanalog+arg⁡Hdigital−ωΔt.\phi_{\mathrm{meas}} =\phi_{\mathrm{true}}+\arg H_{\mathrm{sensor}} +\arg H_{\mathrm{analog}}+\arg H_{\mathrm{digital}} -\omega\Delta t.

这里 Ω=ω/fs\Omega=\omega/f_s，Δt\Delta t 是本通道相对参考的时间延迟。传感器方向、负增益、不同延迟与相位约定也要统一。

整个包含 ADC 与抽取的多速率系统不能不加条件地塞进单一连续时间 H(s)H(s)；要按每段采样率和物理连接建模。频响接近零的频带做逆校正会放大噪声，不能无限补偿。

13. 最小工程核对例子

跳转到“13. 最小工程核对例子”

设已换算为加速度的理想信号为 100 Hz、峰值 1 m/s²；两个通道相差 10 μs：

Δϕ=−2π100×10−5=−0.00628319 rad=−0.36∘.\Delta\phi=-2\pi100\times10^{-5} =-0.00628319\ \mathrm{rad}=-0.36^\circ.

若数字两点平均器工作在 fs=1000f_s=1000 Hz：

Ω=2π100/1000=0.2π,\Omega=2\pi100/1000=0.2\pi, ∣H∣=cos⁡(0.1π)≈0.951057,arg⁡H=−0.1π=−18∘.|H|=\cos(0.1\pi)\approx0.951057,\qquad \arg H=-0.1\pi=-18^\circ.

所以滤波后的稳态峰值约 0.951057 m/s²，相位还加上滤波的 −18∘-18^\circ。通道不同步的 −0.36∘-0.36^\circ 与滤波相位是两项不同误差来源。

14. 本章练习与答案

跳转到“14. 本章练习与答案”
  1. τ=1\tau=1 s、Ts=0.1T_s=0.1 s 的精确 RC 系数：a=e−0.1≈0.904837a=e^{-0.1}\approx0.904837。
  2. 5 点平均器、2000 Hz 采样的群延迟：(5−1)/(2×2000)=1(5-1)/(2\times2000)=1 ms。
  3. 输入 4000 Hz 抽取 8 倍：输出 500 Hz，新 Nyquist 为 250 Hz，先检查抗混叠。
  4. 位移 1 mm、频率 10 Hz 的速度峰值：2π10×0.001≈0.06283192\pi10\times0.001\approx0.0628319 m/s。
  5. I=3I=3、Q=4Q=4 的 IQ 幅度：5，角度约 53.130∘53.130^\circ，单位沿用 IQ 输入单位。
  6. 1200 rpm 的三阶频率：3×1200/60=603\times1200/60=60 Hz。
  7. 100 Hz 信号延迟 1 ms：−36∘-36^\circ；不同通道比较前统一参考方向。