第5篇 傅里叶递推算法
第5篇 傅里叶递推算法
一个以T 为周期的函数()t f T ，若在[]0，T -上满足狄氏条件（电网中的电压、电流满足），那么，在[]0，T -上就可以展成傅氏级数。
在计算电网中的电压、电流的基波时，存在两种算法：一种随截取不同时刻的窗（积分区间），得到不同的初相角；另一种维持初相角不变。
例如，[]11---k k t T t ，的基波值
()tdt t f T a k k t T t T k ωcos 2111?----=
，()tdt t f T b k k t T
t T k ωsin 21
11
?----=。 计算[]k k t T t ，-的基波值 第一种算法
()tdt T t f T a S t T t T k k k ωcos 211+=
?---，()tdt T t f T
b S t T t T k k k ωsin 21
1+=?---。 ()()dt t t f T a S t T t T k k k ?ω-=?-cos 2，()()dt t t f T b S t T t T k k
k ?ω-=?-sin 2。 ()1
第二种算法
()()dt T t T t f T a S S t T t T k k k ++=
?---ωcos 211，()()dt T t T t f T b S S t T t T k k k ++=?---ωsin 21
1
。 ()tdt t f T a k k t T t T k ωcos 2?-=，()tdt t f T
b k
k t T t T k ωsin 2?-=。 ()2
k k k b j a c 2
1
21+=
比较()1式与()2式，初相角差()1--==k k S S t t T ωω?。这是由于被分解函数()t f T 与相关函数t ωcos ，t ωsin 的时间差引起的。被分解函数()t f T 后移S T ，而相关函数t ωcos ，
t ωsin 未移。若相关函数同步后移S T ，就消除了初相角差S ?。
电网的应用中并不关心相量的绝对初相角，只关心它们之间的相对相角（相位差）。因
此，同时刻的相量运算，只要截取相同的窗，采用相同的算法，得到的相位差是正确的。但是，不同时刻的相量运算，也必须坚持正确的相角关系。第一种算法的窗只能相差T n ?，而第二种算法无此要求。例如计算突变量，第一种算法故障前窗超前故障后窗T n ?且随故障后窗同步推移。第二种算法固定故障前窗且靠近故障时刻，故障后窗随时间推移。直观上
()2式比()1式简单、规整，例如采用第二种算法计算
()()[]tdt T t f t f T a a k k t t T T k k ωcos 211?---=
--，()()[]tdt T t f t f T b b k
k t t T T k k ωsin 21
1?---=-- ()3
若采用第一种算法计算就相对复杂。 将()3离散得递推公式
()()[]N k t f t f N a a N k T k T k k π2cos 2
1---+
=，()()[]N
k t f t f N b b N k T k T k k π2sin 21---+= ()4
应用()4式计算故障分量。这里引入一个新概念：当前时刻t ，变化量（t ）=故障后量（t ）
-固定故障前量（譬如，记忆启动前一周波）；突变量（t ）=故障后量（t ）-故障后量（t-T ）。可见，故障后一个周波的变化量=突变量。将()4式两边同减0a ，0b 。得故障分量递推公式
()()()[]N k t f t f N a a N k T k T k k π2cos 2
010-----+?=? ()()()[]N
k t f t f N b b N k T k T k k π
2sin
2010
-----+?=?
傅立叶采取全波递推算法，无论稳态量，还是变化量。我们认为，启动前故障已经发生，约定：
启动点
故障前一个周波的记忆数据
t
故障后的第1点
故障前的最后1点
变化量全波傅氏算法
启动前预计算2个点，启动后进入故障处理程序，接着计算第3点……一直递推下去。 具体算法：令，故障前一个周波的全波傅氏计算值为零作为初值。预递推2个点，第3点在启动后递推。
()()()[]N k t f t f N a a N k T k T k k π2cos 2
010-----+
?=? ()()()[]N
k t f t f N b b N k T k T k k π
2sin
2010-----+?=? 5.2-=k 发生故障，初值=0。2-=k 为故障后的第1点。
稳态量全波傅氏算法
启动前预计算2个点，启动后进入故障处理程序，接着计算第3点……一直递推下去。 具体算法：令，故障前一个周波的全波傅氏计算值作为初值。预递推2个点，第3点在启动后递推。
()()[]N k t f t f N a a N k T k T k k π2cos 2
1---+
= ()()[]N k t f t f N b b N k T k T k k π
2sin
21---+=
可见，变化量与稳态量计算公式完全一样，仅仅是初值不同而已！故障前一个周波的采样值是必须记忆的，假设故障前一个周波的变化量为零，而稳态量是实际值。
递推的傅里叶算法计算变化量的好处在于，在相量形成的过程中，随时间推移逐渐逼近满窗逐渐准确。动作特性的裕度也随之逐渐减小直至为零。
123024232225
采样点
递推公式：
()()()()25
252
1
241
t f t f t f t f i k
i k
-=-∑∑==，()()()()24
2525
25
2
1
252
-=-=-+=∑∑t
f t f t f t f i k i k
。
递推最小二乘法算法
题目： （递推最小二乘法） 考虑如下系统： )()4(5.0)3()2(7.0)1(5.1)(k k u k u k y k y k y ξ+-+-=-+-- 式中，)(k ξ为方差为0.1的白噪声。 取初值I P 610)0(=、00=∧ ）（θ。选择方差为1的白噪声作为输入信号)(k u ，采用PLS 法进行参数估计。 Matlab 代码如下： clear all close all L=400; %仿真长度 uk=zeros(4,1); %输入初值：uk(i)表示u(k-i) yk=zeros(2,1); %输出初值 u=randn(L,1); %输入采用白噪声序列 xi=sqrt(0.1)*randn(L,1); %方差为0.1的白噪声序列 theta=[-1.5;0.7;1.0;0.5]; %对象参数真值 thetae_1=zeros(4,1); %（）θ初值 P=10^6*eye(4); %题目要求的初值 for k=1:L phi=[-yk;uk(3:4)]; %400×4矩阵phi 第k 行对应的y(k-1),y(k-2),u(k-3), u(k-4) y(k)=phi'*theta+xi(k); %采集输出数据 %递推最小二乘法的递推公式 K=P*phi/(1+phi'*P*phi); thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1); P=(eye(4)-K*phi')*P; %更新数据 thetae_1=thetae(:,k); for i=4:-1:2 uk(i)=uk(i-1); end uk(1)=u(k); for i=2:-1:2 yk(i)=yk(i-1);
几种最小二乘法递推算法的小结
一、 递推最小二乘法 递推最小二乘法的一般步骤： 1. 根据输入输出序列列出最小二乘法估计的观测矩阵?： ] )(u ... )1( )( ... )1([)(T b q n k k u n k y k y k ------=? 没有给出输出序列的还要先算出输出序列。 本例中， 2)]-u(k 1),-u(k 2),-1),-y(k -[-y(k )(T =k ?。 2. 给辨识参数θ和协方差阵P 赋初值。一般取0θ=0或者极小的数，取σσ,20I P =特别大，本例中取σ=100。 3. 按照下式计算增益矩阵G ： ) ()1()(1)()1()(k k P k k k P k G T ???-+-= 4. 按照下式计算要辨识的参数θ： )]1(?)()()[()1(?)(?--+-=k k k y k G k k T θ?θθ 5. 按照下式计算新的协方差阵P ： )1()()()1()(---=k P k k G k P k P T ? 6. 计算辨识参数的相对变化量，看是否满足停机准则。如满足，则不再递推；如不满足， 则从第三步开始进行下一次地推，直至满足要求为止。 停机准则：ε???<--) (?)1(?)(?max k k k i i i i 本例中由于递推次数只有三十次，故不需要停机准则。 7. 分离参数：将a 1….a na b 1….b nb 从辨识参数θ中分离出来。 8. 画出被辨识参数θ的各次递推估计值图形。 为了说明噪声对递推最小二乘法结果的影响，程序5-7-2在计算模拟观测值时不加噪 声， 辨识结果为a1 =1.6417，a2 = 0.7148，b1 = 0.3900，b2 =0.3499，与真实值a1 =1.642， a2 = 0.715， b1 = 0.3900，b2 =0.35相差无几。 程序5-7-2-1在计算模拟观测值时加入了均值为0，方差为0.1的白噪声序列，由于噪 声的影响，此时的结果为变值，但变化范围较小，现任取一组结果作为辨识结果。辨识结果为a1 =1.5371， a2 = 0.6874， b1 = 0.3756，b2 =0.3378。 程序5-7-2-2在计算模拟观测值时加入了有色噪声，有色噪声为 E(k)+1.642E(k-1)+0.715E(k-2)，E(k)是均值为0，方差为0.1的白噪声序列，由于有色噪声的影响，此时的辨识结果变动范围远比白噪声时大，任取一组结果作为辨识结果。辨识结果为a1 =1.6676， a2 = 0.7479， b1 = 0.4254，b2 =0.3965。 可以看出，基本的最小二乘法不适用于有色噪声的场合。
递推算法
递推算法典型例题 一、教学目标 1、由浅入深，了解递推算法 2、掌握递推算法的经典例题 二、重点难点分析 1、重点：递推关系的建立 2、难点：如何将所求问题转化为数学模型 三、教具或课件 微机 四、主要教学过程 (一)引入新课 客观世界中的各个事物之间或者一个事物的内部各元素之间，往往存在(隐藏)着很多本质上的关联。我们设计程序前．应该要通过细心的观察、丰富的联想、不断的尝试推理．尽可能先归纳总结出其内在规律，然后再把这种规律性的东西抽象成数学模型，最后再去编程实现。递推关系和递归关系都是一种简洁高效的常见数学模型，我们今天先来深入研究一下递推算法如何实现。 (二)教学过程设计 递推法是一种重要的数学方法，在数学的各个领域中都有广泛的运用，也是计算机用于数值计算的一个重要算法。这种算法特点是：一个问题的求解需一系列的计算，在已知条件和所求问题之间总存在着某种相互联系的关系，在计算时，如果可以找到前后过程之间的数量关系（即递推式），那么，这样的问题可以采用递推法来解决。从已知条件出发，逐步推出要解决的问题，叫顺推；从问题出发逐步推到已知条件，此种方法叫逆推。无论顺推还是逆推，其关键是要找到递推式。这种处理问题的方法能使复杂运算化为若干步重复的简单运算，充分发挥出计算机擅长于重复处理的特点。 递推算法的首要问题是得到相邻的数据项间的关系（即递推关系）。递推算法避开了通项公式的麻烦，把一个复杂的问题的求解，分解成了连续的若干步简单运算。一般说来可以将递推算法看成是一种特殊的迭代算法。（在解题时往往还把递推问题表现为迭代形式，用循环处理。所谓“迭代”，就是在程序中用同一个变量来存放每一次推算出来的值，每一次循环都执行同一个语句，给同一变量赋以新的值，即用一个新值代替旧值，
AIC法定阶的依阶次递推算法程序.
AIC 法定阶的依阶次递推算法程序 依阶次递推算法所得到估计θ ?，再按下式计算残差方差的估计值： ∑==N j j e N 122 )?,(1?θσ 由上式的结果计算AIC ： )(2?lg AIC 2b a n n N ++=σ 在结果中找到AIC 最小的模型（阶次和参数）就是估计的模型。 由输出数据可知当k1=5时aic 的值最小。所以最后的辨识结果取阶次为5，参数为： –1.18394，0.813938，–0.518174，0.348744，–0.116818， 1.07998，–0.74386，0.475444，–0.253022，0.122781 判断的阶次的最小aic 值： aic= – 8981.58发生在阶次为5时。 源程序： #include
实验一用递推公式计算定积分
实验一 用递推公式计算定积分 09信息 符文飞 07 1、实验目的： 由于一个算法是否稳定，十分重要。如果算法不稳定，则数值计算的结果就会严重背离数学模型的真实结果，因此，在选择数值计算公式来进行近似计算时，我们应特别注意选用那些在数值计算过程中不会导致误差迅速增长的公式。体会稳定性在选择算法中的地位．误差扩张的算法是不稳定的，是我们所不期望的；误差衰竭的算法是稳定的．是我们努力寻求的，这是贯穿本课程的目标．通过上机计算，了解舍入误差所引起的数值不稳定性。 2、实验题目： 对n =0,1,2,…,20,计算定积分dx x x y n n ?+=10 5 3、实验原理 由于y(n)= = – 在计算时有两种迭代方法，如下： 方法一： y(n)= – 5*y(n-1),n=1,2,3, (20) 取y(0)= = ln6-ln5 ≈ 0.182322 方法二：
利用递推公式：y(n-1)=-*y(n),n=20,19, (1) 而且，由 = * ≤≤* = 可取：y(20)≈*()≈0.008730. 4、实验内容： 算法1的程序： y0=log(6.0)-log(5.0); y1=0; n=1; while n<=30 y1=1/n-5*y0; fprintf('y[%d]=%-20f',n,y1); y0=y1; n=n+1; if mod(n,1)==0; fprintf('\n') end end 算法2的程序： y0=(1/105+1/126)/2;
y1=0; n=1; while n<=30 y1=1/(5*n)-y0/5; fprintf('y[%d]=%-20f',n,y1) y0=y1; n=n+1; if mod(n,1)==0 fprintf('\n') end end 5、实验结果 对于算法1： y[1]=0.088392 y[2]=0.058039 y[3]=0.043139 y[4]=0.034306 y[5]=0.028468 y[6]=0.024325 y[7]=0.021233 y[8]=0.018837 y[9]=0.016926
递推算法
递推算法 一、学习要点：（经典的五种递推关系） 1.fibonacci（斐波那契）数列 2.hanoi（汉诺）塔 3.平面分割 4.catalan数（卡特兰数） 5.第二类string数 二、重点提示 1.产生catalan数的参考程序： Const max=21; Var c:array[2..max] of longint; n,i,k:integer; total:longint; begin write(‘input n=’);readln(n); c[2]:=1; for i:=3 to n do begin c[i]:=0; for k:=2 to i-1 do c[i]:=c[i]+c[k]*c[i-k+1]; end; writeln(‘catalan=’,c[n]); end. Catalan数，从第二项开始为：1,1,2,5,14,42,132,429,1430,4862,16796,58786,208012, 742900,2674440,9694545,35357670,129644790,477638700,1767263190 2.产生第二类的string数的参考程序： Var n,k:integer; Function s(n,k:integer):longint; Begin If (n 实用算法(基础算法-递推法-01) 有一类试题，每相邻两项数之间的变化有一定的规律性，我们可将这种规律归纳成如下简捷的递推关系式： F n=g(F n-1) 这就在数的序列中，建立起后项和前项之间的关系，然后从初始条件(或最终结果)入手，一步步地按递推关系递推，直至求出最终结果(或初始值)。很多程序就是按这样的方法逐步求解的。如果对一个试题，我们要是能找到后一项与前一项的关系并清楚其起始条件(最终结果)，问题就好解决，让计算机一步步算就是了，让高速的计算机做这种重复运算，可真正起到“物尽其用”的效果。 递推分倒推法和顺推法两种形式。一般分析思路： if求解条件F1 then begin{倒推} 由题意(或递推关系)确定最终结果Fa； 求出倒推关系式F i-1=g'(F i); i=n;{从最终结果Fn出发进行倒推} while 当前结果F i非初始值F1 do由F i-1=g(F1)倒推前项； 输出倒推结果F1和倒推过程； end {then} else begin{顺推} 由题意(或顺推关系)确定初始值F1(边界条件)； 求出顺推关系式F1=g(Fi-1); i=1;{由边界条件F1出发进行顺推} while 当前结果Fi非最终结果Fn do由Fi=g(Fi-1)顺推后项； 输出顺推结果Fn和顺推过程； end; {else} 一、倒推法 所谓倒推法，就是在不知初始值的情况下，经某种递推关系而获知问题的解或目标，再倒推过来，推知它的初始条件。因为这类问题的运算过程是一一映射的，故可分析得其递推公式。然后再从这个解或目标出发，采用倒推手段，一步步地倒推到这个问题的初始陈述。 下面举例说明。 [例1] 贮油点 一辆重型卡车欲穿过1000公里的沙漠，卡车耗油为1升/公里，卡车总载油能力为500公升。显然卡车一次是过不了沙漠的。因此司机必须设法在沿途建立几个储油点，使卡车能顺利穿越沙漠，试问司机如何建立这些储油点？每一储油点应存多 递推算法集锦 一、编写程序求50以内的勾股弦数。即满足c*c=b*b+a*a的三个自然数，要求 b>a。将所有符合要求的a,b,c组合输出至屏幕。 解答： 采用穷举算法，在主函数中实现。 #include 递推最小二乘法的一般步骤： 1. 根据输入输出序列列出最小二乘法估计的观测矩阵?： ] )(u ... )1( )( ... )1([)(T b q n k k u n k y k y k ------=? 没有给出输出序列的还要先算出输出序列。 本例中， 2)]-u(k 1),-u(k 2),-1),-y(k -[-y(k )(T =k ?。 2. 给辨识参数θ和协方差阵P 赋初值。一般取0θ=0或者极小的数，取σσ,20I P =特别 大，本例中取σ=100。 3. 按照下式计算增益矩阵G ： ) ()1()(1)()1()(k k P k k k P k G T ???-+-= 4. 按照下式计算要辨识的参数θ： )]1(?)()()[()1(?)(?--+-=k k k y k G k k T θ?θθ 5. 按照下式计算新的协方差阵P ： )1()()()1()(---=k P k k G k P k P T ? 6. 计算辨识参数的相对变化量，看是否满足停机准则。如满足，则不再递推；如不满足， 则从第三步开始进行下一次地推，直至满足要求为止。 停机准则：ε???<--) (?)1(?)(?max k k k i i i i 本例中由于递推次数只有三十次，故不需要停机准则。 7. 分离参数：将a 1….a na b 1….b nb 从辨识参数θ中分离出来。 8. 画出被辨识参数θ的各次递推估计值图形。 为了说明噪声对递推最小二乘法结果的影响，程序5-7-2在计算模拟观测值时不加噪声， 辨识结果为a1 =1.6417，a2 = 0.7148，b1 = 0.3900，b2 =0.3499，与真实值a1 =1.642， a2 = 0.715， b1 = 0.3900，b2 =0.35相差无几。 程序5-7-2-1在计算模拟观测值时加入了均值为0，方差为0.1的白噪声序列，由于噪声 的影响，此时的结果为变值，但变化范围较小，现任取一组结果作为辨识结果。辨识结果为a1 =1.5371， a2 = 0.6874， b1 = 0.3756，b2 =0.3378。 程序5-7-2-2在计算模拟观测值时加入了有色噪声，有色噪声为 E(k)+1.642E(k-1)+0.715E(k-2)，E(k)是均值为0，方差为0.1的白噪声序列，由于有色噪声的影响，此时的辨识结果变动范围远比白噪声时大，任取一组结果作为辨识结果。辨识结果为a1 =1.6676， a2 = 0.7479， b1 = 0.4254，b2 =0.3965。 可以看出，基本的最小二乘法不适用于有色噪声的场合。 求递推数列通项公式的常用方法归纳 目录 一、概述·································· 二、等差数列通项公式和前n项和公式·································· 1、等差数列通项公式的推导过程································ 2、等差数列前n项和公式的推导过程·································· 三、一般的递推数列通项公式的常用方法·································· 1、公式法·································· 2、归纳猜想法·································· 3、累加法·································· 4、累乘法·································· 5、构造新函数法(待定系数法）·································· 6、倒数变换法·································· 7、特征根法·································· 8、不动点法································· 9、换元法································· 10、取对数法·································· 11、周期法··································实用算法(基础算法-递推法)
递推算法解析集锦
几种最小二乘法递推算法的小结
根据递推公式求数列通项公式的常用方法总结归纳(新)