原载新浪 http://blog.sina.com.cn/s/blog_729a9214010303y6.html
试验程序ExpFitFrfit.sce如图片1.的右半部所示。用Scilab-6.1.0_x64频率响应拟合（frequency response fit）函数frfit的帮助文档中的例子（复制了它的几个程序语句），生成频率响应数据x0。把x0按信（不管大小都不排除均值）噪比要求放大后视为无噪数据，再加复值高斯白噪声，然后直接对含噪序列做极点模型拟合得z0，计算z0与无噪数据之间的相对误差（范数比），记入误差矩阵NR0，矩阵的行和列号即内和外层for循环指针分别对应噪声试验序号和给定的信噪比值的序号。假设x0有L个采样点已足够宽，其均匀采样间隔足够小，允许周期化引入误差，未知道系统传函（transfer function，传递函数；设分子和分母中系数都为实数）以及时域冲激响应。设序列y1的长度大于L+L(本短文测验取最小的2的整数次幂；添零多时，影响噪声统计特性，增加后续计算)；除第2至L+1点的一段等于含噪频响序列之外，其余值全为零。对y1做IFFT（或FFT），抛弃其变换结果的虚部，只截取其实部的2倍做极点模型拟合（第1点值不变），再对拟合结果做FFT（或IFFT）得z1，计算z1的第2至L+1点值与无噪数据之间的相对误差，记入矩阵NR1。最后，显示各信噪比（行）时误差的均值（列）、标准差、最大值。
图片1.左半部的命令窗，显示了运行结果。与最右边一列（原有噪声误差）比，数值越小越好。每一信噪比（120、70、35、25、15和5dB）用200个噪声序列。使用的极点（无再筛选）模型阶数为30、12，在信噪比大到25dB时，从这组结果看起来可行。用越少的待定参数透过现象抓关键越好。当然，人们评估误差的方式，依赖于研究和应用目的。
每次获取一个高斯分布的随机数，当它大于理想均值时，计算就用了IFFT-FFT对更符合习惯，否则用FFT-IFFT对，而其它语句不知这一选择。用了一组频域数据为例演示DFT的应用以及提取变换域内极点，但是不被模拟量和2pi纠缠；当需要量纲时，才确认变元的符号和采样间隔的物理本义以及DFT的长度点数。可以定义频带，挪动数据填充。与此类以及对应DFT中基向量的虚部全为零处的特殊有关，居士也曾提过采样和非整数移位（或插值重抽样）以及磁共振成像和连续小波变换。
若直接用两个复变元多项式之比，能拟合好频响序列，那么可实现很紧凑的数据表达。尤其允许灵活甚至随机化的尽量少的探测频率点摆脱僵化采样，那更有益。然而，观测数据常常是有噪声的，而且未知优化采样位置以及数据足以表达好了的多项式阶数。据微分方程和无阻尼复指数信号输入集，转而解线性方程组估计多项式系数，这仍未如用理想冲激函数输入时的全连续频谱均匀激励那样周全。实践中求尽可能好的近似解。
图片1.中的0.373526是frfit举例时由无噪数据辨识传函后再获得响应的误差，显得较大。其处理对于采样点的改变和噪声还太敏感。测试程序Sl6xFrfit.sce和运行结果如图片2.所示；若用[0.005:0.005:2.5]频率点，误差为0.0457656；用[0.005:0.005:2]频率点时，误差为0.518053；弱尾部后或有大影响。[0.01:0.01:2]是文档中原有的频率点，在120dB信噪比时，误差约为0.372与无噪时接近，但是，在70dB时，程序就出现警告滤波器阶数太低，或算出荒谬结果，甚至崩溃了。从文档中的曲线图已易见幅谱尖峰或谷点附近的差别。用analpf生成些滤波器的频率响应（文档中绿线经不同半功率点）测试frfit，居士也感到它处理有纹波的频率响应较困难（未必改善了设计削弱纹波）。
可选用控制系统（Control Systems-CACSD）模块的函数frep2tf而绕过信号处理（Signal Processing）模块的frfit，从含噪频响估计连续时间系统传函，再用这个传函求频响z，考察z与无噪数据之间的相对误差。为此，紧接着图片2.的程序结束之后，在命令窗中运行：
rand('seed',8); [T0,T1]=TwoFr2tf(25,200,0,1,4);
rand('seed',8); [T0,T1]=TwoFr2tf(120,200,0,1,4);
可见结果如图片3.所示。在120dB时，快速得到了很好（包括传函分母多项式系数的实质序列之间的误差）结果。在25dB时frep2tf用其默认加权迭代设置和repfreq配对，看起来也还好。它作为应用方便的工具的话，程序需些完善。
第12行突显了与第一段中测验不同。保持无噪数据不变，按信噪比要求改变噪声方差，不然可能遇到辨识程序的其它问题。反复运行frep2tf帮助文档中辨识随机系统的例子，也可见很漂亮曲线重合，但是，如图片3.下部所示，用rand('seed',1237023648)设置随机数发生器，并把例中语句[Sys2,err]=frep2tf(frq,rep,10)改为[Sys2,err]=frep2tf(frq,(30*rep),10)即提高了增益30倍而已，“紧接着”运行时，frep2tf因错误而终止了，未能近似处理，且提示：误用inv求逆矩阵，其输入矩阵类型错了。以便顺畅地试验那个传函、几个滤波器和随机系统以及超阶逼近、增益范围和频点分布等，居士试写了点代码以较直接些的方式求解方程组，有趣的是，在信噪较高时，也常见了很小的拟合误差。
FFTW使基础且应用广泛的DFT或FFT的计算复杂化了（《无阻尼单指数序列的FFT频谱泄漏与常值序列的超阶拟合误差》，2017-08-25）。用资源管理器见到两版（默认完全安装6.0.0和6.1.0，无互联网）的bin\fftw目录中的dll文件大小都为2662451字节（2601KB）。若不用默认的fftw_flags值64，易见6.1.0中fft错误地输出了零序列；或在启动之前，干脆先更改bin\fftw的名字以避开其dll，有时就见变换错误地直接输出了输入。处理实数序列，易测到与5.3.3版的差别。这些可能使一些人们对有“fft”这种老东西的程序感到紧张，即使在用默认设置时对有限例子已逐一监测确认了每次调用结果是可取的。
图片1.中装载的函数子程序集SsSvdLfitExp.bin，仍是居士2018年在Scilab-6.0.0制作的，但是，2017年从县城买的那台电脑去年就坏掉了。现在用lenovo笔记本电脑，Windows 7 旗舰版，Service Pack 1，64位操作系统，AMD E1-1200 APU,1400Mhz，4GB RAM。运行format('e',25);%eps,cos(1+%i)，见到了
%eps =
2.220446049250313081D-16
ans =
8.337300251311489108D-01 - 9.888977057628650646D-01i，与《同机运行两种(64位6.0.0及32位5.3.3)Scilab小测兼容误差》（2017-07-08）比较，可知系统又变了。涉及多项式，也易测到变异。5.3.3版称均匀分布随机数区间为(0,1)，6.1.0改用[0,1)。不过，用同机6.0.0平台运行第一个程序，从命令窗里看见相同误差值。
新浪赛特居士SciteJushi-2020-09-18。
图片1.DFT对偶域内极点模型拟合的试验程序和结果
图片2.测试frfit对频率点和噪声的敏感
图片3.增益和随机数或使frep2tf崩溃
转载本文请联系原作者获取授权，同时请注明本文来自唐白玉科学网博客。
链接地址：https://wap.sciencenet.cn/blog-530150-1251104.html?mobile=1
收藏
