MATLAB 动平衡计算与同步测量示例

这个例子验证什么

本页给出五个独立 .m 文件，放在同一目录后运行 balancingDemo。它们验证单面/双面影响系数计算，以及从已校验的每转参考中提取 1X 复数。演示使用明确单位的合成数据，不连接采集卡、不控制电机，也不把计算残差当作实机验收。

代码仅用 MATLAB 基础函数。2026-10-03 在 MATLAB R2026a Update 2 运行；正文代码由测试直接提取，测试范围与结果见文末。数值输入校验后转为 double，避免整数 ADC 码值在除法与阈值插值中被取整；ADC 码值到物理量的校准仍由调用方完成。理论和相位约定先看 影响系数法 与 键相测量。

1. 单平面求解：singlePlaneBalancing.m

输入振动为复数，质量为 g、半径为 mm、角度为度。返回的 correction 是 g·mm 总校正量；保留试重时使用 additionalIfTrialKept，再除以准备安装的半径得到质量。函数拒绝零试重响应；实际应用还须按重复测量散布判断响应是否足够可信。

function result = singlePlaneBalancing(V0, Vtrial, trialMass, trialRadius, trialDeg, correctionRadius)
% Inputs: complex vibration in one consistent unit; mass g; radius mm; deg.
validateattributes(V0, {'numeric'}, {'scalar','finite'});
validateattributes(Vtrial, {'numeric'}, {'scalar','finite'});
validateattributes(trialMass, {'numeric'}, {'real','scalar','finite','positive'});
validateattributes(trialRadius, {'numeric'}, {'real','scalar','finite','positive'});
validateattributes(correctionRadius, {'numeric'}, {'real','scalar','finite','positive'});
validateattributes(trialDeg, {'numeric'}, {'real','scalar','finite'});
V0=double(V0); Vtrial=double(Vtrial);
trialMass=double(trialMass); trialRadius=double(trialRadius);
trialDeg=double(trialDeg); correctionRadius=double(correctionRadius);
Ut = trialMass * trialRadius * exp(1j * deg2rad(trialDeg));
delta = Vtrial - V0;
if abs(delta) == 0
    error('balancing:NoResponse', 'Trial response must be nonzero.');
end
alpha = delta / Ut;
Uc = -V0 / alpha;
result = struct('alpha',alpha, 'trialUnbalance',Ut, 'correction',Uc, ...
    'mass',abs(Uc)/correctionRadius, 'degrees',mod(rad2deg(angle(Uc)),360), ...
    'additionalIfTrialKept',Uc-Ut, 'predictedResidual',V0+alpha*Uc);
end

零原始振动会得到零校正量，此时安装角没有物理意义。近零试重响应虽可能通过浮点非零检查，仍可能得出极大的不可靠校正量；机械上限和信号质量必须在应用层单独检查。

2. 双面求解：twoPlaneBalancing.m

Vtrials(:,1) 是仅加 A 面试重时两个测点的结果，第二列是移除 A 试重后仅加 B 面试重的结果。Utrials 已包含各面质量、半径及角度，因此可以使用不同试重和半径。

function result = twoPlaneBalancing(V0, Vtrials, Utrials, correctionRadii)
% V0: 2x1; Vtrials: 2x2, each column measured with one independent trial.
% Utrials: complex 2x1 in g*mm; correctionRadii: 2x1 in mm.
validateattributes(V0, {'numeric'}, {'size',[2 1],'finite'});
validateattributes(Vtrials, {'numeric'}, {'size',[2 2],'finite'});
validateattributes(Utrials, {'numeric'}, {'size',[2 1],'finite'});
validateattributes(correctionRadii, {'numeric'}, {'size',[2 1],'real','finite','positive'});
V0=double(V0); Vtrials=double(Vtrials);
Utrials=double(Utrials); correctionRadii=double(correctionRadii);
if any(abs(Utrials) == 0)
    error('balancing:ZeroTrial', 'Both trial unbalances must be nonzero.');
end
A = (Vtrials - V0) ./ Utrials.';
if rank(A) < 2
    error('balancing:RankDeficient', 'Trial columns do not identify two planes.');
end
Uc = -(A \ V0);
result = struct('influence',A, 'conditionNumber',cond(A), 'correction',Uc, ...
    'mass',abs(Uc)./correctionRadii, 'degrees',mod(rad2deg(angle(Uc)),360), ...
    'predictedResidual',V0+A*Uc);
end

代码用非共轭转置 .' 让复数试重变为行向量；换成 ' 会把试重角度取反。rank 检查只排除数值上不能独立辨识的两列；可求解不等于测量足够可靠，应同时查看返回的条件数、实测散布和配重约束。多测点与正则化扩展见 理论页。

3. 迟滞边沿：pulseTimesFromTacho.m

该函数在低阈值重新使能，在高阈值记录上升沿，并对相邻样本作线性阈值插值。起始时已经为高电平不会被误记为新一转。

function pulseTimes = pulseTimesFromTacho(signal, fs, low, high)
% Rising crossings of high, rearmed only below low. Times are seconds.
validateattributes(signal, {'numeric'}, {'vector','real','finite','nonempty'});
validateattributes(fs, {'numeric'}, {'real','scalar','finite','positive'});
validateattributes(low, {'numeric'}, {'real','scalar','finite'});
validateattributes(high, {'numeric'}, {'real','scalar','finite','>',low});
signal = double(signal(:));
fs=double(fs); low=double(low); high=double(high);
pulseTimes = zeros(numel(signal),1);
count = 0;
armed = signal(1) <= low;
for k = 2:numel(signal)
    if signal(k) <= low
        armed = true;
    elseif armed && signal(k) >= high
        fraction = (high-signal(k-1))/(signal(k)-signal(k-1));
        count = count+1;
        pulseTimes(count) = (k-2+fraction)/fs;
        armed = false;
    end
end
pulseTimes = pulseTimes(1:count);
end

阈值必须使用输入信号的实际单位。它只实现边沿检测，不自动处理漏脉冲、多个标记或不可能的转速变化；输出进入下一函数前，应用应检查时间间隔与仪器允许范围。

4. 等角度 1X 提取：tacho1x.m

x 是“样本数 × 通道数”的实矩阵，pulseTimes 是同一时基中的每转参考秒数。对每个完整转数按线性角度模型插值，直接投影到 1X，最后平均复数。输出幅值是峰值，采用余弦与相位超前约定。

function result = tacho1x(x, fs, pulseTimes, pointsPerRev)
% x: real samples-by-channels; pulseTimes: one valid reference per revolution.
% Convention: x(theta)=real(V*exp(1j*theta)); output V uses peak amplitude.
validateattributes(x, {'numeric'}, {'2d','real','finite','nonempty'});
validateattributes(fs, {'numeric'}, {'real','scalar','finite','positive'});
validateattributes(pulseTimes, {'numeric'}, {'vector','real','finite','nonempty'});
validateattributes(pointsPerRev, {'numeric'}, {'real','scalar','finite','integer','>=',8});
x=double(x); fs=double(fs); pointsPerRev=double(pointsPerRev);
pulseTimes = double(pulseTimes(:));
if size(x,1)<2 || numel(pulseTimes)<2 || any(diff(pulseTimes)<=0) ...
        || pulseTimes(1)<0 || pulseTimes(end)>(size(x,1)-1)/fs
    error('balancing:InvalidPulseTimes', 'Need complete, ordered, in-range revolutions.');
end
time = (0:size(x,1)-1)'/fs;
theta = (0:pointsPerRev-1)' * (2*pi/pointsPerRev);
kernel = exp(-1j*theta);
perRev = complex(zeros(numel(pulseTimes)-1,size(x,2)));
for k = 1:numel(pulseTimes)-1
    query = pulseTimes(k) + (theta/(2*pi)) * (pulseTimes(k+1)-pulseTimes(k));
    y = interp1(time,x,query,'linear');
    perRev(k,:) = 2/pointsPerRev * sum(y.*kernel,1);
end
meanVector = mean(perRev,1);
result = struct('perRevolution',perRev,'meanVector',meanVector, ...
    'peak',abs(meanVector),'degrees',mod(rad2deg(angle(meanVector)),360), ...
    'rpmPerRevolution',60./diff(pulseTimes));
end

pointsPerRev 应足够覆盖所保留的最高阶次，但增加插值点数不能补回原始采样缺失的带宽。该实现不包含降采样抗混叠滤波、脉冲质量判别或转内加速度估计。小于输入带宽需求的角度网格会使高阶成分混叠到低阶；实机必须先确定原始有效带宽和角度重采样策略。

5. 可复算的演示：balancingDemo.m

仿真明确使用复数叠加：Vtrials = V0 + Atrue*diag(Utrial)。它真正改变了试重后的响应，不会像反复生成同一正弦的示例那样，把纯噪声差异误当作影响系数。

function result = balancingDemo()
% Fixed complex model for teaching only. This is not a hardware acquisition.
polar = @(m,d) m.*exp(1j*deg2rad(d));
result.single = singlePlaneBalancing(polar(10,30),polar(6,45),20,100,90,100);
fprintf('Single plane: %.5f g at %.5f deg\n',result.single.mass,result.single.degrees);

% Model units: displacement um per (g*mm); true unbalance uses g*mm.
Atrue = [polar(0.008,20),polar(0.003,-35); ...
         polar(0.002,70),polar(0.009,-10)];
Utrue = [polar(15*80,45);polar(20*80,135)];
Utrial = [polar(10*80,0);polar(12*100,60)];
V0 = Atrue*Utrue;
Vtrials = V0 + Atrue*diag(Utrial);
result.dual = twoPlaneBalancing(V0,Vtrials,Utrial,[80;100]);
result.actualModelResidual = Atrue*(Utrue+result.dual.correction);
fprintf('Two planes: [%.5f, %.5f] g at [%.5f, %.5f] deg\n', ...
    result.dual.mass,result.dual.degrees);

% 50 Hz signal sampled at 25600 Hz, 11 reference pulses for 10 revolutions.
fs = 25600;
t = (0:5120)'/fs;
pulses = (0:10)'/50;
x = 3*cos(2*pi*50*t+deg2rad(30)) + 0.4*cos(2*pi*100*t)+2;
result.signal = tacho1x(x,fs,pulses,512);
fprintf('1X peak: %.5f; phase: %.5f deg\n',result.signal.peak,result.signal.degrees);
end

结果约为：

Single plane: 44.62230 g at 110.27174 deg
Two planes: [15.00000, 16.00000] g at [225.00000, 315.00000] deg
1X peak: 3.00000; phase: 30.00000 deg

第二面的原始不平衡是 20 g × 80 mm，校正半径改为 100 mm，所以需要 16 g。两个校正向量应分别抵消仿真中 45° 和 135° 的真值。合成波形还叠加了 DC 和 2X；在整圈等角网格上，它们不会进入这里的 1X 投影。

从数值例子扩展到学习实验

实验	改什么	应记录什么
复数几何	原始不平衡、试重角和安装半径	校正矢量、复数合成残差，避免只看克数
试重信噪比	采用固定随机种子，为每次响应加入已知复数噪声	重复试验误差分布，而非单条随机曲线的“最优点”
可辨识性	逐渐使矩阵两列接近共线	实际 cond(A)、奇异值与校正误差
同步测量	脉冲抖动、漏脉冲、双触发、采样/通道延迟	无效标记和 1X 幅相误差
速度变化	不同转长、转内速度变化	线性转角模型的误差与适用范围
执行误差	校正质量/角度偏差、模型变化	重新生成的实测响应与预测残差差异

条件数实验可以用 A=Q1*diag([1,1/k])*Q2' 构造已知奇异值的矩阵，其中 Q1,Q2 为酉矩阵；再核实 cond(A)。不能把 [1,.1;.1,1/k] 的参数 k 直接当作该矩阵的真实条件数。噪声、正则化与误差上界也必须在统一尺度下比较。

图形展示可用极坐标画校正矢量，用曲线画真实残余与迭代次数，用矩阵热图显示交叉影响。轴标签需要注明物理量、单位和“仿真”；人为指定的技能评分、精度曲线或固定残余序列不能充当算法性能证据。完整信号回放设计见 动平衡仿真器。

迁入 Simulink 时的边界

flowchart LR
  Source["已知真值信号 / 实测回放"] --> Frame["同步采样帧与质量标记"]
  Frame --> Ref["参考边沿与完整转数缓存"]
  Ref --> Order["等角重采样与 1X 复数"]
  Order --> State["原始 / 试重 / 校正状态"]
  State --> Solver["影响系数与受约束求解"]
  Solver --> Log["结果、条件数与原始记录"]

这是模型接口设计，不是已经接线或验证的 .slx。实现时应明确采样时间、缓冲区上限、有效长度、更新脉冲、复位和故障状态，新的 1X 结果未产生时保留上次有效值并输出有效标志，不能每个采样点都输出零冒充新测量。

原始数据源可先用可重复信号或记录回放；真实 Analog Input 块取决于所用硬件和支持包，并不存在适用于所有硬件的 simulink/Sources/Analog Input。固定步长/变步长求解器、在线执行、代码生成和实时硬件是不同配置条件，搭出几个空子系统不等于具备实时采集能力。

将算法放进 MATLAB Function 块前，需重新核验所用函数的代码生成支持、可变数组边界和执行时间；模型连线、数据类型、单元仿真与硬件接口均需独立验收。本页未执行 Simulink 或硬件验证。

验证范围

本次 23 个类式单元测试全部通过，覆盖例题数值、半径换算、保留试重、测量量一致转换、双面已知真值、复数试重转置、秩亏与零试重拒绝、起始高电平、迟滞边沿、DC/2X 隔离、跨 0° 平均、不同转长、多通道、整数输入精度，以及边沿数量/顺序/范围错误。

这些测试验证本页所给函数和合成数据，不证明任何设备的机械安全或实际平衡精度。MATLAB API 依据：线性方程求解、同步平均。