A Sliding Goertzel Algorithm

5.9 Goertsel 算法及其应用5.9.1 Goertsel 算法如前所述, N 点时域序列)(n x 的离散付里叶变换式为∑=-=-10)()(N n knN W n x k X ， 1,....,2,1,0-=N k （5.9.1）这N 点频域序列是同时被算出的，不可能只计算其中某一个或几个指定点。

Goetzel 算法是为了解决这个问题而提出的。

这个算法把离散付里叶变换看作一组滤波器，将输入端的时域序列与其中一个滤波器的冲激响应序列进行卷积运算，求滤波器的输出序列，即得)(k X 序列的一点。

这种算法利用旋转因子kN W 的周期性，使DFT 运算化为线性滤波运算。

由于1))(2(==-π--kN Nj kN NeW故式（5.9.1）可化为 )(110)()()(m N k N N m km N N m kNNW m x Wm x Wk X ---=-=-∑=∑=， 1,....,2,1,0-=N k （5.9.2） 定义序列)(n y k 为)(10)()(m n k NN m k W m x n y ---=∑= （5.9.3A ）可见)(n y k 是由两个序列卷积而得到的序列。

)()()(n h n x n y k k ⊗= （5.9.3B ） 其中，)(n x 是输入的N 点序列，另一个序列被看作滤波器的冲激响应序列)()(n u W n h knN k -=。

（5.9.4）对比式(5.9.2）和式(5.9.3），可知：按式（5.9.3A ）进行卷积运算，当N n =时，滤波器的输出)(N y k 就是)(k X ：N n k n y k X ==|)()( （5.9.5）对式(5.9.4)进行Z 变换，可得滤波器的系统函数111)(---=zWz H k k （5.9.6）这是一个一阶系统。

图5.9.1示出这个系统的信号流图，相应的差分方程为)()1()(n x n y W n y k kN k +-=-， 0)1(=-y (5.9.7)按照此式进行递推运算，到了 N n = 时刻，即可依据式（5.9.5）得到)(k X 。

一种用于数字锁相放大器的改进Goertzel算法赵宇红; 王帝; 王艳蓉; 闫江【期刊名称】《《电子世界》》【年(卷),期】2019(000)023【总页数】4页(P5-7,13)【作者】赵宇红; 王帝; 王艳蓉; 闫江【作者单位】北方工业大学【正文语种】中文传统数字锁相放大器必需独立生成参考信号序列，且计算量很大，系统复杂。

针对该问题提出了一种改良的Goertzel算法，以完成对采样序列的下变频与降采样，在现有的基于傅里叶变换的数字锁相放大器基础上，引入迭代的相位校正因子，使理论检测范围扩展至任意频率，可以精确的进行对信号相位的跟踪，并对低信噪比信号的幅值与相位进行准确识别。

1 引言锁相放大器(Lock-In Amplifier，LIA)是微弱信号检测领域极其重要的一类仪器，传统数字正交锁相算法除了对采样序列进行解调和低通滤波运算外，还需要使用数字振荡器生成独立参考信号序列。

常见数字振荡器算法有正弦查表法与正交迭代法等，这些算法均需要消耗额外的运算与存储空间，极大增加系统复杂度，也不利于扩展至多频率/多通道检测，难以满足现代信号检测系统需求。

因此，针对数字锁相放大器的简化算法研究正在被大量开展（赵俊杰，数字锁相放大器的实现研究：大连理工大学，2011；Zhang S,Li G,Lin L,et al.Optimization of a digital lock-in algorithm with a square-wave reference for frequency-divided multi-channel sensor signal detection:Review of Scientific Instruments,2016；蒋国亮，杨晖，李然，基于高调制频率的快速数字锁相算法：仪表技术与传感器，2015；Aguirre J,Garcia-Romeo D,Medrano N,et al.Square-Signal-Based Algorithm for Analog Lock-In Amplifiers:IEEE Transactions on Industrial Electronics,2014）。

A regularizing L-curve Lanczos method forunderdetermined linear systemsS.Morigi,F.Sgallari *,1Department of Mathematics,CIRAM,University of Bologna,Piazza di Porta San Donato 5,I-40127Bologna,ItalyAbstractMany real applications give rise to the solution of underdetermined linear systems of equations with a very ill conditioned matrix A ,whose dimensions are so large as to make solution by direct methods impractical or infeasible.Image reconstruction from projections is a well-known example of such systems.In order to facilitate the com-putation of a meaningful approximate solution,we regularize the linear system,i.e.,we replace it by a nearby system that is better conditioned.The amount of regularization is determined by a regularization parameter.Its optimal value is,in most applications,not known a priori.A well-known method to determine it is given by the L-curve approach.We present an iterative method based on the Lanczos algorithm for inexpensively evaluating an approximation of the points on the L-curve and then determine the value of the optimal regularization parameter which lets us compute an approximate solution of the regularized system of equations.Ó2001Elsevier Science Inc.All rights reserved.Keywords:Ill posed problems;Regularization;L-curve criterion;Lanczos algorithm1.IntroductionIn this article,we describe a new regularizing iterative method for the so-lution of ill conditioned underdetermined linear systems of equationsA x b ;A P R m Ân ;x P R n ;b P R m ;m <n : 1/locate/amcApplied Mathematics and Computation 121(2001)55±73*Corresponding author.E-mail address:sgallari@dm.unibo.it (F.Sgallari).1This research was supported by the University of Bologna,funds for research topics.0096-3003/01/$-see front matter Ó2001Elsevier Science Inc.All rights reserved.PII:S 0096-3003(99)00262-356S.Morigi,F.Sgallari/put.121(2001)55±73Linear systems of this kind arise,for example,in the reconstruction of images from projections.In this application,each column corresponds to one pixel of the image to be reconstructed,and each rowto a``ray''from an emitter to a detector.The matrices that arise are typically large and have many more col-umns than rows.Further,they are generally very ill conditioned,because the linear system is obtained by discretizing an ill posed problem,the inversion of a Radon transform.The ill conditioning makes a solution of(1)very sensitive to perturbations in the right-hand side vector.Such perturbations arise,e.g.,from measurement errors[1,2].In order to be able to compute a meaningful approximate solution of(1),the system has to be regularized,i.e.,the linear system(1)has to be replaced by a linear system of equations with a not very ill conditioned matrix,such that the unique solution of this system yields an acceptable approximate solution of(1). We introduce the matrixB: AA T:The solution of the minimal norm of(1)can be written as x: A T y,where y P R m solves the linear systemB y b: 2 We regularize this system by adding a multiple of the identity matrix to B and obtainB k I y b; 3 where k>0is a regularization parameter.We denote the solution of(3)by y k. The smaller the value of k,the closer the solution y k of(3)is to the minimal norm solution of(2),and the more ill conditioned the matrix B k I.We wish to®nd a value of k in(3),such that the matrix B k I is not very ill conditioned, and the solution of(3)is an acceptable approximate solution of(2).Generally, it is not known a priori for which values of k these requirements are satis®ed. The approximate solution x k A T y k of(1),obtained by®rst solving(3)for y k,can also be computed by standard Tikhonov regularization of system(1). This approach,however,as has been shown in[3],requires more computer storage and arithmetic work when the linear system has many more columns than rows.Let us consider the minimization problemfk B yÀb k2 k k A T y k2g; 4 miny P R mwhere kÁk denotes the Euclidean norm here and in the sequel.The solution y k of(4)satis®es the equationB T B k AA T y B T b: 5S.Morigi,F.Sgallari/put.121(2001)55±7357 The comparison of(3)±(5)shows the equivalence of problem(3)and(5) when B is a nonsingular matrix and this justi®es what follows.Let components of the right-hand side vector b be contaminated by errors, e.g.,due to inaccurate measurements,and b exact and e denote the unperturbed right-hand side and the error,respectively.If the norm of the error e is known, then the value of k can be chosen so that the associated solution satis®es the Morozov discrepancy principle[3].However,in many applications the norm of e is not explicitly known,and other approaches to determine a value of the regularization parameter k have to be employed.We will show that the plot of the curve(k B y kÀb k;k A T y k k)is shaped roughly like an``L'',and the optimal value of the regularization parameter k corresponds to the``vertex''point of the L-curve.This idea of L-curve was observed by Lawson and Hanson[4]to determine the regularizing parameter k in the Tikhonov approach,while the choice of the optimal regularization pa-rameter as the vertex point was suggested recently in[5,6],and successively considered in[7].The computation of the L-curve is quite costly for large problems;in fact, the determination of a point on the L-curve requires that both the norm of the regularized approximate solution and the norm of the corresponding residual vector be available.Therefore,usually only a fewpoints on the L-curve are computed,and these values,which we will refer to as the discrete L-curve,are used to determine a value of the regularization parameter.In[8],a newmethod for approximately determining the points on the L-curve is proposed.They showed how rectangular regions that contain points on the L-curve(L-ribbon)can be determined fairly inexpensively for several values of the regularization parameter without having to solve the corresponding regularized minimization problem but using Gauss and Gauss±Radau quadr-ature rules.When the L-ribbon is narrow,it is possible to infer the approximate location of the vertex of the L-curve from the shape of the L-ribbon.In order to solve problem(4),or equivalently(3),we cannot follow the Calvetti et al.[8]approach because their assumptions are not satis®ed due to the term A T in the regularizing functional.We present an iterative procedure that determines successive encapsulated intervals containing the optimal regularization parameter corresponding to the vertex of the L-curve.At each step,some newpoints on the L-curve are computed by a fewsteps of the Lanczos bidiagonalization algorithm,while a new interval for k is de-termined by a minimum distance from the local origin criterion.The®nal value obtained for the regularization parameter is used to compute an approximation of the solution of the regularized system(3).This article is organized as follows.In Section2,we give a brief descrip-tion of properties of the L-curve and we show how to use it to solve prob-lem(4).Section3explains the relation between the bidiagonalization Lanczosalgorithm and the approximated L-curve and howto update the information extracted from the Lanczos process to compute a newpoint.Section4de-scribes the regularizing algorithm for the solution of the linear system(3). Section5contains several computed examples which illustrate the performance of the algorithm.2.L-curve for the choice of regularization parameter valueA recent method to choose the regularization parameter is the so-called L-curve.The L-curve is a parametric plot of q k ;g k ,where g k and q k measure the size of the regularized solution and the corresponding residual. The underlying idea is that a good method for choosing the regularization parameter for discrete ill posed problems must incorporate information about the solution size in addition to using information about the residual size.This is indeed quite natural,because we are seeking a fair balance in keeping both of these values small.The L-curve has a distinct L-shaped corner located exactly where the solution y k changes in nature from being dominated by regulariza-tion errors(i.e.,by over smoothing)to being dominated by the errors in the right-hand side.Hence,the corner of the L-curve corresponds to a good bal-ance between minimization of the two quantities,and the corresponding reg-ularization parameter k is a good choice.In particular the continuous L-curve for our problem consists of all the points q k ;g k for k P 0;I ,where q k k B y kÀb k and g k k A T y k k. In what follows,we present some important properties of the continuous L-curve and a method to compute a good approximation of the discrete L-curve. We®rst need to introduce the following singular value decomposition for matrix A U R V T.Here U P R mÂm and V P R nÂn are orthogonal matrices, R R m0 P R mÂn,and R m diag f r i g,where the singular values r i are ordered as0<r1T r2TÁÁÁT r m.It is straightforward to see that the solution of problem(4)can be written asy kX mi 1u T i bk r2iu i: 6Using these notations,the L-curve has the following properties:Proposition1.Let y k denote the solution to(4).Then,k A T y k k is a monotonically decreasing function of k B y kÀb k.Proof.The fact that the k A T y k k is a decreasing function of k B y kÀb k follows immediately from the following expressions:k A T y k k2X mi 1r i u T i br2i k2; 758S.Morigi,F.Sgallari/put.121(2001)55±73k B y k Àb k 2X m i 1kr 2i k u T ib 2; 8which are easily derived from (6).ÃIf a fewreasonable assumptions are made about the ill posed problem,then it is possible to characterize the L-curve behaviour analogously to that done in [5,6]for the Tikhonov regularization.The three assumptions made are:1.The coe cients j u T i b exact j on average decay to zero faster than the r i ;2.The perturbation vector e is essentially ``white noise''with standard devia-tion r ;3.The norm of e is bounded by the norm of b exact .Under Assumptions 1±3,the L-curve q k ;g k is shown to exhibit a corner behaviouras a function of k ,and the corner appears approximately at r m p ;k A T y 0k .Here,y 0is the unregularized solution to the unperturbed problem (2).The larger the di erence between the decay rates of j u T i b exact j and j u Ti e j ,the more distinct the corner will appear.Following [6],we summarize the reasoning that leads to the previous characterization.We consider the behaviour of the L-curve for an unperturbedproblem with e 0.For k (r 21,w e have r i = r 2i k %1=r i ,so that y k %y 0,andk B y k Àb exact k 2%X m i 1k r 2i u T i b exact 2T k r 21k b exact k 2: 9 Hence,for small values of k ,the L-curve is approximately a horizontal line atk A T y k k k A T y 0k .As k increases,it follows from (7)that k A T y k k starts to de-crease,while k B y k Àb k still grows.Now,consider the L-curve associated with the perturbation e of the right-hand side.The corresponding solution y k given by (6)with b replaced by e ,satis®es (from Assumption 2)k A T y k k 2 X m i 1r i u T i e r 2i k 2%r 2X m i 1r i r 2i k 2; 10 k B y k Àb k 2X m i 1kr 2i k u T ie 2%r 2X m i 1kr 2i k2:11For k (r 21,this L-curve is approximately a horizontal line at k A T y k k %m p r =r 1,and it starts to decrease towards the abscissa axis for much smaller values of k than the L-curve for b exact (see Assumption 3).Moreover,we see that as k increases,k B y k Àb k becomes almost independent of k ,while k A T y k k isS.Morigi,F.Sgallari /put.121(2001)55±7359dominated by a fewterms and becomes about r2=kp.Hence,this L-curve soonbecomes almost a vertical line at k B y kÀb k%rmpas k3I.The actual L-curve for a given problem,with a perturbed right-hand side b b exact e,is a combination of the above two special L-curves.For small k, the behaviour is completely dominated by contributions from e,whereas,for large k it is entirely dominated by contributions from b exact.In between there is a small region where both b exact and e contribute,and this region contains the L-shaped corner of the L-curve.In actual computations,the values q k and g k are usually not available; hence,we search for good approximations.Let y k k be an approximation of the solution y k of(3),and introduce the residual vectorr k k: bÀ B k I y k k 12 and the quantitiesq k k : k B y k kÀb k; 13g k k : k A T y k k k: 14 The following proposition is concerned with how well the functions q k and g k approximate q and g,respectively.Proposition2.Let q k k and g k k be defined by(13)and(14),respectively. Thenj q k k Àq k j T k r k k k; 15 j g k k Àg k j T k r k k k: 16 Proof.It follows from the de®nition of q and q k thatj q k k Àq k j T k B y k kÀy k k k B y k kÀ B k I À1b kk B B k I À1r k k k T k B B k I À1kk r k k k;and in viewof k B B k I À1k T1,the®rst part of the proposition follows.For the second part,we have from the de®nition of g and g kj g k k Àg k j T k A T y k kÀy k k k A T y k kÀ B k I À1b kk A T B k I À1r k k k T k A T B k I À1kk r k k k:Using the singular value decomposition of the matrix A U R V T and B U R2U T,we obtain that k A T B k I À1k T1,then the last part of the proposition follows.Ã60S.Morigi,F.Sgallari/put.121(2001)55±73Thus,we can obtain a good approximation of the discrete L-curve by computing q k k ;g k k ,and quantifying the errors by the length of the re-sidual vector.3.The Lanczos algorithm and regularizationIn this section,we brie¯y review the bidiagonalization Lanczos process[9], and,we explore its relation with the functions q k k and g k k introduced in Section2.Given a symmetric matrix B AA T and a nonvanishing vector b,if we ex-ecute k iterations of the Lanczos bidiagonalization Algorithm1in exact arithmetic,we get the orthogonal mÂk matrixV k v0;v1;...;v kÀ1 17 and the orthogonal nÂk matrixW k w0;w1;...;w kÀ1 : 18 The matrices V k and W k are connected by the two equationsAW k V k L k d k v k e T k; 19A T V k W k L T k; 20 where L k denotes the kÂk lower bidiagonal matrixL k:c0d1c1......d kÀ2c kÀ2d kÀ1c kÀ12666664377777521and e k is the k th column of the identity matrix.Algorithm1(Lanczos bidiagonalization algorithm).Input:b P R m n f0g,B P R mÂm,0<k<m;Output:f d j g k j 1;f c j g kÀ1j 0;f v j g kÀ1j 0;v0: b=k b k;s0: A T v0;c0: k s0k;w0: s0=c0;for j 2;3;...;k doS.Morigi,F.Sgallari/put.121(2001)55±7361r jÀ1: A w jÀ2Àc jÀ2v jÀ2;d jÀ1: k r jÀ1k;v jÀ1: r jÀ1=d jÀ1;s jÀ1: A T v jÀ1Àd jÀ1w jÀ2;c jÀ1: k s jÀ1k;w jÀ1: s jÀ1=c jÀ1;end jLet T k 1;k denote a tridiagonal matrixT k 1;k:a0b1b1a1b2.........b kÀ2a kÀ2b kÀ1b kÀ1a kÀ1b k2666666666437777777775P R k 1 Âk: 22It is nowstraightforw ard to identify T k;k L k L T k if a j d2j c2j, b j d j c jÀ1;j 1;...;kÀ1,and a0 c20.By combining the two relations(19) and(20),we getBV k V k L k L T k d k c kÀ1v k e T k V k 1T k 1;k; 23 with V k 1 v0;v2;...;v k ,such that V T k 1V k 1 I k 1and V k 1e1 b=k b k. Throughout this article,I k 1denotes the identity matrix of order k 1,I k 1;k its leading principal k 1 Âk submatrix.It follows from Proposition2that in order to have a good approximation of the L-curve points q k ;g k ,we have to choose k su ciently large and in-crease it until k r k k k is under a given tolerance.Now,we want to show how,by means of Algorithm1,it is possible to compute the residual k r k k k for increasing values of k in a straightforward manner.For notational simplicity,we will consider in the sequel the T k 1;k matrix even if,from a computational point of view,only the necessary elements are determined starting from L k.We determine an approximate solution of(3)of the form y k k V k z by solving the minimization problemmin z P R k k B k I V k zÀb k minz P R kk V k 1 T k 1;k k I k 1;k zÀV k 1e1k b kkminz P R kk T k 1;k k I k 1;k zÀe1k b kk: 24Let us introduce the QR-factorization of the tridiagonal matrix~Tk 1;k : T k 1;k k I k 1;k ~Q k 1~R k 1;k; 2562S.Morigi,F.Sgallari/put.121(2001)55±73where ~Q k 1P R k 1 Â k 1 ,~Q T k 1~Q k 1 I k 1and ~R k 1;k P R k 1 Âk has an upper triangular leading principal k Âk submatrix,denoted by ~R k ,and a vanishinglast row.It follows thatmin z P R kk B k I V k z Àb k min z P R kk ~R k 1;k z À~Q T k 1e 1k b kk j e T k 1~Q T k 1;k e 1jk b khas the solutionz k k: ~RÀ1kI T k 1;k ~Q Tk 1e 1k b k :26Thus,the associated approximate solution of (3)is given byy k k : V k ~R À1k I T k 1;k ~Q T k 1e 1k b k ;27and the corresponding residual vector (12)has the normk r k k k j e T k 1~Q Tk 1e 1jk b k :28We note for future reference that,when y k k is given by (27),the quantities q kand g k de®ned by (13)and (14),respectively,can be evaluated asq k k k T k 1;k z k k Àe 1k b kk ;29 g k k k L T k z k k k ;30where z k k is given by (26).When increasing k ,the QR-factorizations of the matrices ~Tk 1;k are updated.We outline the computations required for this.As the matrix ~Tk 1;k is tridiag-onal,it is convenient to compute its QR-factorization by Givens rotations [10,Chapter 5].Let G ik 1be a Givens rotation of order k 1that rotates the i ;i 1 coordinate planes so that the i th subdiagonal element of the matrix G i À1 k 1;...;G 1 k 1~T k 1;k vanishes.Then,the factors in the QR-factorization (25)of ~Tk 1;k are given by ~Q k 1: G k k 1;...;G 1 k 1T ;~R k 1;k : G k k 1;...;G 1 k 1~T k 1;k :31Taking an additional step with the Lanczos algorithm yields the matrixwhose QR-factorization can be determined by updating factors (31)as follows.For 1T i T k ,let G ik 2denote the Givens rotation of order k 2whose leadingprincipal submatrix of order k 1is G ik 1.Then,32where Ãdenotes a matrix element that may be nonvanishing.The upper tri-angularmatrixin the QR-factorization ~Tk 2;k 1 ~Q k 2~R k 2;k 1is obtained by multiplying matrix (32)from the left by a Givens rotation G k 1k 2that annihilates the k 2;k 1 entry b k 2.Thus,andQ T k 2e 1 G k 1 k 2~Q T k 1e 102435:64S.Morigi,F.Sgallari /put.121(2001)55±73We point out that in order to compute ~Rk 2;k 1,given ~R k 1;k ,we only need to apply the Givens rotations G k À1 k 2,G k k 2and G k 1k 2,in this order,to the lastcolumn of ~Tk 2;k 1.We note that q k k and g k k ,for ®xed k ,can be evaluated for di erent values of k without recomputing the Lanczos decomposition (23).Each newvalue of k yields a newtridiagonal matrix ~Tk 1;k ,whose QR-factorization,see (25)and (31),has to be computed.The computation of this factorization is inexpensive;it requires only O k arithmetic operations,as it can be carried with k Givens rotations.4.A regularizing algorithm based on the Lanczos processIn this section,our method for the solution of underdetermined ill condi-tioned linear systems (1)is described.The evaluation of points on the L-curve requires that the exact solution of the linear system (3)is available.However,as already mentioned,computation of the exact or a highly accurate solution of (3)for each value of k is expensive,and therefore,we replace the equations for q k and g k by (29)and (30),respectively.It follows from Proposition 2and (28)that the norm of the re-sidual vector r k k is a computable inexpensive bound for the error introduced when approximating q k and g k by q k k and g k k .Algorithm 2computes the regularization parameter and the corresponding approximate solution of the regularized system (3).Algorithm 2(Regularizing L-curve Lanczos iteration ).Input :A P R m Ân ,b P R m , d , r P R ,[k i ,k f ],k i ;k f P R ,N ;k val P N ;Output :k opt ,regularized approximate solution x kopt k P R n ;%Initialization %.k : k val;j 0;k 1i k i ,k 1f k f ;k 0opt : 0;repeat.j j 1;.k j s : k j f Àk ji = N À1 ,%L-curve computation phase %.Perform k Lanczos steps by applying Algorithm 1;for k `: k j i to k j f step k js do.Compute QR-factorization of ~Tk 1;k ;.Evaluate k r k k ;while (k r k k > r )do .k k 1;.Apply one more Lanczos step;end ;S.Morigi,F.Sgallari /put.121(2001)55±736566S.Morigi,F.Sgallari/put.121(2001)55±73 .Compute point P` q k k` ;g k k` ;end;%Detection corner phase%.Detect the smallest interval that contains the closest point P j to the origin; .Compute the associated k j opt;until(j k jÀ1optÀk j opt j T d j k jÀ1opt j);%Resolution phase%.Evaluate the®nal solution of A x b by solving x k opt k A T y k opt k.The initialization step needs the de®nition of an initial interval k i;k f for the regularizing parameter,the number of required points(N)on the discrete L-curve,and the initial dimension k val of the Krylov space.We will refer to each step j of the repeat loop as re®nement step j.Each of these steps consists of an L-curve computation phase,which calculates N even spaced points on the discrete L-curve belonging to the interval k j i;k j f ,and of a detection phase for determining the corner of the L-curve at step j.Let us consider the L-curve computation phase in more detail.Given an initial k value(k val),k Lanczos steps are preliminarily performed. Using the relations between d and c entries of L k matrix and a and b entries of T k 1;k the~Q~R factorization(25)is then computed and further Lanczos steps are performed until the residual norm,evaluated by relation(28),does not satisfy the given precision r.The approximate point P` q k k` ;g k k` of the discrete L-curve is then given by(29)and(30).To determine the next point P`of the L-curve the algorithm uses for k and L k the previous values,then,using(25),a newQR-factorization is computed and consequently the residual is available.The corner detection phase considers the L-curve resulting from the previ-ous phase and determines the smallest interval containing the k value corre-sponding to the closest point P j opt to a local origin to a given tolerance d.The local origin coordinates are computed in the®rst re®nement step as min`q k k` and min`g k k` .Note that during each re®nement step,we can also verify whether the given j th interval contains the corner point P j opt.In fact,if the k value cor-responding to the returned point P opt j is equal to one of the extreme points of the interval k j i;k j f ,then the k value associated with the corner point is not contained in the present interval.Thus,it will be necessary to expand the interval and consequently to compute some additional points to complete the L-curve.The last step of Algorithm2evaluates the approximate solution y k opt k of(3) by(27),and the regularized approximate solution of the original problem(1) by x k opt k A T y k opt k.As a®nal remark,we observe that each re®nement step can take advantage of the previous Lanczos step without recomputing all the values.5.Numerical resultsThe purpose of this section is to illustrate the performance of Algorithm2in some numerical examples obtained on a Sun Ultra workstation with a unit roundo of 2À52%2:2204Â10À16.We implemented our method to com-pute an approximation of the L-curve in FORTRAN77language,and we used MATLAB commercial software package[11]to obtain the exact computation of points on the L-curve.In order to obtain small and large scale test problems,we generated an m-by-n matrix A de®ned by its singular value decompositionA U R V T; 33 whereU IÀ2uu Tk u k22;V IÀ2vv Tk v k22; 34and the entries of u and v are random numbers.R is a diagonal matrix with the singular values r i eÀ0:2 iÀm .On the right-hand side of the linear system(1)to be solved b exact is a random vector,and e is the noise vector generated by®rst determining a vector with normally distributed random numbers with zero mean and variance one as components,and then scaling it to have the desired norm.We will refer to the quantity k e k=k b k as the noise level.All computations here reported correspond to small ill conditioned test problems m 100;n 200;r m 1:0;r1 2:5175Â10À9 in order to be able to perform qualitative comparisons with the MATLAB results.The results of each re®nement step of Algorithm2are summarized in Table 1which consider,respectively,the three di erent noise levels0:1;0:01and 0:001.The®rst column in the tables indicates the level of re®nement(LR),the second column labelled`` k i;k f ''shows the interval of the regularization pa-rameter k,the columns labelled``k opt e''and``k opt a''contain the optimal regular-ization parameter determined in that re®nement interval,respectively,by MATLAB(exact)and by Algorithm2(approximated).The last two columns report the distance of the corner from the local origin,chosen for a given noise level(distance),and the number of evaluation points in the re®nement step (NP).The precision`` r'',used in Lanczos procedure,is reported in the captions. The dimension k of the Krylov subspace is not reported in the tables because it maintains a value of®ve for all the computations.The only exception is for some extremely small k values(e.g.,1:0Â10À9;1:0Â10À10)when a high pre-cision r is used in the Lanczos procedure.To avoid the use of k greater than 100,we suggest either to consider a smaller interval for k,or to use,for the computation related to the critical extreme value,a higher precision r.S.Morigi,F.Sgallari/put.121(2001)55±7367Table1Test problemLR k i;k f k opt e k opt a Distance NP(a)Noise level0:1, r 0:11 1:0Â10À9;2:0 3:0Â10À13:0Â10À15:963200Â10À1202 4:0Â10À1;2:0Â10À1 2:5Â10À12:5Â10À15:956456Â10À1203 2:4Â10À1;2:6Â10À1 2:5Â10À12:5Â10À15:956456Â10À1201 1:0Â10À9;2:0 2:5Â10À12:5Â10À15:956456Â10À1200(b)Noise level0:01, r 0:051 1:0Â10À9;2:0 1:0Â10À11:1Â10À1 1.211257202 1:0Â10À9;2:0Â10À1 1:0Â10À21:0Â10À28.901137Â10À1203 1:0Â10À9;2:0Â10À2 1:0Â10À31:0Â10À37.939300Â10À1204 1:0Â10À9;2:0Â10À3 1:3Â10À31:0Â10À37.930785Â10À1201 1:0Â10À9;2:0 1:0Â10À21:0Â10À28.901137Â10À12002 1:0Â10À9;2:0Â10À2 9:0Â10À41:0Â10À37:930785Â10À1200(c)Noise level0:001, r 0:0051 1:0Â10À10;2:0 1:0Â10À11:0Â10À1 2.232340202 1:0Â10À10;2:0Â10À1 1:0Â10À21:0Â10À2 1.715200203 1:0Â10À10;2:0Â10À2 1:0Â10À31:0Â10À3 1.306169204 1:0Â10À10;2:0Â10À3 1:0Â10À41:0Â10À49.563796Â10À1205 1:0Â10À10;2:0Â10À4 1:0Â10À51:0Â10À58.587058Â10À1206 1:0Â10À10;2:0Â10À5 7:0Â10À66:0Â10À68.572574Â10À1207 6:0Â10À6;8:0Â10À6 7:2Â10À65:8Â10À68.570727Â10À1201 1:0Â10À10;2:0 1:0Â10À21:0Â10À29.806748Â10À12002 1:0Â10À10;2:0Â10À2 1:0Â10À41:0Â10À48.646498Â10À12003 1:0Â10À10;2:0Â10À4 5:0Â10À67:0Â10À68.572574Â10À120068S. Morigi, F. Sgallari / Appl. Math. Comput. 121 (2001) 55±73。

《激光杂志》2021年第42卷第2期LASER JOURNAL(Vol.42,No.2,2021)147•光电技术与应用•基于YOLOv5网络模型的人员口罩佩戴实时检测谈世磊J别雄波S卢功林J谈小虎''重庆大学自动化学院，重庆400044；2重庆科技发展战略研究院有限责任公司，重庆401123；'太原科技大学应用科学学院，太原030024摘要:近年来，随着硬件算力的提升和人工智能算法的创新发展，使得深度学习算法在目标检测方面有着广泛的应用。

针对现有人工方式查看人员口罩佩戴情况的不足，提出了一种基于深度学习YOLO v5算法实现对口罩佩戴情况的实时检测。

算法首先将数据集进行归一化处理，再将数据接入YOLO v5网络进行迭代训练，并将最优权重数据保存用作测试集测试，算法通过tensorboard可视化显示训练和测试结果。

实验结果表明，所提算法检测的准确性高，实时性强，满足实际使用需求。

关键词:YOLOv5；口罩检测；最优权重;tensorboard中图分类号:TN249文献标识码：A doi：10.14016/ki.jgzz.2021.02.147Real-time detection for mask-wearing of personnel based onYOLOv5network modelTAN Shilei1, BIE Xiongbo2,LU Gonglin1,TAN Xiaohu31 College of automation,Chongqing University,Chongqing400044,China；2Chongqing Academy of Science and Technology far Development,Chongqing400123,China；3College of Applied Sciences,Taiyuan University of science and technology,Taiyuan030024,ChinaAbstract:In recent years,with the improvement of hardware computing power and the innovation and develop­ment of an artificial intelligence algorithm,a deep learning algorithm has been widely used in target detection.Given the deficiency of the existing manual way to check the masks'wearing,this paper proposes a real-time detection meth・od based on deep learning YOLOv5algorithm.Firstly,the data set is normalised,and then the data is connected to the YOLOv5network for iterative training,and the optimal weight data is saved for the test set test.The algorithm vis­ually displays the training and test results through Tensorboard.The experimental results show that the proposed algo­rithm has high accuracy,strong real-time performance,and is able to meet the practical demands.Key words:YOLOv5；mask detection；optimal weight；Tensorboard1引言近年来，人工智能技术被应用于各行各业中，从现有产品的实际效果来看，计算机视觉和自然语言处理算法的效果显著，优于现有的人工或机器学习方法。

名 词英文名称运筹学Operations research(注：美国称：Operations research；英国称为：Operational research;英文缩写：OR）Linear programming线性规划英文缩写LPLinear programming英文缩写LP线性规划数学模型Mathematical models单纯形法Simplex methodSimplex method单纯形法改进单纯形法Revised Simplex method目标函数Objective约束条件Constraints可行解Feasible solutions可行域Feasible region名 词英文名称线性规划图解法Graphical Solution of Linear Programs对偶理论Duality theory对偶单纯Dual Simplex method形法影子价格Shadow price运输问题transportation problemtransportation problem运输问题目标规划法Goal programming表上作业法Tabular method表上作业法Tabular method图上作业法Graphical method灵敏度分析Sensitivity analysis西北角法Northwest corner rule名 词英文名称最小元素The least cost rule法运输论法Transportation闭回路调整法Close circular adjust method名 词英文名称非线性规划Nonlinear programming斐波那契法Fibonacci search0.618法Golden section search（黄金分割法）欧拉回路Euler loop整数规划Integer programming松弛问题Slack problem名 词英文名称割平面法Cutting plane method分枝限界Branch and bound method法整数线性规划Integer linear programming纯整数线性规划Pure Integer linear programming名 词英文名称混合整数Mixed Integer linear prog.线性规划0-1型整数线性规划Zero-one Integer linear programming 隐枚举法Implicit Enumeration马氏决策Markov decision programming规划最小树问Minimal tree problem题名 词英文名称最短路问Shortest-route problems题Dijkstra算法Dijkstra algorithmFloyd算法Floyd algorithm最大流问Maximal-Flow problems题图与网络分析Graph theory and network analysis 网络计划Network program网络Network网络方法Network method and network planning 和网络计划网络分析Network analysis网络技术Network techniques名 词英文名称关键线路Critical path method 简称CPM法计划评审Program eval-法uation & rev-iew technique简称PERT网络图Network graphic多重图和Multiple graph and简单图simple graph连通图Connected graph无向图Non-oriented graph有向图oriented graph名 词英文名称最短路径问题Shortest path problem动态规划Dynamic programming缩写DP决策Decision决策论Decision theory名 词英文名称现代决策Modern decision theory理论古典决策Classical decision theory理论战略决策Strategy decision风险型决Risk decision策益损矩阵Opportunity loss matrix最大可能Maximal probability criterion法名 词英文名称期望值法Expected value method决策树法Decision trees method局中人Player策略Policy马氏决策Markov decision programming规划英文缩写：MDP悲观准则Max-min criterion（max-min准则）乐观准则Max-max criterion名 词英文名称折衷准则Trade-off criterion等可能准则(Laplace准则)Laplace criterion遗憾准则（min-max准则）Regret criterion对策论Game theory合作对策Cooperative games名 词英文名称非合作对Non-cooperative games 策纳什平衡Nash equilibrium帕雷托最Pareto optimality优斯塔克尔贝格对策Stackelberg strategy统筹法Overall planning method 指派问题Assignment problem匈牙利解Hungarian method法存储论Inventory theory排队论Queuing theory名 词英文名称排队系统Queuing system生灭过程Birth-death process支付矩阵The payoffmatrix内 容运筹学是一门运用于管理有组织系统的科学。

DTMF Decoding with a PIC16xxxAuthor: David HedleyMicrochip Technology Inc.Additional Input:Chris MacCallumDes HowlettMicrochip Europe INTRODUCTIONThis application note describes how to decode standard DTMF tones using the minimum number of external discrete components and a PIC. The two examples use a PIC which has an 8 bit timer and either a comparator or an ADC, although it can be modified for use on a PIC which has only digital I/O. The Appendices have example code for the 16C662 (with comparator) and 16F877 (using the ADC).As the majority of the Digital Signal Processing is done in software, little is required in the way of external signal conditioning. Software techniques are used to model the individual elements of a DTMF Decoder IC.∙Data collection mechanism.∙ A time domain to frequency domain conversion (Fourier Transform)∙ A bank of 8 Narrow Bandpass Filters∙Decision Making∙OutputTHEORYTable 1 : Common DTMF Frequencies.The digits 0 through 9 and the characters * and # are used in public telecommunications services. The characters A through D are reserved for use in non-public telecommunications networks. Although the character D, is sometimes used as an identifier in Caller Identification signalling on the public network.For both working examples, the input signals to the pin of the PIC were 0v to 5v positive waves. Some circuitry may be needed to ensure the correct levels are sent to the device. The aim of this Application is to decode the signals without the need for external filtering equipment, or a DTMF decoder IC, and allow enough remaining processor memory and overhead to incorporate other features. i.e. Call Routing, Identification, FSK Demodulation etc.Within the software, the data is collected a sample at a time and acted upon immediately. The sample rate according to Nyquist should be at least twice the highest frequency in the range. In the case of standard DTMF tones, the highest frequency is 1633Hz, the sample frequency should therefore be greater than 3.3KHz (Approx Rate: 300us). Using Timer0 to create the sample trigger via an Interrupt, means we can use the TMR0 Prescaler to set an appropriate sample rate (256uS / 4KHz is easy to achieve). The input is scanned, from the input mechanism (comparator/ ADC /digital I/O) and is stored in a RAM bit for comparison.The comparison is widely known as a “modified Goertzel algorithm”; a number of inputs of 1 bit resolution are provided (the incoming signal), which must be converted from the time domain into a component of frequency. This operation is essentially a 1 bit Discrete Fourier Transform (DFT). In a PIC, this involves building a frequency response table in RAM, as the samples are taken. Each incoming bit is compared against the expected result for each of the 8 frequencies.The expected values are held in program memory in RETLW tables. These tables were produced from a spreadsheet using theoretical SINE and COSINE waves at the given DTMF frequencies. (File Attached *.xls). Each “match” increments an accumulator variable for that particular frequency. A “non-match” will decrement the same counter.Example 1:Incoming SignalExpected Signal+1+1+1+1+1+1+1+1+1+1+1+1+1+1+1+1 = 16Example 2:Incoming SignalExpected Signal1+1+1-1-1-1-1-1-1+1+1+1+1+1-1-1-1-1-1 = -3We can see that Example 1 yields a large sum, compared to Example 2 which is off frequency. Given more samples, Example 2 will continue to move around zero.We must compare each incoming bit with the expected SINE and COSINE bits for each of the 8 frequencies, a total of 16 comparisons for each incoming bit. The reason for comparing the incoming bit with the SINE and COSINE of each frequency is that even if there is an exact frequency match, the phase of the incoming signal affects the magnitude of the sum of products. In the extreme case where the signal is a match, there arepositions of phase that would cause the SINE accumulator to sum to zero. However, if we check against the COSINE as well, this will give a large magnitude of matches at 90 degrees and the sum of SINE and COSINE results will still prove the match.Table 2: Comparison and Accumulation.The comparison performed in software can be shown manually with the simplified table above, which demonstrates the systems immunity to phase differences. There wouldn’t be enough samples in the system for this to work, but it shows that the result will be the same, regardless of the phase of the incoming signal.To get the SINE Total for Sig1, compare A,SINE (0) against A,Sig1 (1). There was a NO MATCH, so the running total is decremented (now –1) B,SINE (1) against B,Sig1 (1) = MATCH so increment the running total (now 0). C,SINE (1) against C,Sig1 (0) = NO MATCH so decrement the running total (now -1). D,SINE (0) against D,Sig1 (0) = MATCH so increment the running total (now 0). Do the same for the COSINE, then cast the variables to be unsigned eg. –x = +x.Summing the results from the SINE and COSINE comparisons for each of the eight frequencies will begin to build a frequency response table in RAM. This has eliminated the potential phase problem. Over time and given enough samples, one winning frequency will emerge from each of the low-band and high-band frequency ranges. If there are not two clear winners, no output will be given.Table 3: Frequency Response Table in RAM.The process of decision making is straightforward and requires the detection of the winning frequencies, and confirmation that the winners are above a set threshold.In the case above, button “*” was pressed, and the output from the PIC can be made to reflect this via LED’s or a serial communication message.The example code for the PIC16c662 shows the output on PORTB, each bit being one of the eight frequencies. A pair will be lit at the end of the conversion if a DTMF tone was detected. The example code for the PIC16F877 sends the ASCII character representationof the dialled button down an RS232 link to a terminating display (in this case HyperTerminal on the PC works well - 9600,8,N,1), using the PIC’s USART at9600baud.IMPLEMENTATIONDuring initialisation, (INIT), user RAM is cleared using indirect addressing. The File Select Register (FSR) is incremented and the pointer contents (INDF) cleared. The loop runs to the end of the specified RAM area. In this case the end is at 0x4F. It’s a good way to ensure your variables are clear at the beginning of your conversion, and only uses 8 lines of code.The input peripheral (comparator or ADC) is set up, as are the timers which need to be used. As the sample time in this application is 256us, the OPTION register is used to ensure an accurate scaling for Timer 0, and the interrupt for the timer peripheral is enabled.The MAIN routine is a simple loop. This is done so that other functions can be performed if the desired application is more complex than the simple decoder shown here. Add your code here, but be aware that during the time spent in this area, the PIC is waiting for Timer interrupts to occur. The ISR saves and restores the W and STATUS registers, but may need to include other registers if more complex code is added. Once the interrupt routine has been executed and a sample taken, a user flag is detected on returning to theMAIN routine. This is the cue to perform the Goertzel comparison. The device runs with a 16MHz oscillator mainly because the comparison needs to be completed before the next interrupt.The routine labelled GOERTZEL has two bytes and one bit passed to it. The two bytes are the expected SINE and COSINE values for that sample, the single bit is the sampled signal level, and is held in a variable called BYTE, bit 0. This level is tested against the other 16 in turn, using rotations, and the 16 variables dedicated to the accumulators are incremented or decremented accordingly. These accumulators can be visualised as the frequency response table. Once all the samples have been acted upon, the SINE accumulator for each frequency is summed unsigned with its COSINE counterpart.This is performed in a routine labelled ACCUM_LOW. Again the FSR is used extensively to address the RAM locations and keep the code to a minimum. Given that the sum totals, now held in the COSINE variable for each frequency can be a positive or a negative, they need to be made positive so that they can be easily compared during the output stage. As there are never more than 128 samples taken in this application, the most significant bit can be used to determine the sign, and if necessary, “flipped” positive around zero.The output uses the FSR again to decide if there was a pair of winning frequencies, and if the winners were above a certain threshold level. The threshold level for a certain match using the full 128 byte tables was set at 0x40 (Decimal 64). This means that the sum ofthe SINE and COSINE comparison must be above a specific level for the match to be proven. In theory, the maximum matched sum (in the case of square waves) would be Decimal 128, so a threshold of half the maximum works well. Equally, if the half tables are used to reduce the conversion time, the threshold should be altered accordingly. It should be noted that this makes the system more susceptible to errors, and the application has not been fully tested in this way, at this time.In order to use the “half-tables” approach, the end point for the counter needs to be changed in the Interrupt Service Routine (ISR). The variable TABLECOUNTER is used to determine the offset used when the table is called. This counter is checked once the table values have been returned, and if the specified number of passes have been completed, a flag is raised to tell the device to finish the process and output the results. Changing the tested value from Decimal 128 to 64, will have the desired effect.TIMING AND STATISTICSNote: All timings taken using ICE2000 emulator system at 16MHz.Using FULL TABLES (128 samples)Theoretical minimum (128 samples * 256uS) - 32.51mSPIC16F877 example:Total conversion time (including output stage) - 32.75mSPIC16C662 example:Total conversion time (including output stage) - 32.74mSUsing HALF TABLES (64 samples)Theoretical minimum (64 samples * 256uS) - 16.384mSPIC16F877 example:Total conversion time (including output stage) - 16.6mSPIC16C662 example:Total conversion time (including output stage) - 16.6mSThe files included are assembled with MPASM and the tables linked in using MPLINK.To set up the project above in MPLAB, open a new project and enter “node properties” for the target file. Under “Language Tool“, select MPLINK and then click, “OK”. This takes you back to the screen above where you can add the nodes listed.Appendix A662 Example Code ;;; DTMF Decoder Software; Ver: LINKED HALF SIZED TABLES V1.1; David Hedley CAE; Microchip Europe; November 1999;; For Example Only;;;;INCLUDE "P16C662.INC"ERRORLEVEL -302UDATASAMPLEBIT EQU 0X40W_TEMP EQU 0X41STATUS_TEMP EQU 0X42TABLECOUNTER EQU 0X43TEMP EQU 0X44LOWTEMPLATE EQU 0X45HIGHTEMPLATE EQU 0X46BYTE EQU 0X47BITVALUE EQU 0X48RESULT EQU H'28'FLAGS EQU H'29'A697SIN EQU H'30'A697COS EQU H'31'A770SIN EQU H'32'A770COS EQU H'33'A852SIN EQU H'34'A852COS EQU H'35'A941SIN EQU H'36'A941COS EQU H'37'A1209SIN EQU 0X38A1209COS EQU 0X39A1336SIN EQU 0X3AA1336COS EQU 0X3BA1447SIN EQU 0X3CA1447COS EQU 0X3DA1633SIN EQU 0X3EA1633COS EQU 0X3FEXTERN HIGHBANDEXTERN LOWBAND#DEFINE DONE FLAGS,0#DEFINE GOT_ONE FLAGS,1#DEFINE HIGHONES FLAGS,2#DEFINE ONELIT FLAGS,3begin code 0x00goto STARTinterrupt_vector code 0x04goto ISRinterrupt code 0x05ISR;****************************;; CONTEXT SAVE;;****************************C_SAVE MOVWF W_TEMPSWAPF STATUS,WBCF STATUS,5 ;ENSURE BANK0 SAVEMOVWF STATUS_TEMPBTFSS INTCON,T0IFGOTO C_RESTORE;LOOK AT THE COMPARATORBCF INTCON,T0IF;*****************************;;;;*****************************BTFSC CMCON,C1OUT ;get compar valueGOTO PIN_HIGHGOTO PIN_LOWPIN_HIGHBSF STATUS,CRLF BYTE,F ;BIT IN LSB OF BYTEGOTO CONTINUEPIN_LOWBCF STATUS,CRLF BYTE,F ;BIT IN LSB OF BYTE CONTINUE;;*****************************;MOVLW LOW LOWBANDADDWF TABLECOUNTER,fMOVLW HIGH LOWBANDBTFSC STATUS,CADDLW 1MOVWF PCLATHMOVF TABLECOUNTER,WCALL LOWBANDMOVWF LOWTEMPLATE;*****************************;; DO SOME CHECKS;;*****************************MOVLW LOW HIGHBANDADDWF TABLECOUNTER,fMOVLW HIGH HIGHBANDBTFSC STATUS,CADDLW 1MOVWF PCLATHMOVF TABLECOUNTER,WCALL HIGHBANDCLRF PCLATHMOVWF HIGHTEMPLATEMOVF TABLECOUNTER,WXORLW D'64'BTFSC STATUS,ZGOTO FINITBSF GOT_ONEINCF TABLECOUNTER,FGOTO C_RESTOREFINITBSF DONEC_RESTOREBCF STATUS,5 ;ENSURE BANK0 RESTORESWAPF STATUS_TEMP,WMOVWF STATUSSWAPF W_TEMP,FSWAPF W_TEMP,WRETFIEcode 0x100INITMOVLW 0X1FMOVWF FSRRAMCLEANINCF FSR,fCLRF INDFMOVF FSR,WXORLW 0X4FBTFSS STATUS,ZGOTO RAMCLEANMOVLW 1MOVWF TABLECOUNTERMOVLW B'11100000'MOVWF INTCONMOVLW 0X02MOVWF CMCON ;comparator setupBANKSEL TRISAMOVLW B'10100010'MOVWF VRCONCLRF TRISBCLRF TRISABSF TRISA,0MOVLW B'11000001' ;256uS SAMPLE TIMEMOVWF OPTION_REG ;START TMR0BANKSEL PORTBRETURNSTARTCLRF PORTBCLRF TMR0STARTPLUSCALL INITNOPMAINMOVF TMR0,WBTFSC DONEGOTO DUNROAMINBTFSC GOT_ONECALL GOERTZELGOTO MAINDUNROAMINBANKSEL OPTION_REGBSF OPTION_REG,5 ;STOP THE CLOCKBANKSEL PORTBCALL ACCUM_LOWCALL WINNERBCF ONELITGOTO STARTPLUSGOERTZELBCF GOT_ONE;*******************************;; LOWTEMPLATE AND HIGHTEMPLATE; CONTAIN THE EXPECTED VALUES; FOR THE INCOMING BIT;;*******************************MOVLW H'30'MOVWF FSR ;ACCUMULATOR ADDRMOVLW 8MOVWF SAMPLEBITCOMPAREBITBTFSS HIGHONESGOTO LABELLOWGOTO LABELHIGH ;BRANCH HERE LABELHIGHRLF HIGHTEMPLATE,FGOTO CARRYONLABELLOWRLF LOWTEMPLATE,FCARRYONBTFSC STATUS,CGOTO WAS1GOTO WAS0WAS1MOVLW 1GOTO CONTWAS0MOVLW 0GOTO CONTCONTMOVWF BITVALUEBTFSC BYTE,0GOTO SAMPLEWAS1GOTO SAMPLEWAS0SAMPLEWAS1MOVLW 1GOTO CONTISAMPLEWAS0MOVLW 0GOTO CONTICONTIXORWF BITVALUE,WBTFSC STATUS,ZGOTO MATCHGOTO NOMATCHNOMATCHDECF INDF,FGOTO CONTINMATCHINCF INDF,FGOTO CONTINCONTIN; move the indf to the beginning of accumulator ram; test the bit for a match and increment thru fsr;; do it 8 times for low accum band;;; reset the indf and do the compare again for the upper band. ;INCF FSR,FDECFSZ SAMPLEBIT,FGOTO COMPAREBIT ;FINISH THE BYTEBTFSC HIGHONESGOTO BOTHDONEBSF HIGHONESMOVLW 8MOVWF SAMPLEBITGOTO COMPAREBIT ;NOW GO DO THE HIGH BAND BYTE BOTHDONEBCF HIGHONES ;DO LOW NEXT TIMECLRF BYTE ;CLEAR IT OUT FOR NEXT PASSRETURNWINNERMOVLW H'30'MOVWF FSRMORERESULTSMOVF INDF,WCLRF INDFINCF FSR,FADDWF INDF,FMOVF FSR,WINCF FSR,FXORLW H'3F'BTFSS STATUS,ZGOTO MORERESULTSMOVLW 0X31 ;FIRST RESULT IS IN 0X31 THEN 0X33...MOVWF FSROUTPUTRESULTMOVF INDF,WSUBLW 0X20 ;THRESHOLDRRF RESULT,fINCF FSR,FINCF FSR,FMOVF FSR,WXORLW 0X41BTFSS STATUS,ZGOTO OUTPUTRESULTCOMF RESULT,fMOVF RESULT,WMOVWF PORTBRETURNACCUM_LOWMOVLW H'30'MOVWF FSRTESTFORNEGBTFSC INDF,7GOTO NEGGOTO NEXTNEG; SUBTRACT aX FROM 0MOVF INDF,WSUBLW 0MOVWF INDFNEXTINCF FSR,FMOVF FSR,WXORLW H'40' BTFSS STATUS,Z GOTO TESTFORNEG RETURNEND662 Map InformationMPLINK v2.00.00, LinkerLinker Map File - Created Mon Feb 14 11:35:28 2000Section InfoSection Type Address Location Size(Bytes) --------- --------- --------- --------- --------- begin code 0x000000 program 0x000002 interrupt_vector code 0x000004 program 0x000002 interrupt code 0x000005 program 0x000060 .code code 0x000100 program 0x0000fe .cinit romdata 0x00017f program 0x000004 sect2 code 0x000200 program 0x000084 sect3 code 0x000400 program 0x000084 .udata udata 0x000020 data 0x000000 Program Memory UsageStart End--------- ---------0x000000 0x0000000x000004 0x0000340x000100 0x0001800x000200 0x0002410x000400 0x000441311 out of 4101 program addresses used, program memory utilization is 7%877 Example Code ;;; DTMF DECODER SOFTWARE; VER: 2.01; FOR 877; WORKY......; DAVID HEDLEY CAE; MICROCHIP EUROPE; NOVEMBER 1999;; FOR EXAMPLE ONLY;;;;INCLUDE "P16F877.INC"ERRORLEVEL -302UDATA;****************************;;User Variables in GPR;;****************************SAMPLEBIT EQU 0X40W_TEMP EQU 0X41STATUS_TEMP EQU 0X42TABLECOUNTER EQU 0X43TEMP EQU 0X44LOWTEMPLATE EQU 0X45HIGHTEMPLATE EQU 0X46BYTE EQU 0X47BITVALUE EQU 0X48COUNT EQU H'27'RESULT EQU H'28'FLAGS EQU H'29'A697SIN EQU H'30'A697COS EQU H'31'A770SIN EQU H'32'A770COS EQU H'33'A852SIN EQU H'34'A852COS EQU H'35'A941SIN EQU H'36'A941COS EQU H'37'A1209SIN EQU 0X38A1209COS EQU 0X39A1336SIN EQU 0X3AA1336COS EQU 0X3BA1447SIN EQU 0X3CA1447COS EQU 0X3DA1633SIN EQU 0X3EA1633COS EQU 0X3FEXTERN HIGHBANDEXTERN LOWBAND#DEFINE DONE FLAGS,0#DEFINE GOT_ONE FLAGS,1#DEFINE HIGHONES FLAGS,2#DEFINE ONELIT FLAGS,3BEGIN CODE 0X00nopGOTO STARTINTERRUPT_VECTOR CODE 0X04GOTO ISRINTERRUPT CODE 0X05ISR;****************************;; CONTEXT SAVE;;****************************C_SAVE MOVWF W_TEMPSWAPF STATUS,WBCF STATUS,5 ;ENSURE BANK0 SAVEMOVWF STATUS_TEMPBTFSS INTCON,T0IF ;Exit ISR if it wasn't a TMR0 Overflow GOTO C_RESTOREBCF INTCON,T0IFBCF INTCON,INTF;*****************************;TEST LEVEL of ADCMOVF ADRESL,WXORLW 0BTFSS STATUS,ZGOTO PIN_HIGHGOTO PIN_LOWPIN_HIGHBSF STATUS,CRLF BYTE,F ;BIT IN LSB OF BYTEGOTO CONTINUEPIN_LOWBCF STATUS,CRLF BYTE,F ;BIT IN LSB OF BYTECONTINUE;*****************************;;Get the expected Lowband Bits;;*****************************MOVLW LOW LOWBANDADDWF TABLECOUNTER,fMOVLW HIGH LOWBANDBTFSC STATUS,CADDLW 1MOVWF PCLATHMOVF TABLECOUNTER,WCALL LOWBANDMOVWF LOWTEMPLATE;*****************************;;Get the expected Highband Bits;;*****************************MOVLW LOW HIGHBANDADDWF TABLECOUNTER,fMOVLW HIGH HIGHBANDBTFSC STATUS,CADDLW 1MOVWF PCLATHMOVF TABLECOUNTER,WCALL HIGHBANDCLRF PCLATHMOVWF HIGHTEMPLATEMOVF TABLECOUNTER,WXORLW D'64' ;This is the number of samples to test against (max 7f)BTFSC STATUS,ZGOTO FINITBSF GOT_ONEINCF TABLECOUNTER,FGOTO C_RESTOREFINITBSF DONE;********************************;; Context Restore;;********************************C_RESTOREBCF STATUS,5 ;ENSURE BANK0 RESTORESWAPF STATUS_TEMP,WMOVWF STATUSSWAPF W_TEMP,FSWAPF W_TEMP,WRETFIECODE 0X100INITMOVLW 0X20MOVWF FSRRAMCLEANINCF FSR,fCLRF INDFMOVF FSR,WXORLW 0X6FBTFSS STATUS,ZGOTO RAMCLEANMOVLW 1MOVWF TABLECOUNTERMOVLW B'11100000'MOVWF INTCONBANKSEL TRISAMOVLW B'10000100'MOVWF ADCON1 ;AN RA0 IO FOR 77CLRF TRISBCLRF TRISABSF TRISA,0MOVLW B'11000001' ;256US SAMPLE TIMEMOVWF OPTION_REG ;START TMR0BANKSEL PORTBRETURNSTARTCLRF PORTBSTARTPLUSCALL INITNOPMAINGETADCMOVLW B'11000001'MOVWF ADCON0BSF ADCON0,GOWAITHEREBTFSC ADCON0,GOGOTO WAITHEREBTFSC DONEGOTO DUNROAMINBTFSC GOT_ONECALL GOERTZELGOTO MAINDUNROAMINBANKSEL OPTION_REGBSF OPTION_REG,5 ;STOP THE CLOCKBANKSEL PORTBCALL ACCUM_LOWCALL WINNERGOTO STARTPLUSGOERTZELBCF GOT_ONE;*******************************;; LOWTEMPLATE AND HIGHTEMPLATE; CONTAIN THE EXPECTED VALUES; FOR THE INCOMING BIT;;*******************************MOVLW H'30'MOVWF FSR ;ACCUMULATOR ADDRMOVLW 8MOVWF SAMPLEBITCOMPAREBITBTFSS HIGHONESGOTO LABELLOWGOTO LABELHIGH ;BRANCH HERE LABELHIGHRLF HIGHTEMPLATE,FGOTO CARRYONLABELLOWRLF LOWTEMPLATE,FCARRYONBTFSC STATUS,CGOTO WAS1GOTO WAS0WAS1MOVLW 1GOTO CONTWAS0MOVLW 0GOTO CONTCONTMOVWF BITVALUEBTFSC BYTE,0GOTO SAMPLEWAS1GOTO SAMPLEWAS0SAMPLEWAS1MOVLW 1GOTO CONTISAMPLEWAS0MOVLW 0GOTO CONTICONTIXORWF BITVALUE,WBTFSC STATUS,ZGOTO MATCHGOTO NOMATCHNOMATCHDECF INDF,FGOTO CONTINMATCHINCF INDF,FGOTO CONTINCONTIN; MOVE THE INDF TO THE BEGINNING OF ACCUMULATOR RAM; TEST THE BIT FOR A MATCH AND INCREMENT THRU FSR;; DO IT 8 TIMES FOR LOW ACCUM BAND;;; RESET THE INDF AND DO THE COMPARE AGAIN FOR THE UPPER BAND. ;INCF FSR,FDECFSZ SAMPLEBIT,FGOTO COMPAREBIT ;FINISH THE BYTEBTFSC HIGHONESGOTO BOTHDONEBSF HIGHONESMOVLW 8MOVWF SAMPLEBITGOTO COMPAREBIT ;NOW GO DO THE HIGH BAND BYTE BOTHDONEBCF HIGHONES ;DO LOW NEXT TIMECLRF BYTE ;CLEAR IT OUT FOR NEXT PASSRETURNWINNERMOVLW H'30'MOVWF FSRMORERESULTSMOVF INDF,WCLRF INDFINCF FSR,FADDWF INDF,FMOVF FSR,WINCF FSR,FXORLW H'3F'BTFSS STATUS,ZGOTO MORERESULTSMOVLW 0X31MOVWF FSROUTPUTRESULT;;OUTPUT HERE;MOVF INDF,WSUBLW 0X20 ;THRESHOLDRRF RESULT,fINCF FSR,FINCF FSR,FMOVF FSR,WXORLW 0X41BTFSS STATUS,ZGOTO OUTPUTRESULTCOMF RESULT,fMOVF RESULT,WXORLW D'17'BTFSC STATUS,ZGOTO GOTCHA1MOVF RESULT,WXORLW D'33'BTFSC STATUS,ZGOTO GOTCHA2MOVF RESULT,WXORLW D'65'BTFSC STATUS,ZGOTO GOTCHA3MOVF RESULT,WXORLW D'18'BTFSC STATUS,ZGOTO GOTCHA4MOVF RESULT,WXORLW D'34'BTFSC STATUS,ZGOTO GOTCHA5MOVF RESULT,WXORLW D'66'BTFSC STATUS,ZGOTO GOTCHA6MOVF RESULT,WXORLW D'20'BTFSC STATUS,ZGOTO GOTCHA7MOVF RESULT,WXORLW D'36'BTFSC STATUS,ZGOTO GOTCHA8MOVF RESULT,WXORLW D'68'BTFSC STATUS,ZGOTO GOTCHA9MOVF RESULT,WXORLW D'40'BTFSC STATUS,ZGOTO GOTCHA0MOVF RESULT,WXORLW D'72'BTFSC STATUS,ZGOTO GOTCHAPOUNDMOVF RESULT,WXORLW D'24'BTFSC STATUS,ZGOTO GOTCHASTARCLRF RESULTGOTO SENDIT GOTCHAPOUNDMOVLW H'23'MOVWF RESULTBSF ONELITGOTO SENDIT GOTCHASTARMOVLW H'2A'MOVWF RESULTBSF ONELITGOTO SENDIT GOTCHA1MOVLW H'31'MOVWF RESULTBSF ONELITGOTO SENDIT GOTCHA2MOVLW H'32'MOVWF RESULTBSF ONELITGOTO SENDIT GOTCHA3MOVLW H'33'MOVWF RESULTBSF ONELITGOTO SENDIT GOTCHA4MOVLW H'34'MOVWF RESULTBSF ONELITGOTO SENDIT GOTCHA5MOVLW H'35'MOVWF RESULTBSF ONELITGOTO SENDIT GOTCHA6MOVLW H'36'MOVWF RESULTBSF ONELITGOTO SENDIT GOTCHA7MOVLW H'37'MOVWF RESULTBSF ONELITGOTO SENDITGOTCHA8MOVLW H'38'MOVWF RESULTBSF ONELITGOTO SENDITGOTCHA9MOVLW H'39'MOVWF RESULTBSF ONELITGOTO SENDITGOTCHA0MOVLW H'30'MOVWF RESULTBSF ONELITGOTO SENDIT;; SEND TO USART;SENDITBANKSEL SPBRGMOVLW D'25'MOVWF SPBRG ;9600 @16MEGMOVLW B'00100000'MOVWF TXSTABANKSEL TXREGBSF RCSTA,SPENMOVF RESULT,WMOVWF TXREGBSF PORTB,0BTFSC ONELITCALL WAITBCF PORTB,0RETURNWAITMOVLW H'5'MOVWF COUNTINTERLOOPCLRF TMR0BCF INTCON,GIEBCF INTCON,T0IFBANKSEL OPTION_REGMOVLW B'00000111'MOVWF OPTION_REGBANKSEL TMR0LOOPAROUNDBTFSS INTCON,T0IFGOTO LOOPAROUNDDECFSZ COUNT,fGOTO INTERLOOPRETURNACCUM_LOWMOVLW H'30'MOVWF FSRTESTFORNEGBTFSC INDF,7GOTO NEGGOTO NEXTNEG; SUBTRACT AX FROM 0 MOVF INDF,WSUBLW 0MOVWF INDFNEXTINCF FSR,FMOVF FSR,WXORLW H'40'BTFSS STATUS,ZGOTO TESTFORNEGRETURNENDAppendix D877 Map InformationMPLINK v2.00.00, LinkerLinker Map File - Created Mon Feb 14 11:26:59 2000Section InfoSection Type Address Location Size(Bytes) --------- --------- --------- --------- --------- BEGIN code 0x000000 program 0x000004 INTERRUPT_VECTOR code 0x000004 program 0x000002 INTERRUPT code 0x000005 program 0x000066 .cinit romdata 0x000038 program 0x000004 .code code 0x000100 program 0x000204 sect3 code 0x000400 program 0x000104 sect2 code 0x000600 program 0x000104 .udata udata 0x000020 data 0x000000 Program Memory UsageStart End--------- ---------0x000000 0x0000010x000004 0x0000390x000100 0x0002010x000400 0x0004810x000600 0x000681574 out of 8197 program addresses used, program memory utilization is 7%。

Machine Vision and Applications (2000) 12: 16–22机器视觉与应用（2000）12:16-22Machine Vision and Applications©Springer-Verlag 2000机器视觉与应用©施普林格出版社2000Andrea Fusiello1, Emanuele Trucco2, Alessandro Verri31 Dipartimento Scientifico e Tecnologico, Universita d i Verona, Ca’ Vignal 2, Strada Le Grazie, 37134 Verona, Italy; e-mail: fusiello@sci.univr.it `2 Heriot-Watt University, Department of Computing and Electrical Engineering, Edinburgh, UK3 INFM, Dipartimento di Informatica e Scienze dell’Informazione, Univ ersita di Genova, Genova, ItalyReceived: 25 February 1999 / Accepted: 2 March 2000收稿日期：1999年2月25日/接受日期：2000年3月2日Abstract. We present a linear rectification algorithm for general, unconstrained stereo rigs. The algorithm takes the two perspective projection matrices of the original cameras,and computes a pair of rectifying projection matrices. It is compact (22-line MATLAB code) and easily reproducible.We report tests proving the correct behavior of our method,as well as the negligible decrease of the accuracy of 3D reconstruction performed from the rectified images directly.摘要：我们在本篇文章中阐述一个用于通用的不加约束的立体视觉设备的线性修正算法。

/algorithm/20080511/8343.html动态规划算法作者：Fox 来源：C++博客酷勤网收集2008-05-11以前在学习非数值算法的时候，曾经了解过动态规划算法（Dynamic programming），以下是对Wikipedia上动态规划的翻译，图也是Wikipedia上的，仓促行文，不到之处，请方家指正。

这篇文章的术语实在是太多了，所以我在文中加入了少量注释，一律以粗斜体注明。

本文的不足之处将随时修正，MIT的《Introduction to Algorithms》第15章是专门讲动态规划的。

_____________________________________________________________动态规划在数学与计算机科学领域，动态规划用于解决那些可分解为重复子问题（overlapping subproblems，想想递归求阶乘吧）并具有最优子结构（optimal substructure，想想最短路径算法）（如下所述）的问题，动态规划比通常算法花费更少时间。

上世纪40年代，Richard Bellman最早使用动态规划这一概念表述通过遍历寻找最优决策解问题的求解过程。

1953年，Richard Bellman将动态规划赋予现代意义，该领域被IEEE纳入系统分析和工程中。

为纪念Bellman的贡献，动态规划的核心方程被命名为贝尔曼方程，该方程以递归形式重申了一个优化问题。

在“动态规划”（dynamic programming）一词中，programming与“计算机编程”（computer programming）中的programming并无关联，而是来自“数学规划”（mathematical programming），也称优化。

因此，规划是指对生成活动的优化策略。

举个例子，编制一场展览的日程可称为规划。

在此意义上，规划意味着找到一个可行的活动计划。

概述图1使用最优子结构寻找最短路径：直线表示边，波状线表示两顶点间的最短路径（路径中其他节点未显示）；粗线表示从起点到终点的最短路径。