科学网—标签

相关帖子	版块	作者	回复/查看	最后发表

热度 1 youmingqing 2014-2-3 10:34

利用初等数学确定铅球抛高、抛程和最佳抛角的关系，介绍了曲线簇的包络线；考虑手臂长度、滑行助跑因素,证明抛体以最佳抛角到达指定距离或高度时沿初始抛掷方向的速度为零。抛角对抛投效果影响较小，投掷铅球时应选择合适的姿态以求得较大的出手速度。图片来自网络,不做商业运用,谨此致谢铅球最佳出手角的初等解法及讨论尤明庆抛点和落点高度相同时中学物理可确定抛角 45 0 时抛射距离最远；而铅球掷远的抛点高于落点，确定最佳抛角通常利用导数知识。龚劲涛等给出一个初等解法，但未给出数学证明。抛体问题较为常见，也有许多实际应用，如飞石系绳，最好能为中学生所了解。 1 抛程、抛高与最佳抛角铅球以速度 V 、角度 θ 从原点 O 抛出，其运动轨迹是 x = Vt cos θ (1) y = Vt sin θ － gt 2 /2 (2) 消去时间参量 t 得到抛物线方程 y = x tan θ － x 2 (1+tan 2 θ ) / 4 H 0 (3) 式中 H 0 = V 2 /2 g 即垂直上抛高度。上式可改写为 x 2 / 4 H 0 = H 0 － y － H 0 (1 － x tan θ / 2 H 0 ) 2 ≤ H 0 － y (4) 给定抛高 y = H （铅球出手高度的负值），抛程 x 达到最大值 S = 2sqrt ( H 0 ( H 0 － H )) (5) 时，最佳抛角为 tan θ = 2 H 0 / S = 1 / sqrt (1 － H / H 0 ) (6) 投掷铅球时抛高 H 小于 0 ，因而最佳抛角小于 45 0 。 2 抛物线簇的包络线在给定铅球出手速度后，不同抛角得到不同的抛物线，即公式 (3) 构成一簇曲线。图 1 给出了抛角 0~90 0 之间每隔 15 0 的抛物线，几何尺度以垂直上抛高度 H 0 无量纲化。最佳状态即公式 (4) 右侧等号成立时构成曲线簇的包络线：给定 y 时曲线簇中所有 x 的最大值，或给定 x 时所有 y 的最大值。该包络线恰巧是在 y = H 0 处以速度 V 平抛的轨迹，在整个抛物线簇 (3) 右上方，且与每一根抛物线都有一点相切（图）。据此可解决坡面上的抛射问题：求出坡面方程与包络线的交点，即最远抛投位置，再由公式 (6) 确定最佳抛角。若铅球的出手速度为 V = 10 m/s, 则在 g = 10 m/s 2 时， H 0 = 5 m ；出手高度 2 m ，即 H = - 2 m ，则最佳抛角为 40.2 0 , 最大抛程为 11.83 m 。若以 45 0 抛掷，相应抛程是 11.71 m, 比最大值小 0.12 m ，相差仅 1% 。如图所示，最佳抛角小于 45 0 之后，抛物线与包络线在很大范围内差别不大，即抛角偏离最佳值对抛高或抛程的影响并不显著。 3 一般情形下铅球的最佳出手角度对于抛掷铅球而言，运动员以速度 U 滑步助跑；且铅球出手高度随抛角 θ 而变化，与臂长 b 相关。以肩关节为坐标原点，铅球的运动轨迹是： x = b cos θ + ( U + V cos θ ) t (7) y = b sin θ + Vt sin θ － gt 2 /2 (8) 包络线与参数方程 (7) 、 (8) 表示的曲线相切，因而 x 对 t 、 y 对 θ 的偏导数之积等于 x 对 θ 、 y 对 t 的偏导数之积，有 ( V + U cos θ ) － gt sin θ = 0 (9) 这表明抛体以最佳抛角达到指定距离或高度，即抛物线与包络线的切点处，沿初始抛掷方向的速度正好减少到零。文献给出的论断是 U = 0 时的特例。又，由公式 (7) 解得 t , 代入公式 (8) ，得到铅球运动的抛物线方程，再依据给定 x 下确定 θ 使 y 达到极值，即 y 对 θ 的导数为零亦可得到公式 (9) 。由公式 (9) 解得 t ，代入公式 (7) 、 (8) ，可得抛物线簇的包络线方程（以最佳抛角 θ 为参变量），进而可以确定给定抛距、抛高或抛距和抛高的关系（如坡面上抛投）的最佳抛角。定性地说，铅球的落点低于抛点，增加了空中飞行时间，抛角稍小于 45 0 可增大水平速度以利用较多的飞行时间；而滑步速度的存在又要求适当增大抛角，从而增加铅球在空中飞行时间以利用较大的水平速度。综合结果是最佳抛角仍在 45 0 左右。具体计算不再给出。 4 结语确定铅球最佳抛角的初等方法及包络线的概念，或许可由高中物理老师向同学介绍。不过，抛角对抛投效果的影响较小，投掷铅球时应选择合适的姿态以求得较大的出手速度。 1 朱照宣，周起钊，殷金生 . 理论力学 ( 上 ) . 北京：北京大学出版社， 1982. 116-117. 2 张培和 . 斜抛运动的最佳角度的选择 . 南平师专学报， 1996 , (4): 74-78. 3 楚安夫 . 关于斜抛运动的分析 . 大学物理， 1997 , 16(9): 46-47. 4 徐月明 . 铅球运动最佳出手角的理论分析和实验模拟研究 . 大学物理 , 2004, 23(2): 21-24. 5 龚劲涛，任全红，冯文林 . 巧解铅球最佳出手角 . 力学与实践， 2011 , 33(5): 61-62. 6 尤明庆 . 飞石系绳的力学分析 . 力学与实践， 2011 , 33(5): 93-95. 尤明庆 . 对铅球最佳抛角的讨论 . 力学与实践， 2013,35(1): 67-68

个人分类: 力学科普|13761 次阅读|5 个评论

39、(八)圆满收官_基于IMF自身的瞬时频率估计

baishp 2012-12-2 15:16

在上篇博文的图38-14中，很容易看到，对应于原信号的每一个切点处，瞬时频率都有一个小小的突变。放大x轴，可以看到突变范围都在切点及切点左右各一点共3点内。从图38-11、图38-15中可以看得很清楚。这是由于反余弦函数在1、-1点处的高度非线性造成的，当然属于信号处理过程中所产生的“虚假信息”，而不是信号本身所含有的信息。本人既然号称要追求连续光滑的瞬时频率估计，那就要一不做二不休，精益更求精。线性插值后，图38-14、图38-15的变化见下面图39-1、图39-2。图39-1 图39-2 在运行各列信号进行瞬时频率估计时，忽然发现有一列信号的瞬时频率右端部有很大的突变，见图39-3、图39-4。这显然是一种端点效应。图39-3 图39-4 我已经做了信号延拓工作，为什么还有端点效应呢？我的信号延拓是在综合振幅（综合切包络线）函数comprehensive_amplitude( x )中进行的，输出时截断至原信号长度。在瞬时频率估计中就没再做信号延拓。估计就是这个原因造成的。于是取消函数comprehensive_amplitude( x )中的信号延拓，改在瞬时频率估计程序中进行。最终的结果见下图，端点效应没有了。图39-5 但是在另一列信号中，信号延拓位置的改变，对瞬时频率的端点效应的影响却不是太大。见图39-6、图39-7、图39-8。图39-6 图39-7 图39-8 这个我分析是信号延拓的“质量”不好，在原信号与延拓信号的连接处不够光滑造成的。在博文《34、(四)信号延拓_基于IMF自身的瞬时频率估计》的 “2、调用上述函数，并将待分析的IMF原函数左、右两个端点分别“混入”上面左、右极值点中”中，只采用了左、右两个端点作为样条插值函数中的已知数据点。现在将原信号左右端部各3点共6点作为样条插值函数中的已知数据点，进行信号延拓。使用此信号延拓函数进行的瞬时频率估计，结果见下图。可见端点效应已经最大程度地被抑制了。下药对症！图39-9 至此，程序运行中所能碰到的所有问题都已解决。当然，无法百分百地保证将来不会碰到任何新问题，但至少目前已经顺利运行了几十列信号，都没有问题。将程序以文件名instantaneous_physical（瞬时物理量）保存为函数文件。调用格式 = instantaneous_physical( x ) 中，输入x是单分量IMF，输出tamp是瞬时振幅，tf是瞬时频率，phase是瞬时相位，y是复信号虚部。需要说明的一点是，这里由x、y组成的复信号x+j*y是所谓基带信号，不再是解析信号了。只要有了一列IMF，等于同时拥有了与之对应的瞬时振幅、瞬时频率、瞬时相位与基带信号。这个方法，从技术层面上来讲，信号信息的损失几近于0。如果结果不够准确，其原因也大多产生在信号处理的其它阶段。将此函数代替《32、(二)时频连续_基于IMF自身的瞬时频率估计》程序中的函数 = afyimf( x )，所得各IMF分量瞬时频率如下4图：图39-10 imfMBrd各分量瞬时频率图39-11 imfMBrd各分量瞬时频率图39-12 imf_TW486012各分量瞬时频率图39-13 imf_TW486012各分量瞬时频率可见瞬时频率曲线较《32、(二)时频连续_基于IMF自身的瞬时频率估计》中曲线光滑得多。基于IMF自身的瞬时频率估计现在完成了，心中有点感叹。这个方法，道理很简单，为什么没有人做？我觉得主要是人们思维习惯造成的。几十年来，Hilbert先生对人们的影响实在是太大了，一提到瞬时频率，几乎无人不立刻想到Hilbert变换、解析信号等。连思想极其活跃开放的黄院士，一发明EMD，便立刻将其与之挂上了钩，提出了HHT。虽然信号处理基本理论已经说明，解析信号的主要作用是用来作数学分析，实际得到它是很困难的。matlab中的hilbert()函数，我觉得主要是用来求高频窄带信号的解析信号，对低频信号，所求得的解析信号是很不准确的。而EMD，要将原信号分解到最低频的频率分量。因此，用解析信号法求各IMF的瞬时频率，先天理论上就有缺陷。有了本人所编的instantaneous_physical程序，HHT这个名词可以送进历史博物馆了。要不就改成HBT吧，让俺“蚂蟥叮上鹭鸶腿”，也尝尝名扬天下的滋味，哈！（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de01017b3r.html 首发时间：2012-07-11 23:22:57）

个人分类: 斤斤计较|3580 次阅读|0 个评论

38、(七)相应问题_基于IMF自身的瞬时频率估计

baishp 2012-12-2 13:57

将综合切包络线代替《31、(一)根本思路_基于IMF自身的瞬时频率估计》的（分段）极点包络线，将遮掩码c12、c34原来由IMF差分函数求取，改为由络基信号（IMF减去包络线所得信号）差分函数求取。所有的信号都可以顺利运行。见图37-1、图37-2。后者第3个子图瞬时频率，已经没有因为IMF与极点包络线相割造成的突变了。分析信号为imfMBrd第11分量。图38-1 图38-2 用切包络线代替极点包络线，产生的一个最大问题是：在两个相邻过零点之间，可能存在两个（或以上。理论如此，实际未见过）切点，而不像极点包络线那样只存在一个极值点。两个切点之间必定有一个络基极值点（络基信号除切点外的极值点）。以这3个点为界限，形成信号源在1、4或2、3两个象限交替出现的现象。见图37-3。我称之为“象限交访”现象吧。图38-3 即使只有一个切点，但仍然可能存在络基极值点。见图37-4、图37-5。它对信号处理的影响，与两个切点是一样的。图38-4 图38-5 上面3图中，络基极值点处瞬时频率的突变，见下图。图38-6 处理方法：通过对遮掩码的特别赋值，将上述情况下，两个相邻过零点之间划分为仅仅两个象限。切点数为两个时，以络基极值点为界限；切点数为一个时，以切点为界限。见图37-7、图37-8、图37-9。图38-7 图38-8 图38-9 经过上述处理后，包括上面3个络基极值点在内的各络基极值点处的瞬时频率突变基本上消失了。见图37-10。上面图37-7中的络基极值点处的瞬时频率见图37-11。图38-10 图38-11 但是在某些特别情况下，络基极值点处的瞬时频率突变还是很大的。见图37-12、图37-13。不过经过许多信号的运行，发现只有一个切点时，络基极值点处瞬时频率突变都很小，基本可以忽略。图38-12 图38-13 处理：强令信号与切包络线在两个切点之间的值相等，即令余弦值等于1（1、4象限）或-1（2、3象限）。此举是将信号源看成是在沿实轴运动，相位恒等于0（2*pi）或pi，瞬时频率恒等于0。见图37-14。图38-14 图37-7中络基极值点处、两个切点之间的瞬时频率如下图：图38-15 至此，因使用综合切包络线作振幅所产生的最大的“问题”处理完毕。不过我对这样的处理方式是存有疑虑的。从道理上来讲，信号处理，需要尽量保存信号本身所含有的信息，而尽量减少因信号处理方式不当所造成的虚假信息。但实际上这二者往往难以区分。在物理实际中，信号源的“象限交访”现象肯定是存在的。在本文中，即使它是“因信号处理方式不当所造成的虚假信息”，那也是上一阶段的信号处理方式（EMD）所遗留下来的。是不是应该在上一阶段的信号处理中采取某种统一的措施来消除它，而不应该在本阶段中采取“头痛医头脚痛医脚”的方式进行处理呢？还有几个次要问题下篇再谈。（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de010178yi.html 首发时间：2012-07-09 20:33:18）

个人分类: 斤斤计较|3929 次阅读|0 个评论

36、(六)四种包络_基于IMF自身的瞬时频率估计

baishp 2012-11-30 16:37

上篇博文谈“切包络线”，文末用此“切包络线”求了信号“切包络线综合振幅”。但那个方法是未经过深思熟虑的方法、经过几次实际运行，就知道是不行的。原因在于：当信号为低频IMF时，信号取绝对值后，原极大值仍是极大值，原极小值变成极大值；但当信号为高频IMF、特别是那种在两个极值点之间振荡的IMF时，情况就不一样了。信号取绝对值后，极大值不一定还是极大值，极小值也不一定变成极大值。见图36 -1。图36 -1 既然极值点的极值属性都不能得到保障，那么以极点包络线为基础求切包络线，当然就行不通了。那个方法既然已经写成博文，就不删除了，作为研究历程保留放那吧。下面是求“综合切包络线”(即“切包络线综合振幅”)可行的方法步骤： 1、将信号两端延拓； 2、求延拓后信号极大、极小值点； 3、将上述极小值点取绝对值，并与极大值“混编、排序”，得信号“综合极值点”。 4、对上述“综合极值点”作样条插值，得“综合极点包络线”，它也是信号“上综合极点包络线”。本博中二名词含义相同。 5、将上述“综合极点包络线”取负，得“下综合极点包络线”。 6、将信号分别减去“上、下综合极点包络线”，分别得上、下“络基信号”； 7、分别求上下“络基信号”的极大值与极小值； 8、分别去除上述极大值与极小值点中，小于、大于“上、下综合极点包络线”的极大值与极小值，剩下来的极大值与极小值，就是阶段性上下切点； 9、将上述阶段性下切点代替上面步骤“3”中的极小值点，上切点代替极大值点，重复至步骤“8”； 10、将“步骤8”所得阶段性上下切点，与上次所得阶段性上下切点进行比较。如果不完全相同，则重复步骤“9－3－4－5－6－7－8”；如果完全相同，则停止循环； 11、将步骤“8”最后所得上、下切点对应的“综合切点包络线”，截去两端延拓部分，就是最后得到的信号的“综合切包络线”。 12、去除最后步骤“8”所得上、下切点中，位于信号两端延拓部分中的上下切点，剩下来的就是原分析信号上下切包络线的上、下切点。根据上面算法编制程序 = comprehensive_amplitude( x ) x是待分析输入IMF，输出tenvb,tenva,tpoib,tpoia分别为“下综合切包络线”、“上综合切包络线”、“下切点”、“上切点”。作图看看： x=imfMBrd(1,:)， = comprehensive_amplitude( x ); （为作图方便，上面输出全部临时改为截去信号延拓部分之前对应的变量。） %下面是为了绘图而重复运行上面函数内指令 = signal_extension(x); lxr= ; ltr= ; llxr = length(lxr); ttt=1:llxr; zs=zeros(1,llxr); = extr(lxr);% %上面是为了绘图而重复运行 fg1=figure(1); set(fg1,'position', ) set(gca,'position', ); plot(t,x) hold on% plot(tl,xl,'r-')% plot(tr,xr,'r-')% plot(ltr,tenva,'k--')% plot(ltr,tenvb,'k--')% plot(indmin+tl(1)-1,lxr(indmin),'ko')%% plot(indmax+tl(1)-1,lxr(indmax),'ko')%% plot(tpoib+tl(1)-1,lxr(tpoib),'.k')% plot(tpoia+tl(1)-1,lxr(tpoia),'.k')% plot(ltr,zs,'k-.') axis( )% title('imfMBrd之1分量延拓后综合切包络线及其切点、极值点') 运行，得：图36 -2 x=imfMBrd(i,:)，分别再令i=4、7、10、13，并改变坐标轴显示范围，重复运行上面的程序，作图如下。这次将功折罪，多运行几列信号，多截几张图吧。图36 -3 图36 -4 图36 -5 图36 -6 可见包络线都没有问题，既与信号相切，又严格对称。至此，我们已经知道了，对信号可以作出4种包络线。分别简述如下： 1、极点包络线，或称“极值点包络线”，简称“包络线”，就是传统的、EMD当中用来进行迭代运算求取IMF的包络线。适用于任意信号。 2、综合极点包络线，或称“综合极值点包络线”，就是在我的博文“33、(三)综合极点_基于IMF自身的瞬时频率估计”文末求取过的包络线。它将大小极值点统一取绝对值后，再作样条插值求包络线。当时把它叫“综合振幅”或“综合包络线”，现在看来这样叫不够准确。它只适用于IMF。 3、切包络线，或称“切点包络线”。就是上篇博文中所求包络线。适用于任意信号。 4、综合切包络线，或称“综合切点包络线”。由于它就是我所要求的IMF瞬时振幅，，所以我以后也可以把它叫做“切包络线综合振幅”，或简称“综合振幅”、“IMF瞬时振幅”。它只适用于IMF。关于信号包络线问题就研究到这里。后面该用它来求信号瞬时频率了。（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de0101700a.html 首发时间：2012-06-30 13:47:40）

个人分类: 斤斤计较|4754 次阅读|0 个评论

35、(五)切包络线_基于IMF自身的瞬时频率估计

热度 1 baishp 2012-11-30 15:32

各位博友：端午节快乐！粽子飘香，亲情满堂！言归正传。求IMF“切包络线”方法步骤如下： 1、将信号两端延拓； 2、求延拓信号极大、极小值点； 3、根据上述极大、极小值点求信号上下极点包络线； 4、将信号分别减去上下极点包络线，分别得上下“络基信号（我创立的名词，表示以包络线为0基线的信号）”； 5、分别求上下“络基信号”的极大值与极小值； 6、分别去除上述极大值与极小值点中，小于、大于上下包络线的极大值与极小值，剩下来的极大值与极小值，就是阶段性上下切点； 7、分别对上述阶段性上下切点作样条插值，所得插值函数为阶段性上下切包络线； 8、将上述阶段性上下切包络线代替上面步骤“4”中的上下极点包络线，重复至步骤“6”； 9、将步骤“6”所得阶段性上下切点，与上次所得阶段性上下切点进行比较。如果不完全相同，则重复步骤“7－8－4－5－6”；如果完全相同，则停止循环； 10、将步骤“7”最后所得插值函数，截去两端延拓部分，就是最后得到的信号的上下切包络线。 11、去除步骤“6”最后所得上下切点中，位于信号两端延拓部分中的上下切点，剩下来的就是原分析信号上下切包络线的上下切点。上面只是大概的步骤，有些细节问题就不列出了。如果分析信号不是IMF，而是任意信号，那就还有更多的细节问题需要处理。由于IMF切包络线在接近“0”时，其对称性不总是太好（见图35-1、图35-2），当时我还没有想到信号取绝对值的方法，就总是想到用切包络线代替EMD中的极点包络线进行迭代运算，以期得到更加对称的切包络线（真是一念之差，绕城十匝）。而且感觉用切包络线代替极点包络线似乎更加合情合理：包络线嘛，就应该尽量将信号包络起来才对，否则包络一部分，留一部分在外面，那就不是完全的包络线了。用切包络线代替极点包络线，最终所得的IMF可能更加合符实际情况，产生虚假IMF的可能性也许会更小，IMF的正交指数也可能会更小。所以我编了求任意信号上下切包络线的程序t angent_e nvelope（实际效果见图35-3），并且将它代替emd程序中的极点包络线进行了迭代运算。虽然现在好几列信号的运行都已经通过，但运行时间需增加几十倍，而且这些信号序列数据都较短（10000点以下）。所以我现在还无法评价它的性能较既有EMD程序的优劣性。绘图看看。所编tangent_envelope函数格式： =tangent_envelope(x) 输出tenvb,tenva,tpoib,tpoia分别是分析信号x的下、上切包络线；下、上切点。tangent_envelope内部未对x作延拓处理。下面是IMF上下切包络线某局部接近0值时的情况： x=imfMBrd(4,:); =tangent_envelope(x); = extr(x);% lx=length(x); zs=zeros(1,lx); fg1=figure(1);% set(fg1,'position', ) set(gca,'position', ); plot(x,'r-') hold on% plot(tenvb,'b--')%下切包络线 plot(tenva,'b--')%上切包络线 plot(zs,'k-.')% plot(tpoib,x(tpoib),'k.')%下切点 plot(tpoia,x(tpoia),'k.')%上切点 plot(indmin,x(indmin),'ko')%极小值点 plot(indmax,x(indmax),'ko')%极大值点 xlim( ) title('imfMBrd的4分量上下切包络线及其切点、极值点') 运行，得：图35-1 图35-2 可见IMF切包络线在接近“0”时，其对称性不总是太好。图中是选了比较极端的情况，一般情况下其对称性当然没有这么差。将上面程序中的x=imfMBrd(4,:)换成任意信号x=MB2917_rd，修改一下绘图句柄、x轴显示范围，运行，得：图35-3 从上图来看，在EMD中用“切包络线”代替“极点包络线”进行求取IMF的迭代运算，似乎更加“合情合理”。上面两图只是看看数据中部局部情况，故分析信号未作延拓。下面看看求整个IMF“切包络线综合振幅”的情况： x=imfMBrd(13,:); = signal_extension(x); lxr= ; ltr= ; abslxr=abs(lxr); absxl=abs(xl); absx=abs(x); absxr=abs(xr); = extr(abslxr); =tangent_envelope(abslxr); tenval=tenva(1:length(tl)); tenvam=tenva(length(tl)+1:length(tl)+length(t)); tenvar=tenva(length(tl)+length(t)+1:length(tl)+length(t)+length(tr)); llxr = length(lxr); zs=zeros(1,llxr); fg3=figure(3); set(fg3,'position', ) set(gca,'position', ); plot(t,absx) hold on% plot(tl,absxl,'r-')% plot(tr,absxr,'r-')% plot(ltr,tenva,'k--') plot(indmax+tl(1),abslxr(indmax),'ko')%% plot(tpoia+tl(1),abslxr(tpoia),'.k') plot(ltr,zs,'k-.') axis( ) title('imfMBrd之13分量延拓并取绝对值后上切包络线及其切点、极值点') 运行，得：图35-4 “切包络线综合振幅”情况良好。续谈。（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de01016two.html 首发时间：2012-06-22 22:49:20）

个人分类: 斤斤计较|4317 次阅读|2 个评论

34、(四)信号延拓_基于IMF自身的瞬时频率估计

baishp 2012-11-30 15:24

关于IMF原信号与其包络线相割的问题，解决之道当然就是求取与与原信号相切的包络线(以后可称之为“切包络线”)，这是能够做到的。为了同时解决瞬时振幅在信号过零点处的阶梯突变，可以将原信号取绝对值之后，再求其切包络线(只求信号上部一条包络线)。为了解决端点效应，求极点包络线时需要先将信号极值点向两端外面延拓。求切包络线时则需先将整个信号数据向两端外面延拓。这个工作比极点延拓稍微麻烦一点，因为延拓部分既要包含延拓极值点，又要与原信号的端点连接是连续、大致光滑的(至少不要变成新极值点了)。步骤如下： 1、将极值点延拓函数 = boundary_conditions(indmin,indmax,t,x,z,nbsym) 的输出tmin,tmax,zmin,zmax改为： = boundary_conditions_lr(indmin,indmax,t,x,z,nbsym)。上面输出tlmin,tlmax是左延极小值、极大值点时间序号，zlmin,zlmax是相应函数值；trmin,trmax是右延极小值、极大值点时间序号，zrmin,zrmax是相应函数值。将原函数名改一下（boundary_conditions_lr）保存。 2、调用上述函数，并将待分析的IMF原函数左、右两个端点分别“混入”上面左、右极值点中； 3、将包括上述“混入”点在内的左延极小值、极大值点时间序号，统一按从小到大顺序排序，并保持对应的函数值跟随排序。同样原则处理右延部分； 4、将上述左延点作样条插值，输出的插值函数时间区段，为左延点最小时间序号点更左一个采样点，至原信号左端点时间序号更左一个采样点（即“0”）。同样原则处理右延部分； 5、令上述所得左延插值函数最左边一点的函数值，等于左边第3点函数值。此举是为了保证左边第2点（原函数极值点延拓而来）仍为极值点。同样原则处理右延部分；信号延拓完成。特别说明一下，本延拓方法适用于任意信号，不局限于IMF。根据上述步骤所编程序，函数文件名存为signal_extension，调用格式为 = signal_extension(x)。其 tl,t,tr分别为左延信号、原信号、右延信号时间序号。 xl,x,xr分别为左延信号、原信号、右延信号函数值。绘图： x=imfMBrd(11,:); = signal_extension(x); lxxx = length(xl)+length(x)+length(xr); zs=zeros(1,lxxx); %下面是为了绘图而重复运行 = extr(x); t=1:length(x); z=x; nbsym=2; = boundary_conditions_lr(indmin,indmax,t,x,z,nbsym); %上面是为了绘图而重复运行 fg1=figure(1); set(fg1,'position', ) set(gca,'position', ); plot(t,x) hold on% plot(tl,xl,'r-') plot(tr,xr,'r-') plot( ,zs,'k--') plot(indmin,x(indmin),'ko') plot(indmax,x(indmax),'ko') plot(tlmin,zlmin,'ko')% plot(tlmax,zlmax,'ko')% plot(trmin,zrmin,'ko')% plot(trmax,zrmax,'ko')% grid title('imfMBrd之11分量左右延拓及其极值点 ') 运行，得：图34-1 将x换成一个IMF高频分量再看看： x=imfMBrd(3,:); = signal_extension(x); lxxx = length(xl)+length(x)+length(xr); zs=zeros(1,lxxx); %下面是为了绘图而重复运行 = extr(x); t=1:length(x); z=x; nbsym=2; = boundary_conditions_lr(indmin,indmax,t,x,z,nbsym); %上面是为了绘图而重复运行 fg2=figure(2); set(fg2,'position', ) li=1; lj=2; for i=1:li for j=1:lj n=(i-1)*lj+j; s(n)=subplot(li,lj,n); set(s(n),'position', ); plot(t,x) hold on% plot(tl,xl,'r-')% plot(tr,xr,'r-')% plot( ,zs,'k--')% plot(indmin,x(indmin),'ko')% plot(indmax,x(indmax),'ko')% plot(tlmin,zlmin,'ko')% plot(tlmax,zlmax,'ko')% plot(trmin,zrmin,'ko')% plot(trmax,zrmax,'ko')% grid if j==1 xlim( ) title( )% else j==2 xlim( ) title( )% end end end 运行，得：图34-2 将x换成一个非IMF分量MB2917_rd看看。只将上图程序绘图句柄及时轴显示区限xlim改一下就可以了。运行，得：图34-3 信号边界延拓完成。后篇续谈。（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de01016qx7.html 首发时间：2012-06-19 22:35:46）

个人分类: 斤斤计较|2883 次阅读|0 个评论

33、(三)综合极点_基于IMF自身的瞬时频率估计

baishp 2012-11-30 14:59

将脉搏序列基本模式函数imfMBrd第10分量的瞬时振幅、瞬时频率、复信号虚部再作几张放大、对比图：（1）IMF原信号及其复信号虚部： x=imfMBrd(10,:); = extr(x);%求极值点、 = cenvelope(x,1);%求包络线 envb=env(1,:);%将负包络线提取出来 enva=env(2,:);%将正包络线提取出来 lx=length(x); zs=zeros(1,lx); = afyimf( x ); fg1=figure(1); set(fg1,'position', ) set(gca,'position', ); plot(x,'r-') hold on plot(envb,'k--') plot(enva,'k--') plot(envmoy,'k:') plot(y,'b-') plot(zs,'k-') grid title('imfMBrd 10分量及其复信号之虚部 ') xlim( ) 运行，得：图33-1 fg2=figure(2); set(fg2,'position', ) set(gca,'position', ); plot(x,'r-') hold on plot(envb,'k--') plot(enva,'k--') plot(f*150,'b-')%f*150是为了放大好对比观察 plot(zs,'k-') grid title('imfMBrd 10分量瞬时频率突变位置分析') xlim( ) 运行，得：图33-2 fg3=figure(3); set(fg3,'position', ) set(gca,'position', ); tt=6500:7500; stem3(tt,f(tt),a(tt));%hold on view(-15,30) grid on title('imfMBrd 10分量局部三维时频图') 运行，得：图33-3 从上面几张图，可以更加清楚地看到，虽然没有了以前所说的瞬时频率跳变，但所求得的瞬时振幅、瞬时频率、复信号虚部，都还存在着突变，没有完全达到我希望的结果。其中瞬时振幅的求取方式，是造成问题的根源。因此我打算针对其过零点与割点问题，采取两个改进措施。下面先看第一个措施。此前的瞬时振幅，是由IMF上下包络线取绝对值后，在信号过零点处分段取值、连接而成。除了造成瞬时振幅自身在信号过零点处的突变，也造成瞬时频率与复信号虚部在过零点处突变。我已经想到了一个很好的方法解决这个问题，就是先将正负极值点统一取绝对值后，再作样条插值（只得到一个插值函数）。具体如下： lx=length(x); t=1:lx; z=x; nbsym=2; INTERP='spline'; = boundary_conditions_emd;(indmin,indmax,t,x,z,nbsym); zmin=abs(zmin); TZ= ; TZ=TZ'; TZTZ=sortrows(TZ,1); TZTZ=TZTZ'; tmm=TZTZ(1,:); zmm=TZTZ(2,:); amp = interp1(tmm,zmm,t,INTERP);%振幅此振幅综合由上下极值点插值而成，因此可以叫它“综合振幅”。它也是原imf信号的上包络线，取负值后又是原信号的下包络线，可以统称为“综合包络线”。将上述amp的计算方法与计算结果代替函数 = afyimf( x )中a的计算方法与计算结果，再重复本文前面的计算过程： x=imfMBrd(10,:); = extr(x);%求极值点、 = cenvelope(x,1);%求包络线 envb=env(1,:);%将负包络线提取出来 enva=env(2,:);%将正包络线提取出来 lx=length(x); zs=zeros(1,lx); = afyimf( x ); fg4=figure(4); set(fg4,'position', ) set(gca,'position', ); plot(x,'r-') hold on plot(envb,'k--') plot(-a,'k:')%综合包络线 plot(enva,'k--') plot(a,'k:')%综合包络线 plot(envmoy,'k--') plot(y,'b-') plot(zs,'k-') grid title('综合包络线 imfMBrd 10分量及其复信号之虚部') xlim( ) 运行，得：图33-4 fg5=figure(5); set(fg5,'position', ) set(gca,'position', ); plot(x,'r-') hold on plot(envb,'k--') plot(-a,'k:')%综合包络线 plot(enva,'k--') plot(a,'k:')%综合包络线 plot(f*150,'b-') plot(zs,'k-') grid title('综合包络线 imfMBrd 10分量瞬时频率突变位置分析') xlim( ) 运行，得：图33-5 fg6=figure(6); set(fg6,'position', ) set(gca,'position', ); tt=6500:7500; stem3(tt,f(tt),a(tt));%hold on view(-15,30) grid on title('综合包络线 imfMBrd 10分量局部三维时频图') 运行，得：图33-6 综合振幅自身在过零点处、对瞬时频率、复信号虚部在过零点处的突变的消除作用，图中已经表达的很清楚了，不再赘述。关于信号与包络线相割的问题，后面继续研究。（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de01016odk.html 首发时间：2012-06-17 12:08:03）

个人分类: 斤斤计较|3202 次阅读|0 个评论

31、(一)根本思路_基于IMF自身的瞬时频率估计

baishp 2012-11-30 14:37

前文说了我对学界既有的瞬时频率估计方法的失望，虽然我并没有将所有的那些方法都试遍，但我已经没有信心与耐心仔细去研究它们了，以后有时间、有心情时再去慢慢详细了解吧。现在我打算抛开一切“高深的”理论，从EMD当中所得到的基本物理量，根据瞬时频率最原始的含义，即各时间点上信号源转动的角度，来估计瞬时频率。众所周知，EMD之所以叫EMD（经验模式分解），是因为它最初并不是基于一种理论推导，而是基于这种分解过程所具有的物理意义被世人所公认，而每一组分解结果都能够经受理论的检验(各IMF之间的正交性)。因此上面我打算采用的瞬时频率估计法，是EMD的一以贯之、很自然的延伸。姑且可称之为“基本模式估计法” 这种方法的关键，在于要明白、承认，IMF的包络线，其实就是信号的瞬时振幅(振幅函数)；而IMF是瞬时振幅在实轴上的投影；某时刻IMF值与IMF包络线值之比，就是该时刻瞬时相位的余弦值。有了瞬时相位，通过相位差分即得到瞬时频率；通过振幅、相位，还可以求出振幅在虚轴上的投影，从而得到IMF对应的复数函数。 IMF的包络线有两条，哪一条是IMF的瞬时振幅呢？从理论上来说，IMF的两条包络线的均值为0，两条包络线是完全对称的，因此两条包络线的绝对值是相等的，都等于IMF的瞬时振幅。但在实际处理时，由于两条包络线不可能做到绝对对称，所以还是应该分别处理。在复平面上，当信号源运行在1、4象限时，IMF由上包络线在实轴上的投影形成；当信号源运行在2、3象限时，IMF由下包络线在实轴上的投影形成。因此，瞬时振幅应该根据IMF的情况由上下包络线分段组成。下面就根据这个原则来求IMF的瞬时振幅、瞬时频率等。下面imfMBrd含义见前两篇博文。 x=imfMBrd(11,:);%先任取imfMBrd的一个分量 lx=length(x); zs=zeros(1,lx); = extr(x);%求极值点、 = cenvelope(x,1);%求包络线及其均值 envb=env(1,:);%将负包络线提取出来 enva=env(2,:);%将正包络线提取出来 fg1=figure(1); set(fg1,'position', ); w=0.04; W=0.94; h=0.04; d=0.003; h1=h+0.75-3*d; h2=h+0.50-2*d; h3=h+0.25-d; h4=h; H=0.18; sb1=subplot(4,1,1); set(sb1,'position', ); plot(x,'r-') hold on plot(envb,'g--') plot(enva,'b--') plot(envmoy,'m--') plot(zs,'k-') grid title('分析信号及其包络线') xlim( ) c41=(x0); c23=(x0); dx=diff(x); dx1= ; dx2= ; dx=(dx1+dx2)/2; c12=(dx0); c34=(dx0); c1=c41.*c12; c2=c23.*c12; c3=c23.*c34; c4=c41.*c34; sb2=subplot(4,1,2); set(sb2,'position', ); stem(c1,'.b-') hold on stem(c2,'.r-') stem(c3,'.g-') stem(c4,'.m-') grid title('四象限掩码') xlim( ) x1=x.*c1; x2=x.*c2; x3=x.*c3; x4=x.*c4; sb3=subplot(4,1,3); set(sb3,'position', ); stem(x1,'.b-') hold on stem(x2,'.r-') stem(x3,'.g-') stem(x4,'.m-') grid title('经掩码作用后的分析信号') xlim( ) aenvb=abs(envb); aenva=abs(enva); tampb=c23.*aenvb; tampa=c41.*aenva; tamp=tampa+tampb;%瞬时振幅 sb4=subplot(4,1,4); set(sb4,'position', ); stem(tampb,'.g-') hold on stem(tampa,'.b-') plot(tamp,'k-') grid title('经掩码作用后的瞬时振幅') xlim( ) 运行，得：图31-1 下面求瞬时相位、瞬时频率与复信号虚部。运行，得：图31-2 从图31-2的第3个子图可以看到，所估计的“瞬时频率”没有前几篇博文中所提到的“跳变”。图中也有一些类似“线谱”一样，或者说象两根“门柱”一样高高耸立的频率成分，但它们是处在同一数量级下的变化，与以前所讲的由接近“0Hz”的负频变成接近“0.5Hz”的正频，性质截然不同的。为了与以前的“频率跳变”相区别，我就将此处的频率变化称之为“频率突变”吧。此“频率突变”是由于分析信号与它的包络线相割造成的。包络线以外部分，由于信号、包络线之比值绝对值大于1，反余弦无意义，所以被强制等于1与-1，反余弦值(瞬时相位)为0与0.5*2*pi，两割点之间部分，瞬时频率等于0；两割点之外靠近割点处，瞬时频率即形成“突变”。虽然“频率突变”是频率在同一数量级下的变化，但它毕竟还是人为造成的，不是信号本身的属性，因此也是不完全符合我预期的目标的。我希望的瞬时频率估计值，更加连续光滑。这一篇就暂到这里，相关问题后面继续分析。（本文首发于： http://blog.sina.com.cn/s/blog_6ad0d3de01016f0a.html 首发时间：2012-06-10 19:35:44）

个人分类: 斤斤计较|3702 次阅读|0 个评论

Origin如何求振荡曲线的包络线

热度 2 microtao 2011-10-27 14:47

Origin最低版本: 8.1SR0 1. 导入Origin安装目录下的文件 \Samples\Signal Processing\Sine Signal with High Frequency Noise.dat. 2. 选中Worksheet中两列，菜单选择 Plot: Line: Line 作一直线图 3. 菜单选择 Analysis: Signal Processing: Envelope 出现X-Function对话框. 4 . 在 Envelope type 下拉菜单中选择 Both Envelopes . 设置 Smooth Points 为1. 5. 按 OK 执行. 图形会出现上下的包络线. 可得到如下结果

个人分类: OriginPro|37812 次阅读|4 个评论

帐号		自动登录	找回密码
密码			注册

关闭安全验证

标签: 包络线

相关帖子

相关日志

关闭 安全验证

标签: 包络线

相关帖子

相关日志

关闭安全验证