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

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

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

3、P神經(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算法的基本思想是:學(xué)習(xí)過程由信號的正向傳播與誤差的反向傳播兩個過程組成。正向傳播時,輸入樣本從輸入層傳入,經(jīng)各隱藏層逐層處理后，傳向輸出層。若輸出層的實際輸出與期望的輸出(教師信號)不符,則轉(zhuǎn)入誤差的反向傳播階段。誤差反傳是將輸出誤差以某種形式通過隱藏層向輸入層反傳,并將誤差分?jǐn)偨o各層的所有單元,從而獲得各層單元的誤差信號,此誤差信號即作為修正各單元權(quán)值的依據(jù)。這種信號正向傳播與誤差反

4、向傳的各層權(quán)值調(diào)整過程,是周而復(fù)始地進(jìn)行的。權(quán)值不斷調(diào)整的過程,也就是網(wǎng)絡(luò)的學(xué)習(xí)訓(xùn)練過程。此過程一直進(jìn)行到網(wǎng)絡(luò)輸出的誤差減少到可接受到的程度,或進(jìn)行到預(yù)先設(shè)定的學(xué)習(xí)次數(shù)為此。3系統(tǒng)實現(xiàn)思想字母識別系統(tǒng)一般分為預(yù)處理、特征提取和分類器三部分。其中，預(yù)處理包括將圖片從模擬圖像到進(jìn)行二值化、歸一化等過程；特征提取和分類器的設(shè)計是整個系統(tǒng)的核心部分。通過對各個部分分別進(jìn)行編程處理來實現(xiàn)，將每部分編為可調(diào)用的函數(shù)，最后統(tǒng)一對函數(shù)進(jìn)行調(diào)用，清晰方便。3.1字母識別整體框圖a BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練過程b BP神經(jīng)網(wǎng)絡(luò)識別過程圖2.1 BP神經(jīng)網(wǎng)絡(luò)識別系統(tǒng)3.2 預(yù)處理及其特征提取方法本文使用傅里葉描述符及其反變

5、換進(jìn)行圖片的二值化、字母輪廓提取，之后進(jìn)行歸一化，將其特征變成1*120的矩陣，之后選取里面的六十個點變?yōu)?*60的矩陣。特征提取程序：function FD=Feature_Building(RGB)%RGB=imread(d:A.bmp);%figure(1),inshow(RGB)B=outline(RGB);%figure(2)%subplot(221),draw_outline(B);%title(outline of object);m,n=size(B);FD=fsd(B,30,m,4);其中outline、fsd為傅里葉描述及其反變換程序。outline程序：%Function

6、 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);fudgeFactor = .5;BWs = edge(I,sobel, threshold * fudgeFactor);%Step 3: Dilate the imagese90 = strel(line, 3,

7、90);se0 = strel(line, 3, 0);BWsdil = imdilate(BWs, se90 se0);%Step 4: Fill interior gapsBWdfill = imfill(BWsdil, 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

8、(BWfinal,seD);bw = bwareaopen(BWfinal,30); % % fill a gap in the pens capB,L = bwboundaries(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ù)，即160=60個輸入神經(jīng)元。3.3.2 隱含層神經(jīng)元個數(shù)的確定隱層節(jié)點數(shù)對網(wǎng)絡(luò)的學(xué)習(xí)和計算特性具有非常重要的影響，是該網(wǎng)絡(luò)結(jié)構(gòu)成敗的關(guān)鍵。若隱層節(jié)點數(shù)過少，則網(wǎng)絡(luò)難以處理復(fù)雜的問題；但若隱層節(jié)點數(shù)過多

9、，則將使網(wǎng)絡(luò)學(xué)習(xí)時間急劇增加，而且還可能導(dǎo)致網(wǎng)絡(luò)學(xué)習(xí)過度，抗干擾能力下降。本文根據(jù)實際的實驗，確定隱含層神經(jīng)元的個數(shù)為15個。3.3.3 輸出層神經(jīng)元個數(shù)的確定因為要識別26個英文大寫字母，因此輸出選擇261的矩陣，即輸出層神經(jīng)元的個數(shù)為26個。當(dāng)26個字母輸入神經(jīng)網(wǎng)絡(luò)后，在對應(yīng)的位置上輸出1，其他位置上輸出零。當(dāng)網(wǎng)絡(luò)進(jìn)入識別過程時，哪個位置上輸出的期望值最大，認(rèn)為識別出的是這個位置上的字母。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.L

10、W2,1*0.01;net.b2=net.b2*0.01;其中minma(P)為神經(jīng)網(wǎng)絡(luò)的對它的60個輸入元素的最大值和最小值的限制。P為訓(xùn)練樣本集合。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ù)。訓(xùn)練函數(shù)采用trainlm 。3.4 BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練3.4.1 訓(xùn)練樣本集合和目標(biāo)值集合字母圖片歸一化后的圖像為601的矩陣，用6026的矩陣形成一個訓(xùn)練樣本；目標(biāo)矢量是希望每一個數(shù)字輸入神經(jīng)網(wǎng)絡(luò)后在輸出神經(jīng)元對應(yīng)的位置上為1,其他的位置為0。為此取目標(biāo)矢量為對角線上為1的2626的單位陣 ,用matl

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

12、P網(wǎng)絡(luò)訓(xùn)練流程圖：使用第一組樣本進(jìn)行訓(xùn)練的結(jié)果：TRAINLM, Epoch 0/5000, SSE 169.303/0.05, Gradient 39.0748/1e-010TRAINLM, Epoch 20/5000, SSE 9.07917/0.05, Gradient 0.647529/1e-010TRAINLM, Epoch 40/5000, SSE 5.45171/0.05, Gradient 0.465044/1e-010TRAINLM, Epoch 60/5000, SSE 3.85999/0.05, Gradient 1.13736/1e-010TRAINLM, Epoch

13、80/5000, SSE 3.37108/0.05, Gradient 0.970379/1e-010TRAINLM, Epoch 100/5000, SSE 1.43394/0.05, Gradient 0.27961/1e-010TRAINLM, Epoch 120/5000, SSE 1.13878/0.05, Gradient 0.661835/1e-010TRAINLM, Epoch 140/5000, SSE 0.561939/0.05, Gradient 0.497918/1e-010TRAINLM, Epoch 160/5000, SSE 0.537153/0.05, Grad

14、ient 0.0963243/1e-010TRAINLM, Epoch 180/5000, SSE 0.518194/0.05, Gradient0.00990168/1e-010TRAINLM, Epoch 200/5000, SSE 0.461637/0.05, Gradient 11.4576/1e-010TRAINLM, Epoch 206/5000, SSE 0.0350697/0.05, Gradient 0.265104/1e-010TRAINLM, Performance goal met.可見經(jīng)過206次訓(xùn)練后，網(wǎng)絡(luò)誤差達(dá)到要求，誤差曲線如下圖：使用第二組樣本進(jìn)行訓(xùn)練的結(jié)果：

15、TRAINLM, Epoch 0/5000, SSE 168.635/0.05, Gradient 33.7987/1e-010TRAINLM, Epoch 20/5000, SSE 3.28669/0.05, Gradient 40.5407/1e-010TRAINLM, Epoch 32/5000, SSE 0.0441687/0.05, Gradient 0.0844925/1e-010TRAINLM, Performance goal met.可見經(jīng)過26次訓(xùn)練之后，網(wǎng)絡(luò)誤差達(dá)到要求。誤差曲線如下圖所示：3.5 字母的識別以上所介紹為網(wǎng)絡(luò)的學(xué)習(xí)期，學(xué)習(xí)過程結(jié)束后，網(wǎng)絡(luò)進(jìn)入工作期，即可以

16、進(jìn)行字母的識別。單一字母識別程序如下：RGB=imread(D:Program FilesMATLAB71work新建文件夾 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；識別為B，則b為2，以此類推。識別結(jié)果如下：a=0.8316b=1可見能夠正確識別字母A。RGB=imread(

17、D:Program FilesMATLAB71work新建文件夾 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可見也能夠正確的識別字母B。本文使用兩組樣本進(jìn)行BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練，一組樣本進(jìn)行字母的識別。識別程序及其結(jié)

18、果如下：load(D:Program FilesMATLAB71work新建文件夾 1index1.mat) for i=1:26 RGBi=imread(D:Program FilesMATLAB71work新建文件夾 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 =4a =0.6114b =5a =0.9755b =6a

19、 =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由以上結(jié)果可知：識別了22個字母，有四個字母未被正確識

20、別（N P S U）。為了使識別準(zhǔn)確率更高，訓(xùn)練更多的樣本，盡量選擇一些略帶有噪聲的圖片，識別時準(zhǔn)確率更高。程序清單BP網(wǎng)絡(luò)訓(xùn)練程序：clcclearload(D:Program FilesMATLAB71work新建文件夾 1index.mat) for i=1:52 RGBi=imread(D:Program FilesMATLAB71work新建文件夾 1,indexi); FDi=Feature_Building(RGBi); FDi=reshape(FDi,1,120); FDi=FDi(1:2:120); endP1=(FD1) (FD2) (FD3) (FD4) (FD5) (F

21、D6) (FD7) (FD8) (FD9) (FD10) (FD11) (FD12) (FD13) (FD14) (FD15) (FD16) (FD17) (FD18) (FD19) (FD20) (FD21) (FD22) (FD23) (FD24) (FD25) (FD26);P2=(FD27) (FD28) (FD29) (FD30) (FD31) (FD32) (FD33) (FD34) (FD35) (FD36) (FD37) (FD38) (FD39) (FD40) (FD41) (FD42) (FD43) (FD44) (FD45) (FD46) (FD47) (FD48) (F

22、D49) (FD50) (FD51) (FD52);%P=P1;P2;T=eye(26);S1=15;S2=26;for n=1:2 %學(xué)習(xí)期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,P

23、n,T);end識別程序：load(D:Program FilesMATLAB71work新建文件夾 1index1.mat) for i=1:26 RGBi=imread(D:Program FilesMATLAB71work新建文件夾 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),

24、inshow(RGB) B=outline(RGB);%figure(2)%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 thr

25、eshold = 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, holes); %Step 5: Remove connected objects on borderB

26、Wnobord = 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 pens capB,L = bwboundaries(bw,noholes);outline = B1;function rFSDs = fsd(outline,H,b,bN)% Forward

27、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 with% iNoOfHarmonicsAnalyse elements.% The output FSDs will be normalised

28、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 %H = iNoOfHarmonicsAnalyse %b = bNormaliseSizeState %m n = size(outline), b

29、= m; %bN = bNormaliseOrientationState rTwoNPi = (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(iNoOfPoints+1,1); % pre-allocate some arraysrDeltaY = zeros(iNoOfPoints+1,1)

30、;rDeltaT = zeros(iNoOfPoints+1,1); for iCount = 2 : iNoOfPoints + 1 rDeltaX(iCount-1) = outline(iCount,1) - outline(iCount-1,1); rDeltaY(iCount-1) = outline(iCount,2) - outline(iCount-1,2);end % Calculate time differences from point to point - actually distances, but we are% carrying on the fiction

31、of a point running around the closed figure at constant speed. % We are analysing the projections on to the x and y axes of this points path around the figurefor iCount = 1 : iNoOfPoints rDeltaT(iCount) = sqrt(rDeltaX(iCount)2) + (rDeltaY(iCount)2);endcheck = (rDeltaT = 0); % remove zeros from rDelt

32、aT, 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 any pointrTime(1) = 0;for iCount = 2 : iNoOfPoints + 1 rTime(iCount) = rTime(iCount - 1) +

33、rDeltaT(iCount-1);end rPeriod = rTime(iNoOfPoints+1); % rPeriod defined for readability % calculate the A-sub-0 coefficientrSum1 = 0;for iP = 2 : iNoOfPoints + 1 rSum2 = 0; rSum3 = 0; rInnerDiff = 0; % calculate the partial sums - these are 0 for iCount = 1 if iP 1 for iJ = 2 : iP-1 rSum2 = rSum2 +

34、rDeltaX(iJ-1); rSum3 = rSum3 + rDeltaT(iJ-1); end rInnerDiff = rSum2 - (rDeltaX(iP-1) / rDeltaT(iP-1) * rSum3); end rIncr1 = (rDeltaX(iP-1) / (2*rDeltaT(iP-1)*(rTime(iP)2-rTime(iP-1)2) + rInnerDiff*(rTime(iP)-rTime(iP-1); rSum1 = rSum1 + rIncr1;end rFSDs(1,1) = (1 / rPeriod) * rSum1) + outline(1,1);

35、 % store A-sub-0 in output FSDs array - this array will be 4 x iNoOfHarmonicsAnalyse % calculate the a-sub-n coefficientsfor iHNo = 2 : H rSum1 = 0; for iP = 1 : iNoOfPoints rIncr1 = (rDeltaX(iP) / rDeltaT(iP)*(cos(rTwoNPi(iHNo-1)*rTime(iP+1)/rPeriod) - cos(rTwoNPi(iHNo-1)*rTime(iP)/rPeriod); rSum1

36、= rSum1 + rIncr1; end rFSDs(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 coefficientsfor iHNo = 2 : H rSum1 = 0; for iP = 1 : iNoOfPoints rIncr1 = (rDeltaX(iP) / rDeltaT(iP)*(sin(rTwoNPi(iHNo-1)

37、*rTime(iP+1)/rPeriod) - sin(rTwoNPi(iHNo-1)*rTime(iP)/rPeriod); rSum1 = rSum1 + rIncr1; end rFSDs(2,iHNo) = (rPeriod / rTwoNSqPiSq(iHNo-1) * rSum1;end % foriHNo = 1 :. % calculate the C-sub-0 coefficientrSum1 = 0;for iP = 2 : iNoOfPoints + 1 rSum2 = 0; rSum3 = 0; rInnerDiff = 0; % calculate the part

38、ial sums - these are 0 for iCount = 1 if iP 1 for iJ = 2 : iP-1 rSum2 = rSum2 + rDeltaY(iJ-1); rSum3 = rSum3 + rDeltaT(iJ-1); end rInnerDiff = rSum2 - (rDeltaY(iP-1) / rDeltaT(iP-1) * rSum3); end rIncr1 = (rDeltaY(iP-1) / (2*rDeltaT(iP-1)*(rTime(iP)2-rTime(iP-1)2) + rInnerDiff*(rTime(iP)-rTime(iP-1)

39、; rSum1 = rSum1 + rIncr1;end rFSDs(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 coefficientsfor iHNo = 2 : H rSum1 = 0; for iP = 1 : iNoOfPoints rIncr1 = (rDeltaY(iP) / rDeltaT(iP)*(cos(rTwoN

40、Pi(iHNo-1)*rTime(iP+1)/rPeriod) - cos(rTwoNPi(iHNo-1)*rTime(iP)/rPeriod); rSum1 = rSum1 + rIncr1; end rFSDs(3,iHNo) = (rPeriod / rTwoNSqPiSq(iHNo-1) * rSum1;end % foriHNo = 1 :. rFSDs(4,1) = 0; % there is no 0th order sine coefficient % calculate the D-sub-n coefficientsfor iHNo = 2 : H rSum1 = 0; f

41、or iP = 1 : iNoOfPoints rIncr1 = (rDeltaY(iP) / rDeltaT(iP)*(sin(rTwoNPi(iHNo-1)*rTime(iP+1)/rPeriod) - sin(rTwoNPi(iHNo-1)*rTime(iP)/rPeriod); rSum1 = rSum1 + rIncr1; end rFSDs(4,iHNo) = (rPeriod / rTwoNSqPiSq(iHNo-1) * rSum1;end % foriHNo = 1 :. % the non-normalised coefficients are now in rFSDs%

42、if we want the normalised ones, this is where it happensif (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 axis of % the first harmonic ellipse rFSDsTemp = rFSDs; rTheta1 = 0.5 * atan(2 * (rFSDs

43、Temp(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 point for iHNo = 1 : H rStarFSDs(1,iHNo) = cos(iHNo-1) * rTheta1) * rFSDsTemp(1,iHNo) + sin(iHNo-1) * rTheta1) * rFSDsTemp(2,iHNo); rStarFSDs(2,