2025年Safari与iOS版本检测新技术解析
2026/10/31 11:37:33
// ============================================================================// damp.cpp —— MATLAB damp 函数的 C++ 复刻(单文件整合版)//// 由 damp.hpp + damp.cpp 合并而成,自包含、无第三方依赖,可直接编译运行。//// 对应 MATLAB: [wn, zeta] = damp(sys)// - wn : 极点自然频率 (rad/s),wn = |p|// - zeta: 阻尼比,zeta = -cos(angle(p)) = -real(p)/|p|// - tau : 时间常数,tau = 1/(zeta*wn) = -1/real(p)(连续域)// - 离散系统按 s = log(z)/Ts 映射到连续域后再计算(与 MATLAB 一致)// - 传递函数输入时用 Durand-Kerner 复根迭代求极点(自包含)//// 用法:// // 1) 直接编译运行(默认含 main 演示与验证)// g++ -std=c++17 -O2 damp_single.cpp -o damp_demo && ./damp_demo//// // 2) 作为库使用(-DDAMP_LIB 去掉 main,接口见下)// g++ -std=c++17 -O2 -DDAMP_LIB -c damp_single.cpp -o damp_single.o// auto res = damp_ns::dampFromPoles(poles, Ts); // Ts>0 离散,Ts<=0 连续// auto res = damp_ns::dampFromTransferFunction(num, den, Ts); // 由传递函数求极点// damp_ns::printDampTable(res);//// 验证用例(main 中自动断言):// 极点 -1±2i -> wn=sqrt(5)=2.2361, zeta=0.4472, tau=1.0// 实极点 -5 -> wn=5, zeta=1.0, tau=0.2// 离散 z=0.9, Ts=0.1 -> s=log(0.9)/0.1=-1.0536, wn=1.0536, zeta=1.0// 积分极点 0 -> zeta=-1, wn=0(MATLAB 惯例)// ============================================================================#include<complex>#include<vector>#include<cmath>#include<cstdio>#include<iomanip>#include<iostream>#include<sstream>namespacedamp_ns{// ---------------------------------------------------------------------------// 接口声明(原 damp.hpp 内容)// ---------------------------------------------------------------------------// 单个极点的阻尼分析结果structDampResult{std::complex<double>pole;// 极点(离散时已用 s=log(z)/Ts 换算)doublewn;// 自然频率 (rad/s)doublezeta;// 阻尼比(实轴极点 zeta=1;原点极点按 MATLAB 惯例 zeta=-1)doubletau;// 时间常数 (s),zeta==0 或 wn==0 时为 inf};// 由极点直接计算 wn/zeta/tau// Ts > 0 : 视为离散极点 z,先 s = log(z)/Ts 再计算// Ts <= 0: 视为连续极点 s,直接计算std::vector<DampResult>dampFromPoles(conststd::vector<std::complex<double>>&poles,doubleTs=0.0);// 由传递函数 num(s)/den(s) 计算极点(den 最高次系数在前),再计算阻尼指标std::vector<DampResult>dampFromTransferFunction(conststd::vector<double>&num,conststd::vector<double>&den,doubleTs=0.0);// 按 MATLAB 格式打印结果表格voidprintDampTable(conststd::vector<DampResult>&results);// ---------------------------------------------------------------------------// 实现(原 damp.cpp 内容)// ---------------------------------------------------------------------------namespace{// 把一个极点换算到 s 域并计算阻尼指标(核心,与 MATLAB 相同公式)DampResultanalyzePole(std::complex<double>p,doubleTs){DampResult r;r.pole=p;if(Ts>0.0){// 离散:s = log(z)/Tsif(std::abs(p)==0.0){// 原点极点:log(0) -> -inf,视为无穷大频率、临界阻尼r.pole=std::complex<double>(-std::numeric_limits<double>::infinity(),0.0);r.wn=std::numeric_limits<double>::infinity();r.zeta=1.0;r.tau=0.0;returnr;}p=std::log(p)/Ts;r.pole=p;}constdoublemag=std::abs(p);if(mag==0.0){// 连续原点极点:MATLAB 惯例显示 Damping=-1, Frequency=0r.wn=0.0;r.zeta=-1.0;r.tau=std::numeric_limits<double>::infinity();returnr;}r.wn=mag;r.zeta=-p.real()/mag;// = -cos(angle(p))r.tau=(std::abs(r.zeta)<1e-12)?std::numeric_limits<double>::infinity():1.0/(r.zeta*r.wn);returnr;}// Durand-Kerner 复根迭代:求首一多项式 roots,系数 a[0] 为最高次项std::vector<std::complex<double>>polyRoots(conststd::vector<double>&a){constintn=static_cast<int>(a.size())-1;std::vector<std::complex<double>>roots;if(n<=0)returnroots;// 去首零std::vector<double>coeff=a;while(!coeff.empty()&&std::abs(coeff.front())<1e-15){coeff.erase(coeff.begin());}if(coeff.empty())returnroots;constintdeg=static_cast<int>(coeff.size())-1;if(deg==1){// 一次:-c1/c0roots.push_back(-coeff[1]/coeff[0]);returnroots;}// 初始猜测:均匀分布在半径 R 的圆上constdoubleR=1.0+std::max(1.0,std::abs(coeff[1]/coeff[0]));roots.resize(deg);for(intk=0;k<deg;++k){doubleang=2.0*M_PI*(k+0.5)/deg;roots[k]=std::complex<double>(R*std::cos(ang),R*std::sin(ang));}// 迭代:r_i = r_i - P(r_i) / prod_{j!=i}(r_i - r_j)for(intiter=0;iter<500;++iter){doublemaxMove=0.0;std::vector<std::complex<double>>next(roots);for(inti=0;i<deg;++i){std::complex<double>num=0.0,den=1.0;// 多项式求值:P(z) = coeff[0]*z^deg + ...std::complex<double>pz=coeff[0];for(intj=1;j<=deg;++j)pz=pz*roots[i]+coeff[j];num=pz;for(intj=0;j<deg;++j){if(j!=i)den*=(roots[i]-roots[j]);}if(std::abs(den)<1e-300)continue;// 避免除零next[i]=roots[i]-num/den;maxMove=std::max(maxMove,std::abs(next[i]-roots[i]));}roots=next;if(maxMove<1e-12)break;}returnroots;}}// namespacestd::vector<DampResult>dampFromPoles(conststd::vector<std::complex<double>>&poles,doubleTs){std::vector<DampResult>out;out.reserve(poles.size());for(constauto&p:poles)out.push_back(analyzePole(p,Ts));returnout;}std::vector<DampResult>dampFromTransferFunction(conststd::vector<double>&num,conststd::vector<double>&den,doubleTs){(void)num;// damp 只关心极点(分母根),分子不影响 wn/zetaautoroots=polyRoots(den);returndampFromPoles(roots,Ts);}voidprintDampTable(conststd::vector<DampResult>&results){std::cout<<"\n"<<" Pole Damping Frequency Time Constant\n"<<" (rad/seconds) (seconds)\n";std::cout<<std::scientific<<std::setprecision(2);for(constauto&r:results){// 极点格式化:-1.00e+00 + 2.00e+00istd::ostringstream oss;constdoublere=r.pole.real(),im=r.pole.imag();if(std::isinf(re)){oss<<"-Inf";}elseif(std::abs(re)<1e-15){oss<<"0.00e+00";}else{oss<<std::scientific<<std::setprecision(2)<<re;}if(!std::isinf(re)){oss<<(im>=0?" + ":" - ")<<std::scientific<<std::setprecision(2)<<std::abs(im)<<"i";}constchar*wn_s=std::isinf(r.wn)?" Inf":nullptr;constchar*tau_s=std::isinf(r.tau)?" Inf":nullptr;std::cout<<std::left<<std::setw(18)<<oss.str()<<" "<<std::right<<std::setw(10)<<r.zeta<<" ";if(wn_s)std::cout<<wn_s;elsestd::cout<<std::setw(12)<<r.wn;std::cout<<" ";if(tau_s)std::cout<<tau_s;elsestd::cout<<std::setw(12)<<r.tau;std::cout<<"\n";}std::cout<<std::defaultfloat<<std::setprecision(6);}}// namespace damp_ns// ============================================================================// 演示与验证(main,默认编译;-DDAMP_LIB 时排除 main 作库使用)// ============================================================================#ifndefDAMP_LIBintmain(){usingnamespacedamp_ns;constdoubleeps=1e-4;std::cout<<"===== 例1:连续系统极点 -1 ± 2i(G(s)=5/(s^2+2s+5))=====";std::vector<std::complex<double>>p1={std::complex<double>(-1,2),std::complex<double>(-1,-2)};autor1=dampFromPoles(p1);printDampTable(r1);// 期望:wn=sqrt(5)≈2.2361, zeta=0.4472, tau=1.0boolok1=std::abs(r1[0].wn-std::sqrt(5.0))<eps&&std::abs(r1[0].zeta-1.0/std::sqrt(5.0))<eps&&std::abs(r1[0].tau-1.0)<eps;std::cout<<"\n===== 例2:实极点 -5(一阶惯性)=====";autor2=dampFromPoles({std::complex<double>(-5,0)});printDampTable(r2);// 期望:wn=5, zeta=1, tau=0.2boolok2=std::abs(r2[0].wn-5.0)<eps&&std::abs(r2[0].zeta-1.0)<eps&&std::abs(r2[0].tau-0.2)<eps;std::cout<<"\n===== 例3:传递函数 den = s^2 + 2s + 5(应同例1)=====";autor3=dampFromTransferFunction({5.0},{1.0,2.0,5.0});printDampTable(r3);boolok3=r3.size()==2&&std::abs(r3[0].wn-std::sqrt(5.0))<eps;std::cout<<"\n===== 例4:离散极点 z=0.9, Ts=0.1(s=log(0.9)/0.1≈-1.0536)=====";autor4=dampFromPoles({std::complex<double>(0.9,0)},0.1);printDampTable(r4);doubles_exp=std::log(0.9)/0.1;// ≈ -1.0536boolok4=std::abs(r4[0].wn-std::abs(s_exp))<eps&&std::abs(r4[0].zeta-1.0)<eps;std::cout<<"\n===== 例5:积分极点 0(MATLAB 惯例 zeta=-1)=====";autor5=dampFromPoles({std::complex<double>(0,0)});printDampTable(r5);boolok5=std::abs(r5[0].wn-0.0)<eps&&std::abs(r5[0].zeta+1.0)<eps;std::cout<<"\n===== 验证结论 =====\n";std::cout<<(ok1?"[PASS] 例1 -1±2i : wn=sqrt5, zeta=0.4472, tau=1.0\n":"[FAIL] 例1\n");std::cout<<(ok2?"[PASS] 例2 -5 : wn=5, zeta=1, tau=0.2\n":"[FAIL] 例2\n");std::cout<<(ok3?"[PASS] 例3 传递函数极点与例1一致\n":"[FAIL] 例3\n");std::cout<<(ok4?"[PASS] 例4 离散 s=log(z)/Ts 换算正确\n":"[FAIL] 例4\n");std::cout<<(ok5?"[PASS] 例5 原点极点 zeta=-1, wn=0\n":"[FAIL] 例5\n");return(ok1&&ok2&&ok3&&ok4&&ok5)?0:1;}#endif// DAMP_LIB