实测说明2026-10-17本文代码已在 Windows 11构建 26200、MATLAB 25.2.0.2998904R2025b、PCWIN64 中实际执行主实验及独立边界补测均通过。表格和配图来自真实 MATLAB 输出补测用 QR 求解与双边残差化交叉核对。requiredFilesAndProducts 实际分析仅返回基础 MATLAB不依赖 Curve Fitting Toolbox、Statistics and Machine Learning Toolbox 或 Signal Processing Toolbox。先期 Python 参考模型结果另行保留不混作 MATLAB 实测数据。最短运行将后文完整代码保存为detrend_joint_fit_demo.m在文件所在目录执行results detrend_joint_fit_demo;。真实关键结果先去线性趋势再拟合的幅值为 0.686491473联合拟合幅值为 0.8。测量一个已知频率的低频正弦数据中还有直流偏置与慢漂移。一个很自然的流程是先detrend再对剩余信号测幅相。本例中真实幅值为 0.8频率为 0.25 Hz。先去线性趋势再拟合正弦得到的幅值约为 0.68649低了 14.19%初相位也从 0.6 rad 变成约 0.26571 rad。这个偏差在无噪声输入上就存在。原因在于短记录里的正弦分量与常量和线性趋势通常不正交。去趋势时真实信号的一部分也被分给了“漂移”。1. 先明确已知量和待估计量适用场景是注入测试、锁定参考频率的传感器测量、已知激励频率的系统响应分析。频率f0来自可信参考待估计的是幅值、相位、偏置和漂移斜率。输入模型为x ( t ) A cos ( 2 π f 0 t ϕ ) c d ( t − t ˉ ) x(t)A\cos(2\pi f_0t\phi)cd(t-\bar t)x(t)Acos(2πf0tϕ)cd(t−tˉ)用三角恒等式展开x ( t ) a cos ( 2 π f 0 t ) b sin ( 2 π f 0 t ) c q τ x(t)a\cos(2\pi f_0t)b\sin(2\pi f_0t)cq\taux(t)acos(2πf0t)bsin(2πf0t)cqτ其中tau(t-mean(t))/timeScaletimeScalemax(abs(t-mean(t)))所以qd*timeScale。系数关系为a A cos ϕ , b − A sin ϕ aA\cos\phi,\quad b-A\sin\phiaAcosϕ,b−Asinϕ于是A a 2 b 2 , ϕ atan2 ( − b , a ) A\sqrt{a^2b^2},\qquad\phi\operatorname{atan2}(-b,a)Aa2b2,ϕatan2(−b,a)这里使用余弦相位约定。若原始定义采用正弦转换公式也要一起修改不能只把图标题改成“正弦相位”。atan2接收的参数顺序是纵坐标、横坐标返回四象限角度。MathWorks atan2 文档2. 只对数据去趋势会改变后续模型定义正弦基矩阵S[cos(2*pi*f0*t), sin(2*pi*f0*t)]趋势基矩阵Z[ones(N,1), tau]。线性去趋势可以写成yDetrendx-Z*(Z\x);betaSequentialS\yDetrend;第一行把x在常量和直线方向上的最小二乘投影去掉第二行却仍然用原来的正弦、余弦基去解释已经变形的信号。这两步一般不等于联合估计。MathWorks 将detrend描述为去除最小二乘意义下的最佳拟合直线detrend(x,0)则去均值。是否应去掉趋势要由分析目标决定。MathWorks detrend 文档本例采用Fs100 Hz、N400、f00.25 Hz标称记录时长N/Fs4 s对应一个周期的离散样本最后一个样本时刻是 3.99 s。正弦与常量项在这组采样上可以正交但它与线性时间项仍会相关。因此“已经采了整周期”不能作为任意去线性趋势后幅相不变的保证。偏差还取决于起始相位、频率和记录长度。本文的 14.19% 是指定数据的结果不是detrend的固定衰减系数更不是可直接补偿的滤波器增益。3. 把目标分量与漂移放进同一个最小二乘问题一次性建立四列矩阵X[S,Z];betaJointX\x;使用四个系数同时解释原始数据。在频率正确、模型匹配、矩阵满列秩且数值条件合理的前提下无噪声数据可在浮点精度范围内恢复参数。MATLAB 的反斜杠运算可以对长方形矩阵求最小二乘解。MathWorks mldivide 文档工程代码不必显式计算inv(X*X)*X*x。对满列秩矩阵形成法方程会把 2 范数条件数平方还增加一次显式求逆。直接使用X\x保留求解器选择稳定路径的空间。时间列先居中再缩放到约[-1,1]可以避免直接放入很大的绝对时间戳。但这只能改善列尺度不能解决“只观察到极小一段正弦时它本身就很像常量加直线”的可辨识性问题。4. 已经去趋势的数据模型基也要做同样处理如果处理流程必须先做去趋势可以对正弦基也施加相同操作SResidualS-Z*(Z\S);betaPartialSResidual\yDetrend;也就是同时对数据和模型去除同一个趋势子空间。这个结果与联合拟合的前两个系数相同适用条件是参与计算的矩阵满足所需的秩条件。从分块最小二乘可以直接推导对任意正弦系数最佳趋势系数会消去趋势方向的残差余下的优化只需要比较残差化之后的S与x。因此这里是一项有独立数学依据的交叉检查。代码没有构造N×N的投影矩阵避免为了一个四参数问题浪费二次方规模的内存。MATLAB 独立补测在 6 种周期数、4 种初相位共 24 组配置下用 QR 分解单独求解并与反斜杠、已知参数和双边残差化结果交叉核对三项最大差分别为约 1.68×10⁻¹⁴、9.44×10⁻¹⁵ 和 8.99×10⁻¹⁵。先期 Python 参考检查单独保留。5. 本次 MATLAB 实测结果检查项本次 MATLAB 结果真实幅值 / 相位0.8 / 0.6 rad仅数据去趋势后拟合0.686491473 / 0.265707073 rad顺序处理的幅值相对误差-14.1886%联合拟合幅值 / 相位约 0.8 / 0.6 rad联合拟合直流 / 斜率约 1.7 / 0.12 输入单位每秒联合模型条件数3.37292双边残差化与联合系数最大差值1.11×10⁻¹⁶显式投影与内置 detrend 最大差值4.44×10⁻¹⁶无噪声模型最大残差检查小于 1×10⁻¹²断言通过加入幅度为 0.01 的 17 Hz 正弦和幅度为 0.006 的 23 Hz 余弦后估计幅值为 0.800155757相位为 0.600289319 rad残差 RMS 为 0.008244689。它们是可重复的确定性干扰不能把这一次结果解释为白噪声统计实验或置信区间。把拟合频率故意提高 5%即使用0.2625 Hz即便原始输入没有噪声幅值也变为约 0.768046944残差 RMS 为 0.018586864。频率是这个方法的前提联合拟合不会自动修复错误的频率参数。上图为本次 MATLAB 实际导出。左下图在固定相位 0.6 rad 下改变记录内的标称周期数它展示本组数据的顺序处理偏差不应当作普遍的幅频响应。右下图展示联合模型也有自己的短窗口限制。6. 只采到 0.05 个周期联合拟合也会变得脆弱把记录内的标称周期数f*N/Fs从 1 减小到 0.05条件数从约 3.37 增加到约 14647.38。此时正弦的局部形状与低阶趋势很接近小的干扰和建模误差都可能造成较大的参数变化。cond(X)的默认 2 范数条件数由最大、最小奇异值的比值得到可用来观察这种敏感性。MathWorks cond 文档独立 MATLAB 补测还向两组模型各自最小奇异值对应的左奇异向量方向加入相对范数为1e-6的确定性输入扰动。0.05 周期时系数向量变化的范数约为 0.02367估计幅值从约 0.8 变为 0.8162351 周期时系数向量变化范数约为 6.12×10⁻⁶幅值约为 0.799998。这是专门构造的敏感性反例不是随机噪声统计结论。不存在一个对所有测量系统都合适的固定条件数门限。应结合允许的幅相误差、信噪比、参考频率稳定度和模型残差评估。条件数较小也不代表模型一定正确。若f00余弦列与直流列重复正弦列全零四个参数无法独立识别。本例还显式检查了这个降秩边界。均匀实采样在奈奎斯特频率处同样不能按普通双正交基方式任意估计幅相。实际应用应避开这些退化点。7. 完整 MATLAB 示例与运行方式保存为detrend_joint_fit_demo.m在文件目录执行resultsdetrend_joint_fit_demo;函数包含顺序处理反例、联合拟合、残差化对照、错误频率对照、短窗扫描和结果导出。本次实际日志已打印JOINT MATLAB PASS11 处断言各执行一次主程序批次进程退出码为 0。与同批 FIR 实验合计主程序有 20 次成功断言执行独立补测有 103 个通过检查另有两次 MAT 保存回读核验不混入前述计数。运行前用version和ver保存实际环境图和 MAT 文件输出到matlab_results建议每次使用新的输出目录。functionresultsdetrend_joint_fit_demo(outDir)% Base-MATLAB joint sinusoid and linear-trend least-squares demo.% Compares sequential detrending, joint fitting, and conditioning limits.ifnargin1outDirfullfile(fileparts(mfilename(fullpath)),matlab_results);endif~exist(outDir,dir),mkdir(outDir);endFs100;N400;f00.25;A0.8;phi0.6;dc1.7;slopePerSecond0.12;t(0:N-1)/Fs;tct-mean(t);timeScalemax(abs(tc));tautc/timeScale;S[cos(2*pi*f0*t),sin(2*pi*f0*t)];Z[ones(N,1),tau];X[S,Z];signalA*cos(2*pi*f0*tphi);trenddcslopePerSecond*tc;xsignaltrend;% Wrong comparison: remove the best-fit line from y but leave S unchanged.yDetrendx-Z*(Z\x);betaSequentialS\yDetrend;aSequentialhypot(betaSequential(1),betaSequential(2));pSequentialatan2(-betaSequential(2),betaSequential(1));% Joint model: cosine, sine, offset, and centered/scaled time together.betaJointX\x;aJointhypot(betaJoint(1),betaJoint(2));pJointatan2(-betaJoint(2),betaJoint(1));yJointSignalS*betaJoint(1:2);residualx-X*betaJoint;% Independent algebraic check: remove Z from BOTH y and signal columns.SResidualS-Z*(Z\S);betaPartialSResidual\yDetrend;partialErrormax(abs(betaPartial-betaJoint(1:2)));detrendErrormax(abs(detrend(x)-yDetrend));phaseErroratan2(sin(pJoint-phi),cos(pJoint-phi));assert(abs(aJoint-A)1e-12,JOINT:Amplitude,Joint amplitude mismatch.);assert(abs(phaseError)1e-12,JOINT:Phase,Joint phase mismatch.);assert(abs(betaJoint(3)-dc)1e-12,JOINT:DC,DC mismatch.);assert(abs(betaJoint(4)/timeScale-slopePerSecond)1e-12,...JOINT:Slope,Slope mismatch.);assert(partialError1e-12,JOINT:Partial,Residualized model differs.);assert(detrendError1e-12,JOINT:Detrend,Explicit detrend differs.);assert(abs(aSequential/A-1)0.1,JOINT:NegativeControl,Bias not exposed.);assert(max(abs(residual))1e-12,JOINT:Residual,Exact model mismatch.);% Reproducible extra tones, not random white noise or an uncertainty model.interference0.01*sin(2*pi*17*t)0.006*cos(2*pi*23*t);xDisturbedxinterference;betaDisturbedX\xDisturbed;aDisturbedhypot(betaDisturbed(1),betaDisturbed(2));pDisturbedatan2(-betaDisturbed(2),betaDisturbed(1));rmsDisturbedsqrt(mean((xDisturbed-X*betaDisturbed).^2));% Frequency is an INPUT to this method. A wrong frequency changes the result.fWrong1.05*f0;XWrong[cos(2*pi*fWrong*t),sin(2*pi*fWrong*t),Z];betaWrongXWrong\x;aWronghypot(betaWrong(1),betaWrong(2));rmsWrongsqrt(mean((x-XWrong*betaWrong).^2));assert(rmsWrong1e-3,JOINT:WrongFrequency,Frequency mismatch not exposed.);% Number of nominal cycles f*N/Fs. Short windows approach a line.cycles[0.050.10.250.512410];conditionNumberszeros(size(cycles));sequentialRatioszeros(size(cycles));fork1:numel(cycles)fcycles(k)/(N/Fs);Sk[cos(2*pi*f*t),sin(2*pi*f*t)];Xk[Sk,Z];xkA*cos(2*pi*f*tphi)trend;ckSk\(xk-Z*(Z\xk));conditionNumbers(k)cond(Xk);sequentialRatios(k)hypot(ck(1),ck(2))/A;endassert(conditionNumbers(1)1e4,JOINT:Condition,Short-window risk missing.);assert(rank([ones(N,1),zeros(N,1),Z])4,...JOINT:ZeroFrequency,DC must be rank deficient.);figfigure(Visible,off,Color,w,Position,[1001001050760]);subplot(2,2,1);plot(t,x,k-,t,trend,b--);grid on;xlabel(Time (s));ylabel(Input units);legend(Tone offset drift,Known drift);title(0.25 Hz tone in a 4 s nominal record);subplot(2,2,2);plot(t,signal,k-,t,yDetrend,r-,t,yJointSignal,b--);grid on;xlabel(Time (s));ylabel(Input units);legend(True tone,After linear detrend,Joint fitted tone);title(Detrending also removes part of the tone);subplot(2,2,3);semilogx(cycles,sequentialRatios,o-);hold on;plot([cycles(1),cycles(end)],[1,1],k--);grid on;xlabel(Nominal cycles f N / Fs);ylabel(Sequential amplitude / true amplitude);title(Bias depends on window and phase; phi 0.6 rad);subplot(2,2,4);loglog(cycles,conditionNumbers,o-);grid on;xlabel(Nominal cycles f N / Fs);ylabel(cond([cos sin 1 tau]));title(Joint fitting also has a short-window limit);print(fig,fullfile(outDir,detrend_joint_fit_matlab.png),-dpng,-r150);close(fig);resultsstruct(Fs,Fs,N,N,f0,f0,trueAmplitude,A,truePhase,phi,...sequentialAmplitude,aSequential,sequentialPhase,pSequential,...jointAmplitude,aJoint,jointPhase,pJoint,jointDC,betaJoint(3),...jointSlope,betaJoint(4)/timeScale,jointCondition,cond(X),...partialError,partialError,detrendError,detrendError,...disturbedAmplitude,aDisturbed,disturbedPhase,pDisturbed,...disturbedResidualRMS,rmsDisturbed,wrongFrequency,fWrong,...wrongAmplitude,aWrong,wrongResidualRMS,rmsWrong,...cycles,cycles,conditionNumbers,conditionNumbers,...sequentialRatios,sequentialRatios);save(fullfile(outDir,detrend_joint_fit_matlab.mat),results);fprintf(JOINT sequential_A%.12g; sequential_phi%.12g\n,aSequential,pSequential);fprintf(JOINT A%.12g; phi%.12g; dc%.12g; slope%.12g; cond%.12g\n,...aJoint,pJoint,betaJoint(3),betaJoint(4)/timeScale,cond(X));fprintf(JOINT partial_error%.16g; detrend_error%.16g\n,partialError,detrendError);fprintf(JOINT disturbed_A%.12g; disturbed_phi%.12g; residual_rms%.12g\n,...aDisturbed,pDisturbed,rmsDisturbed);fprintf(JOINT wrong_f%.12g; wrong_A%.12g; residual_rms%.12g\n,...fWrong,aWrong,rmsWrong);fprintf(JOINT MATLAB PASS; low_cycle_condition%.12g\n,conditionNumbers(1));end8. 接入测量脚本时的检查顺序从可信的激励源或采样参考得到f0。若频率未知或随时间变化应另做频率估计与模型选择本文不覆盖这一阶段。用真实采样时间建立正弦基和趋势基。对不等间隔时间戳显式矩阵仍可建立但默认按样本序号处理的detrend(x)不再等于按真实时间去趋势不能照搬该对照。对缺失值做一致的行筛选同时更新x、t及所有基函数。不要只删数据、不删矩阵中的对应行。检查矩阵秩、条件数、残差以及观测窗内的周期数。增大窗口常常比增加拟合阶数更有用但要确认较长窗口内频率和漂移模型仍成立。把相位参考时刻写进结果。本例相位相对于t0若平移时间原点相位也要相应换算。幅值接近噪声水平时即使atan2返回数字相位也可能没有实际可信度。保留原始数据和拟合出的漂移。真实慢变化可能就是要测的物理量不能因为“看着不平”便默认删掉。高阶去趋势会进一步侵占低频子空间必须重新验证目标信号误差。这个例子的重点是让预处理与后续模型保持一致。在有限观测窗口内先删去一部分数据方向再用未经修改的模型解释它会造成可量化的系统偏差。已知频率条件下把目标分量和干扰项放进同一个最小二乘问题能把这种偏差显式地暴露并控制起来。延伸阅读MATLAB Hann 相干增益与补零标度对照。如果后续再用 FFT 检查拟合残差仍需要独立处理窗口、幅值标度及非相干采样问题。工程问题如果你的观测窗只能覆盖 0.05 个周期而频率参考又存在小幅误差你会先延长窗口还是仍直接报告联合拟合的幅值可以沿用本文的条件数、频率误设和扰动测试先检查结果对这些变化是否稳定。