基于的神經(jīng)網(wǎng)絡(luò)地印刷體字母識別

1、實用標準文案基于 BP神經(jīng)網(wǎng)絡(luò)的印刷體字母識別1 背景隨著社會的發(fā)展， 英語作為國際通用語言得到了日益廣泛的應(yīng)用，因此有大量的英文文檔整理、 查詢、統(tǒng)計的工作需要完成， 而英文字母識別系統(tǒng)可以輕而易舉地完成很多以前難以想象的工作。智能控制作為一門新興的交叉學科，在許多方面都優(yōu)于傳統(tǒng)控制， 而智能控制中的人工神經(jīng)網(wǎng)絡(luò)由于模仿人類的神經(jīng)網(wǎng)絡(luò)，具有感知識別、學習、聯(lián)想、記憶、推理等智能，更是有著廣闊的發(fā)展前景。人工神經(jīng)網(wǎng)絡(luò)理論的應(yīng)用主要在人工智能, 自動控制 , 模式識別 , 機器人 , 信息處理 ,CAD/CAM等方面。如 :(1) 空間科學。航空飛行器及汽車的自動駕駛導航系統(tǒng) , 飛行路徑模擬，

2、 飛行器制導和飛行程序優(yōu)化管理等。(2) 控制和優(yōu)化。機器人運動控制 , 各種工業(yè)過程控制和制造過程控制， 如集成電路布線設(shè)計 , 生產(chǎn)流程控制等等。(3) 模式識別和圖像處理。如人臉識別 , 語言識別 , 指紋識別 , 簽名識別 , 手寫體和印刷體字符識別 , 目標檢測與識別 , 圖像復原 , 圖像壓縮等等。(4) 智能信息管理系統(tǒng)。如股價預測 , 不動產(chǎn)價格預測 , 外匯 , 黃金等大宗產(chǎn)品價格預測 , 公司財務(wù)分析 , 地震及各種自然災害預報等等。其中最核心的是反向傳播網(wǎng)絡(luò)(Back Propagation Network),簡稱 BP 網(wǎng)絡(luò)。本文介紹了運用matlab 工具箱確定隱層神

3、經(jīng)元的個數(shù)和構(gòu)造BP神經(jīng)網(wǎng)絡(luò)，并用兩組樣本對該神經(jīng)網(wǎng)絡(luò)進行訓練，然后運用訓練后的網(wǎng)絡(luò)對字母進行識別。2 BP 網(wǎng)絡(luò)介紹BP神經(jīng)網(wǎng)絡(luò)又稱誤差反向傳遞神經(jīng)網(wǎng)絡(luò)。它是一種依靠反饋值來不斷調(diào)整節(jié)點之間的連接權(quán)值而構(gòu)建的一種網(wǎng)絡(luò)模型。它的整個體系結(jié)構(gòu)分為輸入層、隱藏層和輸出層 , 其中隱藏層根據(jù)具體情況的需要, 可以是一層結(jié)構(gòu)也可為多層結(jié)構(gòu)。BP 算法的基本思想是: 學習過程由信號的正向傳播與誤差的反向傳播兩個過程組成。正向傳播時 , 輸入樣本從輸入層傳入 , 經(jīng)各隱藏層逐層處理后， 傳向輸出層。若輸出層的實際輸出與期望的輸出 ( 教師信號 ) 不符 , 則轉(zhuǎn)入誤差的反向傳播階精彩文檔實用標準文案段。誤

4、差反傳是將輸出誤差以某種形式通過隱藏層向輸入層反傳, 并將誤差分攤給各層的所有單元 , 從而獲得各層單元的誤差信號, 此誤差信號即作為修正各單元權(quán)值的依據(jù)。這種信號正向傳播與誤差反向傳的各層權(quán)值調(diào)整過程 , 是周而復始地進行的。權(quán)值不斷調(diào)整的過程 , 也就是網(wǎng)絡(luò)的學習訓練過程。此過程一直進行到網(wǎng)絡(luò)輸出的誤差減少到可接受到的程度, 或進行到預先設(shè)定的學習次數(shù)為此。3 系統(tǒng)實現(xiàn)思想字母識別系統(tǒng)一般分為預處理、 特征提取和分類器三部分。 其中，預處理包括將圖片從模擬圖像到進行二值化、 歸一化等過程； 特征提取和分類器的設(shè)計是整個系統(tǒng)的核心部分。 通過對各個部分分別進行編程處理來實現(xiàn)， 將每部分編為可

5、調(diào)用的函數(shù)，最后統(tǒng)一對函數(shù)進行調(diào)用，清晰方便。3.1 字母識別整體框圖a BP 神經(jīng)網(wǎng)絡(luò)訓練過程b BP 神經(jīng)網(wǎng)絡(luò)識別過程圖 2.1 BP 神經(jīng)網(wǎng)絡(luò)識別系統(tǒng)3.2 預處理及其特征提取方法本文使用傅里葉描述符及其反變換進行圖片的二值化、 字母輪廓提取， 之后進行歸一化，將其特征變成 1*120 的矩陣，之后選取里面的六十個點變?yōu)?1*60 的矩陣。特征提取程序：function FD=Feature_Building(RGB)%RGB=imread('d:A.bmp');精彩文檔實用標準文案%figure(1),inshow(RGB)B=outline(RGB);%figure(

6、2)%subplot(221),draw_outline(B);%title('outline of object');m,n=size(B);FD=fsd(B,30,m,4);其中 outline、 fsd 為傅里葉描述及其反變換程序。outline程序：%Function for extracting outline of object; Q.K., 2008.4.29 %Deaprtment of Automation, Tsinghua Univ. Beijing 100084, China. function outline=outline(RGB)I=rgb2gra

7、y(RGB);junk threshold = edge(I, 'sobel');fudgeFactor = .5;BWs = edge(I,'sobel', threshold * fudgeFactor);%Step 3: Dilate the imagese90 = strel('line', 3, 90);se0 = strel('line', 3, 0);BWsdil = imdilate(BWs, se90 se0);%Step 4: Fill interior gapsBWdfill = imfill(BWsdil,

8、 'holes');%Step 5: Remove connected objects on borderBWnobord = imclearborder(BWdfill, 4);%Step 6: Smoothen the objectseD = strel('diamond',1);BWfinal = imerode(BWnobord,seD);BWfinal = imerode(BWfinal,seD);bw = bwareaopen(BWfinal,30); % % fill a gap in the pen's cap B,L = bwbound

9、aries(bw,'noholes');outline = B1;fsd 程序見程序清單。3.3 BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)3.3.1輸入層神經(jīng)元個數(shù)的確定將圖像的特征向量作為神經(jīng)網(wǎng)絡(luò)的輸入， 所以神經(jīng)網(wǎng)絡(luò)的輸入層神經(jīng)元個數(shù)等于特征向量的維數(shù)，即 1×60=60個輸入神經(jīng)元。3.3.2隱含層神經(jīng)元個數(shù)的確定隱層節(jié)點數(shù)對網(wǎng)絡(luò)的學習和計算特性具有非常重要的影響， 是該網(wǎng)絡(luò)結(jié)構(gòu)成敗的關(guān)鍵。若隱層節(jié)點數(shù)過少，則網(wǎng)絡(luò)難以處理復雜的問題；但若隱層節(jié)點數(shù)過多，則將使網(wǎng)絡(luò)學習時間急劇增加， 而且還可能導致網(wǎng)絡(luò)學習過度， 抗干擾能力下降。本文根據(jù)實際的實驗，確定隱含層神經(jīng)元的個數(shù)為15個。精彩

10、文檔實用標準文案3.3.3輸出層神經(jīng)元個數(shù)的確定因為要識別 26個英文大寫字母，因此輸出選擇 26×1的矩陣，即輸出層神經(jīng)元的個數(shù)為 26個。當 26個字母輸入神經(jīng)網(wǎng)絡(luò)后，在對應(yīng)的位置上輸出 1，其他位置上輸出零。當網(wǎng)絡(luò)進入識別過程時， 哪個位置上輸出的期望值最大， 認為識別出的是這個位置上的字母。3.3.4 BP 神經(jīng)網(wǎng)絡(luò)的構(gòu)造建立一個前向 BP神經(jīng)網(wǎng)絡(luò)函數(shù) newff:net=newff(minmax(P),S1,S2, logsig , logsig , trainlm ); net.LW2,1=net.LW2,1*0.01;net.b2=net.b2*0.01;其中 minm

11、a(P)為神經(jīng)網(wǎng)絡(luò)的對它的 60個輸入元素的最大值和最小值的限制。 P為訓練樣本集合。 S1、S2分別為該神經(jīng)網(wǎng)絡(luò)的隱含層和輸出層的神經(jīng)元個數(shù)。 logsi g , logsig 為神經(jīng)網(wǎng)絡(luò)的各層的轉(zhuǎn)移函數(shù) , 均設(shè)置為對數(shù) S型激活函數(shù)。訓練函數(shù)采用 trainlm 。3.4 BP 神經(jīng)網(wǎng)絡(luò)的訓練3.4.1訓練樣本集合和目標值集合字母圖片歸一化后的圖像為60× 1的矩陣，用 60× 26的矩陣形成一個訓練樣本；目標矢量是希望每一個數(shù)字輸入神經(jīng)網(wǎng)絡(luò)后在輸出神經(jīng)元對應(yīng)的位置上為 1, 其他的位置為 0。為此取目標矢量為對角線上為 1的26×26的單位陣 , 用mat

12、lab 命令實現(xiàn)為 :T=eye(26);3.4.2網(wǎng)絡(luò)訓練隱含層神經(jīng)元的傳遞函數(shù)采用s型對數(shù)函數(shù) logsig ，輸出層神經(jīng)元傳遞函數(shù)也采用 s型對數(shù)函數(shù)，訓練函數(shù)采用 trainlm ，性能函數(shù)采用 sse，訓練步數(shù)設(shè)置為最大 5000，性能目標值為 0.05, 。BP訓練程序：net.performFcn='sse'%設(shè)置目標性能函數(shù)net.trainParam.goal=0.05;%性能目標值net.trainParam.show=20;%顯示間隔次數(shù)net.trainParam.epochs=5000;%最大訓練次數(shù)net.trainParam.mc=0.95;net

13、,tr=train(net,P,T);BP網(wǎng)絡(luò)訓練流程圖：精彩文檔實用標準文案使用第一組樣本進行訓練的結(jié)果：TRAINLM, Epoch 0/5000, SSE 169.303/0.05, Gradient 39.0748/1e-010 TRAINLM, Epoch 20/5000, SSE 9.07917/0.05, Gradient 0.647529/1e-010 TRAINLM, Epoch 40/5000, SSE 5.45171/0.05, Gradient 0.465044/1e-010 TRAINLM, Epoch 60/5000, SSE 3.85999/0.05, Gradi

14、ent 1.13736/1e-010 TRAINLM, Epoch 80/5000, SSE 3.37108/0.05, Gradient 0.970379/1e-010 TRAINLM, Epoch 100/5000, SSE 1.43394/0.05, Gradient 0.27961/1e-010 TRAINLM, Epoch 120/5000, SSE 1.13878/0.05, Gradient 0.661835/1e-010 TRAINLM, Epoch 140/5000, SSE 0.561939/0.05, Gradient 0.497918/1e-010 TRAINLM, E

15、poch 160/5000, SSE 0.537153/0.05, Gradient 0.0963243/1e-010 TRAINLM, Epoch 180/5000, SSE 0.518194/0.05, Gradient0.00990168/1e-010 TRAINLM, Epoch 200/5000, SSE 0.461637/0.05, Gradient 11.4576/1e-010 TRAINLM, Epoch 206/5000, SSE 0.0350697/0.05, Gradient 0.265104/1e-010 TRAINLM, Performance goal met.可見

16、經(jīng)過 206次訓練后，網(wǎng)絡(luò)誤差達到要求，誤差曲線如下圖：精彩文檔實用標準文案使用第二組樣本進行訓練的結(jié)果：TRAINLM, Epoch 0/5000, SSE 168.635/0.05, Gradient 33.7987/1e-010 TRAINLM, Epoch 20/5000, SSE 3.28669/0.05, Gradient 40.5407/1e-010 TRAINLM, Epoch 32/5000, SSE 0.0441687/0.05, Gradient 0.0844925/1e-010 TRAINLM, Performance goal met.可見經(jīng)過 26次訓練之后，網(wǎng)絡(luò)誤

17、差達到要求。誤差曲線如下圖所示：3.5 字母的識別精彩文檔實用標準文案以上所介紹為網(wǎng)絡(luò)的學習期， 學習過程結(jié)束后， 網(wǎng)絡(luò)進入工作期， 即可以進行字母的識別。單一字母識別程序如下：新建文件夾 1A11.bmp');% 工作期 A11 為大寫字母 A 略帶噪聲的圖片F(xiàn)DB=Feature_Building(RGB);% 提取字母輪廓特征FDB=reshape(FDB,1,120);FDB=FDB(1:2:120);%歸一化處理a,b=max(sim(net,(FDB)')% 字母識別 a 為網(wǎng)絡(luò)工作后輸出層輸出的最大值，b為所識別字母的行數(shù) 即如果識別為 A，則 b為 1；識別為

18、B，則 b為 2，以此類推。識別結(jié)果如下：a=0.8316b=1可見能夠正確識別字母 A。新建文件夾 1B11.bmp');% 工作期 B11 為大寫字母 B 略帶噪聲的圖片F(xiàn)DB=Feature_Building(RGB);% 提取字母輪廓特征FDB=reshape(FDB,1,120);FDB=FDB(1:2:120);%歸一化處理a,b=max(sim(net,(FDB)')% 字母識別 a 為網(wǎng)絡(luò)工作后輸出層輸出的最大值，b為所識別字母的行數(shù) 即如果識別為 A，則 b為 1；識別為 B，則 b為 2，以此類推。識別結(jié)果如下：a=0.9741b=2可見也能夠正確的識別字母

19、B。本文使用兩組樣本進行 BP神經(jīng)網(wǎng)絡(luò)的訓練， 一組樣本進行字母的識別。 識別程序及其結(jié)果如下：新建文件夾 1index1.mat')for i=1:26RGBi=imread('D:Program新建文件夾1',index1i);FDi=Feature_Building(RGBi);FDi=reshape(FDi,1,120);FDi=FDi(1:2:120);a,b=max(sim(net,(FDi)')end結(jié)果：a =0.8316b =1a =0.9741b =2a =0.8805b =3a =0.9315b =4精彩文檔實用標準文案a =0.6114b

20、 =5a =0.9755b =6a =0.9715b =7a =0.9780b =8a =0.9770b =9a =0.9958b =10a =0.8759b =11a =0.9610b =12a =0.9695b =13a =0.8119b =8a =0.9718b =15a =0.9752b =5a =0.9039b =17a =0.5457b =18a =0.8177b =6a =0.9728b =20a =0.2953b =18a =0.8482b =22a =0.9092b =23a =0.9743b =24a =0.9534b =25a =0.9764b =26精彩文檔實用標準文案

21、由以上結(jié)果可知：識別了 22個字母，有四個字母未被正確識別（ N P S U）。為了使識別準確率更高， 訓練更多的樣本， 盡量選擇一些略帶有噪聲的圖片， 識別時準確率更高。精彩文檔實用標準文案程序清單BP網(wǎng)絡(luò)訓練程序：clcclear新建文件夾 1index.mat')for i=1:52RGBi=imread('D:Program新建文件夾1',indexi);FDi=Feature_Building(RGBi);FDi=reshape(FDi,1,120);FDi=FDi(1:2:120);endP1=(FD1)' (FD2)' (FD3)'

22、 (FD4)' (FD5)' (FD6)' (FD7)'(FD8)'(FD9)'(FD10)'(FD11)'(FD12)'(FD13)'(FD14)'(FD15)' (FD16)' (FD17)' (FD18)' (FD19)' (FD20)' (FD21)'(FD22)' (FD23)' (FD24)' (FD25)' (FD26)'P2=(FD27)'(FD28)'(FD29)'(FD

23、30)'(FD31)'(FD32)'(FD33)' (FD34)' (FD35)' (FD36)' (FD37)' (FD38)' (FD39)'(FD40)' (FD41)' (FD42)' (FD43)' (FD44)' (FD45)' (FD46)'(FD47)' (FD48)' (FD49)' (FD50)' (FD51)' (FD52)'%P=P1;P2;T=eye(26);S1=15;S2=26;fo

24、r n=1:2 %學習期net=newff(minmax(Pn),S1 S2,'logsig' 'logsig','trainlm');net.LW2,1=net.LW2,1*0.01;net.b2=net.b2*0.01;net.performFcn='sse'net.trainParam.goal=0.05;net.trainParam.show=20;net.trainParam.epochs=5000;net.trainParam.mc=0.95;net,tr=train(net,Pn,T);end識別程序：新建文件夾 1

25、index1.mat')for i=1:26RGBi=imread('D:Program新建文件夾1',index1i);FDi=Feature_Building(RGBi);FDi=reshape(FDi,1,120);FDi=FDi(1:2:120);a,b=max(sim(net,(FDi)')精彩文檔實用標準文案end傅里葉變換程序：function FD=Feature_Building(RGB)%RGB=imread('d:A.bmp');%figure(1),inshow(RGB)B=outline(RGB);%figure(2)%

26、subplot(221),draw_outline(B);%title('outline of object');m,n=size(B);FD=fsd(B,30,m,4);%Function for extracting outline of object; Q.K., 2008.4.29 %Deaprtment of Automation, Tsinghua Univ. Beijing 100084, China.function outline=outline(RGB)I=rgb2gray(RGB);junk threshold = edge(I, 'sobel&#

27、39;);fudgeFactor = .5;BWs = edge(I,'sobel', threshold * fudgeFactor);%Step 3: Dilate the imagese90 = strel('line', 3, 90);se0 = strel('line', 3, 0);BWsdil = imdilate(BWs, se90 se0);%Step 4: Fill interior gapsBWdfill = imfill(BWsdil, 'holes');%Step 5: Remove connected

28、objects on borderBWnobord = imclearborder(BWdfill, 4);%Step 6: Smoothen the objectseD = strel('diamond',1);BWfinal = imerode(BWnobord,seD);BWfinal = imerode(BWfinal,seD);bw = bwareaopen(BWfinal,30); % % fill a gap in the pen's cap B,L = bwboundaries(bw,'noholes');outline = B1;fun

29、ction rFSDs = fsd(outline,H,b,bN)% Forward elliptical Fourier transform - see Kuhl FP and Giardina CR% "Elliptic Fourier features of a closed contour" Computer Graphics and精彩文檔實用標準文案% Image Processing 18:236-258 1982 for theory.% Returns a shape spectrum of input x,y data "outline&quo

30、t; with% iNoOfHarmonicsAnalyse elements.% The output FSDs will be normalised for location, size and orientation% if bNormaliseSizeState and bNormaliseOrientationState are TRUE% Pre-calculate some constant arrays% n * 2 * pi% n2 * 2* pi2% where n is the number of harmonics to be used in the analysis%

31、H = iNoOfHarmonicsAnalyse%b = bNormaliseSizeState %m n = size(outline), b = m; %bN = bNormaliseOrientationStaterTwoNPi = (1:1:H)* 2 * pi;rTwoNSqPiSq = (1:1:H) .* (1:1:H)* 2 * pi * pi;iNoOfPoints = size(outline,1) - 1; % hence there is 1 more data point in outline than iNoOfPointsrDeltaX = zeros(iNoO

32、fPoints+1,1); % pre-allocate some arrays rDeltaY = zeros(iNoOfPoints+1,1); rDeltaT = zeros(iNoOfPoints+1,1);for iCount = 2 : iNoOfPoints + 1rDeltaX(iCount-1) = outline(iCount,1) - outline(iCount-1,1); rDeltaY(iCount-1) = outline(iCount,2) - outline(iCount-1,2);end%Calculate 'time' difference

33、s from point to point - actually distances, but we are% carrying on the fiction of a point running around the closed figure at constant speed.%Weare analysing the projections on to the x and y axes of this point's path around the figurefor iCount = 1 : iNoOfPointsrDeltaT(iCount) = sqrt(rDeltaX(i

34、Count)2) + (rDeltaY(iCount)2); endcheck = (rDeltaT = 0); % remove zeros from rDeltaT, rDeltaX.rDeltaT = rDeltaT(check);rDeltaX = rDeltaX(check);rDeltaY = rDeltaY(check);iNoOfPoints = size(rDeltaT,1) - 1; % we have removed duplicate points精彩文檔實用標準文案% now sum the incremental times to get the time at a

35、ny point rTime(1) = 0;for iCount = 2 : iNoOfPoints + 1rTime(iCount) = rTime(iCount - 1) + rDeltaT(iCount-1); endrPeriod = rTime(iNoOfPoints+1); % rPeriod defined for readability% calculate the A-sub-0 coefficient rSum1 = 0;for iP = 2 : iNoOfPoints + 1rSum2 = 0;rSum3 = 0;rInnerDiff = 0;% calculate th

36、e partial sums - these are 0 for iCount = 1 if iP > 1for iJ = 2 : iP-1rSum2 = rSum2 + rDeltaX(iJ-1);rSum3 = rSum3 + rDeltaT(iJ-1);endrInnerDiff= rSum2 -(rDeltaX(iP-1)/rDeltaT(iP-1)*rSum3);endrIncr1=(rDeltaX(iP-1)/(2*rDeltaT(iP-1)*(rTime(iP)2-rTime(iP-1)2)+rInnerDiff*(rTime(iP)-rTime(iP-1);rSum1 =

37、 rSum1 + rIncr1;endrFSDs(1,1) = (1 / rPeriod) * rSum1) + outline(1,1); % store A-sub-0 in output FSDs array - this array will be 4 x iNoOfHarmonicsAnalyse% calculate the a-sub-n coefficients for iHNo = 2 : HrSum1 = 0;for iP = 1 : iNoOfPointsrIncr1=(rDeltaX(iP)/rDeltaT(iP)*(cos(rTwoNPi(iHNo-1)*rTime(

38、iP+1)/rPeriod)-cos(rTwoNPi(iHNo-1)*rTime(iP)/rPeriod);rSum1 = rSum1 + rIncr1;endrFSDs(1,iHNo) = (rPeriod / rTwoNSqPiSq(iHNo-1) * rSum1; end % "foriHNo = 1 :."精彩文檔實用標準文案rFSDs(2,1) = 0; % there is no 0th order sine coefficient% calculate the b-sub-n coefficients for iHNo = 2 : HrSum1 = 0;for

39、 iP = 1 : iNoOfPointsrIncr1=(rDeltaX(iP)/rDeltaT(iP)*(sin(rTwoNPi(iHNo-1)*rTime(iP+1)/rPeriod)-sin(rTwoNPi(iHNo-1)*rTime(iP)/rPeriod);rSum1 = rSum1 + rIncr1;endrFSDs(2,iHNo) = (rPeriod / rTwoNSqPiSq(iHNo-1) * rSum1; end % "foriHNo = 1 :."% calculate the C-sub-0 coefficient rSum1 = 0;for iP

40、 = 2 : iNoOfPoints + 1rSum2 = 0;rSum3 = 0;rInnerDiff = 0;% calculate the partial sums - these are 0 for iCount = 1 if iP > 1for iJ = 2 : iP-1rSum2 = rSum2 + rDeltaY(iJ-1); rSum3 = rSum3 + rDeltaT(iJ-1);endrInnerDiff = rSum2 - (rDeltaY(iP-1) / rDeltaT(iP-1) * rSum3); endrIncr1=(rDeltaY(iP-1)/(2*rD

41、eltaT(iP-1)*(rTime(iP)2-rTime(iP-1)2)+rInnerDiff*(rTime(iP)-rTime(iP-1);rSum1 = rSum1 + rIncr1;endrFSDs(3,1) = (1 / rPeriod) * rSum1) + outline(1,2); % store C-sub-0 in output FSDs array - this array will be 4 x iNoOfHarmonicsAnalyse% calculate the C-sub-n coefficients for iHNo = 2 : HrSum1 = 0;for

42、iP = 1 : iNoOfPointsrIncr1=(rDeltaY(iP)/rDeltaT(iP)*(cos(rTwoNPi(iHNo-1)*rTime(iP+1)/rPeriod) - cos(rTwoNPi(iHNo-1)*rTime(iP)/rPeriod);rSum1 = rSum1 + rIncr1;精彩文檔實用標準文案endrFSDs(3,iHNo) = (rPeriod / rTwoNSqPiSq(iHNo-1) * rSum1; end % "foriHNo = 1 :."rFSDs(4,1) = 0; % there is no 0th order s

43、ine coefficient% calculate the D-sub-n coefficients for iHNo = 2 : HrSum1 = 0;for iP = 1 : iNoOfPointsrIncr1=(rDeltaY(iP)/rDeltaT(iP)*(sin(rTwoNPi(iHNo-1)*rTime(iP+1)/rPeriod) - sin(rTwoNPi(iHNo-1)*rTime(iP)/rPeriod);rSum1 = rSum1 + rIncr1;endrFSDs(4,iHNo) = (rPeriod / rTwoNSqPiSq(iHNo-1) * rSum1; e

44、nd % "foriHNo = 1 :.% the non-normalised coefficients are now in rFSDs% if we want the normalised ones, this is where it happens if (b = 1) | (bN = 1)% rTheta1 is the angle through which the starting position of the first % harmonic phasor must be rotated to be aligned with the major axisof% th

45、e first harmonic ellipse rFSDsTemp = rFSDs;rTheta1 = 0.5 * atan(2 * (rFSDsTemp(1,2) * rFSDsTemp(2,2) + rFSDsTemp(3,2) * rFSDsTemp(4,2) / .(rFSDsTemp(1,2)2 + rFSDsTemp(3,2)2 - rFSDsTemp(2,2)2 - rFSDsTemp(4,2)2);% calculate the partially normalised coefficients - normalised for% starting pointfor iHNo = 1 : HrStarFSDs(1,iHN