// zpkdatd.cpp// 独立实现 MATLAB 的 zpkdata 函数无第三方依赖// 编译: g -stdc17 zpkdata.cpp -o zpkdata#includeiostream#includevector#includecomplex#includecmath#includestdexcept#includealgorithm#includeiomanipusingComplexstd::complexdouble;usingComplexVectorstd::vectorComplex;usingRealVectorstd::vectordouble;// ---------- 结果结构体 ----------structZpkResult{ComplexVector zeros;// 零点ComplexVector poles;// 极点doublegain;// 增益};// ---------- Horner 法计算多项式在复数点的值 ----------// coeffs[0] 为最高次项系数降幂staticComplexpolyEval(constRealVectorcoeffs,constComplexx){Complexr(0.0,0.0);for(doublec:coeffs)rr*xc;returnr;}// ---------- 求多项式全部复根Durand-Kerner 方法 ----------staticComplexVectorpolynomialRoots(RealVector coeffs){ComplexVector roots;// 去掉前导零size_t start0;while(start1coeffs.size()coeffs[start]0.0)start;coeffs.erase(coeffs.begin(),coeffs.begin()start);intnstatic_castint(coeffs.size())-1;// 多项式次数if(n1)returnroots;// 归一化使首项系数为 1constdoubleleadcoeffs[0];for(autoc:coeffs)c/lead;// 初始猜测几何分布 (0.40.9i)^kroots.resize(n);constComplexseed(0.4,0.9);Complexcur(1.0,0.0);for(inti0;in;i){roots[i]cur;cur*seed;}// 迭代constintmaxIter2000;constdoubletol1e-13;for(intiter0;itermaxIter;iter){doublemaxDelta0.0;for(inti0;in;i){Complex numpolyEval(coeffs,roots[i]);Complexden(1.0,0.0);for(intj0;jn;j)if(j!i)den*(roots[i]-roots[j]);if(std::abs(den)1e-300)continue;Complex deltanum/den;roots[i]-delta;maxDeltastd::max(maxDelta,std::abs(delta));}if(maxDeltatol)break;}// 将极小虚部置零for(autor:roots)if(std::abs(r.imag())1e-8)rComplex(r.real(),0.0);// 在 polynomialRoots 中return roots; 之前加入// 按实部降序、虚部降序排序以匹配 MATLAB 输出顺序std::sort(roots.begin(),roots.end(),[](constComplexa,constComplexb){if(std::abs(a.real()-b.real())1e-9)returna.real()b.real();returna.imag()b.imag();});returnroots;}// ---------- zpkdata 主逻辑 ----------staticZpkResultzpkdata(constRealVectornum,constRealVectorden){if(num.empty()||den.empty())throwstd::invalid_argument(Numerator/denominator must be non-empty.);ZpkResult r;r.zerospolynomialRoots(num);r.polespolynomialRoots(den);r.gainnum.front()/den.front();returnr;}// ---------- 打印 ----------staticvoidprintComplex(constComplexc){doublerec.real();doubleimc.imag();if(std::abs(re)1e-12)re0.0;if(std::abs(im)1e-12)im0.0;std::cout(re,im);}intmain(){// G(s) 10*(s-1)*(s2) / ((s1)*(s^24s5))RealVector num{10.0,10.0,-20.0};RealVector den{1.0,5.0,9.0,5.0};ZpkResult rzpkdata(num,den);std::coutstd::setprecision(6);std::coutZeros: ;for(constautoz:r.zeros){printComplex(z);std::cout ;}std::cout\nPoles: ;for(constautop:r.poles){printComplex(p);std::cout ;}std::cout\nGain: r.gainstd::endl;return0;}