This appendix contains the program list to compute the Fourier coefficient of a periodic nonsinusoidal function. An output example is given for a rectangular wave.
Following is the program list (in FORTRAN) to compute the Fourier terms (DC offset, even and odd harmonics) of a periodic nonsinusoidal function:
C | PROGRAM ANAFOR |
C | test, square wave |
COMMON NW,DEG(361),RAD(361),XV(361),YVALUE(361) | |
COMMON NPTS,NPS,INT,INC,NH,GA,PI,WW,OM,NC,RF,ER,EL | |
COMMON ERR,ITMAX | |
DIMENSION Y(30,365),SCALE(30),YMAX(30),X(365) | |
DIMENSION CURVT(100),CURVC(100),CURVS(100),CURVP(100) | |
DIMENSION AMP(30,100), ANG(30,100) | |
C | READ DATA FROM DATA FILE. |
OPEN(5,FILE=’INDUCTOR.DAT’,STATUS=’OLD’, | |
ACCESS=’SEQUENTIAL’) | |
C | REWIND 5 |
C******************* | |
NSETS = 1 | |
NPTS = 42 | |
C*************** | |
SCALE(1) = 1.0 | |
C | CLOSE(UNIT = 5) |
C | OPEN(3,FILE=’FODA.DAT’,STATUS=’OLD’, ACCESS=’SEQUENTIAL’) |
C | REWIND 3 |
C | DO 5 J = 1,NSETS |
C 5 | READ(3,132) (Y(J,I),I = 1,NPTS + 4) |
132 | FORMAT(E11.4,2X,E11.4,2X,E11.4,2X,E11.4,2X,E11.4) |
C | Input of data |
Do 1320 J = 1,20 | |
1320 | Y(1,J) = 1.0 |
Do 1321 J = 21,41 | |
1321 | Y(1,J) = -1.0 |
Y(1,42) = 1.0 | |
C | ECHO THE DATA ON THE SCREEN. |
DO 10 J = 1,NSETS | |
10 | WRITE(0,1000) J,SCALE(J),(Y(J,I),I = 1,NPTS) |
1000 | FORMAT(//‘ SET NO.’,I4,‘ SCALE FACTOR:’ |
*,F16.10/ ‘Y-COORDINATES BEFORE SCALING:’ (/5E11.4)) | |
C | SUBTRACT DC COMPONENT AND MULTIPLY THE SCALE FACTORS |
DO 30 J = 1,NSETS | |
YMAX(J) = 0. | |
DC = 0.0 | |
DO 37 K = 1,NPTS-1 | |
37 | DC = DC + Y(J,K) |
DC = DC/(NPTS-1) | |
WRITE(0,*) ‘DC=’, DC | |
DO 20 I = 1,NPTS | |
Y(J,I) = SCALE(J)*(Y(J,I)-DC) | |
20 | IF (ABS(Y(J,I)) .GT. YMAX(J)) YMAX(J) = ABS(Y(J,I)) |
30 | CONTINUE |
C | WRITE THE SCALED DATA |
DO 40 J = 1,NSETS | |
40 | WRITE(0,1010) J,YMAX(J),(Y(J,I),I = 1,NPTS) |
1010 | FORMAT(// ‘SET NO.’,I4,‘ PEAK ABSOLUTE VALUE:’, |
*F16.10/‘ Y-COORDINATES AFTER SCALING:’(/5E11.4)) | |
DO 50 I = 1,NPTS | |
50 | X(I) = FLOAT(I-1) |
NSPACE = NPTS-1 | |
DEGMULT = 360.0/NSPACE | |
DO 19 J = 1,NPTS | |
X(J) = X(J)*DEGMULT | |
19 | IF (J.EQ.NPTS) X(J) = X(J) + 0.001 |
C | OPEN(10,FILE = ‘FODA01.DAT’,STATUS = ‘OLD’) |
WRITE (0,*) ‘THE NUMBER OF KNOWN POINT (NPTS) IS’,NPTS | |
WRITE (0,*) ‘THE NUMBER OF DATA SETS (NSETS) IS’ ,NSETS | |
WRITE (0,*) ‘THE KNOWN POINTS OF X AND F(X) ARE : ’ | |
C | THIS PART IS BASED ON THE ALGORITHMS ON PAGE 119 |
C | OF VANDERGRAFT. |
OM = 314.1592654 | |
GA = 0.017453292 | |
PI = 3.141593 | |
INTF = 42 | |
INT = 42 | |
INC = 5 | |
NW = 1 | |
NH = 29 | |
RADIN = PI/20.5 | |
DEGIN = 180./20.5 | |
T1 = -RADIN | |
T2 = -DEGIN | |
DO 22 I = 1,INT | |
T1 = T1 + RADIN | |
RAD(I) = T1 | |
T2 = T2 + DEGIN | |
DEG(I) = T2 | |
22 | CONTINUE |
Write(0,*) (RAD(I), I = 1,NPTS) | |
Write(0,*) (DEG(I), I = 1,NPTS) | |
DO 1 J = 1,NH + 1 | |
CURVC(J) = 0.0 | |
CURVS(J) = 0.0 | |
CURVT(J) = 0.0 | |
CURVP(J) = 0.0 | |
1 | CONTINUE |
DO 17 I = 1,361 | |
XV(I) = FLOAT(I)-1.0 | |
17 | CONTINUE |
DO 12 J = 1,NSETS | |
NC = J | |
DO 13 I = 1,361 | |
YVALUE(I) = Y(J,I) | |
13 | CONTINUE |
YMAX(J) = 0.0 | |
DO 137 I = 1,361 | |
IF(YVALUE(I) .GE. YMAX(J)) YMAX(J) = YVALUE(I) | |
137 | CONTINUE |
WRITE(0,*) ‘J=’,J | |
WRITE(0,*) ‘YMAX(J)=’,YMAX(J) | |
WRITE(0,1013) J,YMAX(J),(YVALUE(I),I = 1,NPTS) | |
1013 | FORMAT(// ‘SET NO.’,I4, ‘PEAK ABSOLUTE VALUE:’, |
*F16.10/ ‘Y-COORDINATES AFTER SCALING:’(/5E11.4)) | |
CALL HARM (YVALUE,CURVC,CURVS,CURVT,CURVP) | |
DO 133 I = 1,NH | |
133 | AMP(J,I) = CURVT(I) |
DO 144 I = 1,NH | |
144 | ANG(J,I) = CURVP(I) |
DO 14 K = 1,NH | |
CURVP(K) = CURVP(K)/GA | |
14 | CONTINUE |
WRITE(0,*) ‘********MAIN**************’ | |
WRITE(0,*) ‘CURVT(‘,NH,J,’)=:’ | |
WRITE(0,100)(CURVT(K),K = 1,NH) | |
WRITE(0,*) ‘THE COSINE BASED PHASOR ANGLES ARE:’ | |
WRITE(0,100)(-CURVT(K)-90.,K = 1,NH) | |
WRITE(0,*) ‘THE MAGNITUDE OF THE HARMONICS AS A’, | |
* ‘PERCENTAGE OF THE FUNDAMENTAL:’ | |
IF (CURVT(1) .EQ.0.00000) THEN | |
WRITE(0,100) (0.0, K = 1,NH) | |
ELSE | |
WRITE(0,100)(100.*CURVT(K)/CURVT(1),K = 1,NH) | |
END IF | |
12 | CONTINUE |
100 | FORMAT (6(1X,1PE11.4)) |
ENDFILE(10) | |
REWIND 10 | |
CLOSE(10) | |
THDI = 0.0 | |
Do 1111 K = 2,NH | |
1111 | THDI = THDI+(CURVT(K))*(CURVT(K)) |
THDIFINAL = (SQRT(THDI))/CURVT(1) | |
WRITE(0,*) THDIFINAL | |
RMSIWITHOUT = 0 | |
Do 1112 K = 1,NH | |
1112 | RMSIWITHOUT = RMSIWITHOUT+(CURVT(K))*(CURVT(K)) |
RMSIWITHOUTF = (SQRT(RMSIWITHOUT))/10. | |
WRITE(0,*) RMSIWITHOUTF | |
RMSIWITH = (DC)*(DC) | |
Do 1113 K = 1,NH | |
1113 | RMSIWITH = RMSIWITH+(CURVT(K))*(CURVT(K)) |
RMSIWITHF = (SQRT(RMSIWITH))/10. | |
WRITE(0,*) RMSIWITHF | |
STOP | |
END | |
SUBROUTINE HARM (VALUE,COMPC,COMPS,COMPT,COMPP) | |
COMMON NW,DEG(361),RAD(361),XV(361),YVALUE(361) | |
COMMON NPTS,NPS,INT,INC,NH,GA,PI,WW,OM,NC,RF,ER,EL | |
COMMON ERR,ITMAX | |
DIMENSION VALUE(361), COMPC(100),COMPS(100),COMPT(100) | |
DIMENSION COMPP(100),VALUEC(361) | |
CHECK = 1.0 | |
WRITE(0,1013) (VALUE(I),I = 1,NPTS) | |
1013 | FORMAT(// ‘ Y-COORDINATES AFTER SCALING:’(/5E11.4)) |
P = SQRT(0.6) | |
A1 = 5./9. | |
A2 = 8./9. | |
FK1 = 0.5*P*(1. + P) | |
FK2 = (1. + P)*(1.-P) | |
FK3 = 0.5*P*(P-1.) | |
DO 44 J = 1,NH | |
DL1 = 0.0 | |
DL2 = 0.0 | |
H = FLOAT(J) | |
DO 33 I = 1,INT-2,2 | |
T1 = H*(RAD(I)) | |
T2 = H*(RAD(I + 1)) | |
T3 = H*(RAD(I + 2)) | |
X1 = FK1*T1 + FK2*T2 + FK3*T3 | |
X2 = T2 | |
X3 = FK3*T1 + FK2*T2 + FK1*T3 | |
C1 = VALUE(I) | |
C2 = VALUE(I + 1) | |
C3 = VALUE(I + 2) | |
C4 = FK1*C1 + FK2*C2 + FK3*C3 | |
C6 = FK3*C1 + FK2*C2 + FK1*C3 | |
D1 = C4*COS(X1) | |
D2 = C2*COS(X2) | |
D3 = C6*COS(X3) | |
D4 = C4*SIN(X1) | |
D5 = C2*SIN(X2) | |
D6 = C6*SIN(X3) | |
DL1 = DL1 + A1*(D1 + D3) + A2*D2 | |
DL2 = DL2 + A1*(D4 + D6) + A2*D5 | |
33 | CONTINUE |
COMPC(J) = DL1*2./41. | |
COMPS(J) = DL2*2./41. | |
COMPT(J) = SQRT(COMPC(J)**2 + COMPS(J)**2) | |
If (COMPS(J).EQ.0.00000) THEN | |
COMPP(J) = 0.0 | |
ELSE | |
COMPP(J) = -ATAN(COMPC(J)/COMPS(J)) | |
END IF | |
IF (COMPS(J) .LT. 0.0) COMPP(J) = COMPP(J) + PI | |
44 | CONTINUE |
WRITE(0,*) ‘*******SUBROUTINE*********’ | |
WRITE(0,*) ‘COMPT(‘,NH,J,’)=:’ | |
WRITE(0,100)(COMPT(K),K = 1,NH) | |
WRITE(0,*) ‘THE COSINE BASED PHASOR ANGLES ARE:’ | |
WRITE(0,100)(-COMPT(K)-90.,K = 1,NH) | |
WRITE(0,*) ‘THE MAGNITUDE OF THE HARMONICS AS A’, | |
* ‘PERCENTAGE OF THE FUNDAMENTAL:’ | |
100 | FORMAT (6(1X,1PE11.4)) |
NWW = 1*NW | |
IF ((CHECK.EQ.1) .AND. (NC .EQ. NWW)) THEN | |
EMDIF = 0 | |
DO 51 I = 1,INT | |
VALUEC(I) = 0. | |
DO 52 J = 1,NH,1 | |
H = FLOAT(J) | |
VALUEC(I) = VALUEC(I) + COMPC(J)*COS(H*RAD(I))+ | |
*COMPS(J)*SIN(H*RAD(I)) | |
52 | CONTINUE |
DIF = ABS(VALUE(I)-VALUEC(I)) | |
IF (DIF.GT.EMDIF) THEN | |
EMDIF = DIF | |
DIFN = FLOAT(I) | |
ENDIF | |
51 | CONTINUE |
WRITE(0,*) ‘***CHECKING FOURIER ANALYSIS*’ | |
WRITE(0,*) | |
WRITE(0,*) ‘VALUE(361) AT NC=’,NC | |
WRITE(0,1001) (VALUE(I), I = 1,42) | |
WRITE(0,1001) (VALUEC(I), I = 1,42) | |
1001 | FORMAT (10(1X,1PE11.4)) |
WRITE(0,*) ‘THE MAXIMUM DIFFERENCE IS :’,EMDIF, ‘AT I=’, DIFN | |
ENDIF | |
RETURN | |
END |
The output of the Fourier program for a square wave is as follows:
SET NO. 1 SCALE FACTOR: 1.0000000000Y-COORDINATES BEFORE SCALING:0.1000E + 01 0.1000E + 01 0.1000E + 010.1000E + 01 0.1000E + 010.1000E + 01 0.1000E + 01 0.1000E + 01 0.1000E + 01 0.1000E + 010.1000E + 01 0.1000E + 01 0.1000E + 01 0.1000E + 01 0.1000E + 010.1000E + 01 0.1000E + 01 0.1000E + 01 0.1000E + 01 0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01-0.1000E + 01 0.1000E + 01 DC = -2.43902E-02SET NO. 1 PEAK ABSOLUTE VALUE: 1.0243902206 Y-COORDINATES AFTER SCALING:0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 010.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 010.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 010.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00 0.1024E + 01THE NUMBER OF KNOWN POINT (NPTS) IS 42THE NUMBER OF DATA SETS (NSETS) IS 1THE KNOWN POINTS OF X AND F(X) ARE :0. 0.153248 0.306497 0.459745 0.612994 0.766242 0.919491 1.07 1.22599 1.37924 1.53248 1.68573 1.83898 1.99223 2.14548 2.29873 2.45198 2.60522 2.75847 2.91172 3.06497 3.21822 3.37147 3.52472 3.67796 3.83121 3.98446 4.13771 4.29096 4.44421 4.59745 4.75070 4.90395 5.05720 5.21045 5.36370 5.51694 5.67019 5.82344 5.97669 6.12994 6.283190. 8.78049 17.5610 26.3415 35.1220 43.9024 52.6829 61.4 70.2439 79.0244 87.8049 96.5854 105.366 114.146 122.927 131.707 140.488 149.268 158.049 166.829 175.610 184.390 193.171 201.951 210.732 219.512 228.293 237.073 245.854 254.634 263.415 272.195 280.976 289.756 298.537 307.317 316.098 324.878 333.659 342.439 351.220 360.000 J = 1YMAX(J) = 1.02439SET NO. 1 PEAK ABSOLUTE VALUE: 1.0243902206 Y-COORDINATES AFTER SCALING:0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 010.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 010.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 010.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00 0.1024E + 01Y-COORDINATES AFTER SCALING:0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 010.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 010.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 010.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01 0.1024E + 01-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00-0.9756E + 00 0.1024E + 01*************SUBROUTINE***********COMPT( 29 30)=:1.2717E + 00 3.4488E-02 4.1969E-01 3.6450E-02 2.4608E-01 3.8742E-021.6870E-01 4.0489E-02 1.2299E-01 4.1067E-02 9.1722E-02 4.0183E-026.8664E-02 3.7860E-02 5.1330E-02 3.4436E-02 3.8801E-02 3.0685E-023.1712E-02 3.8456E-02 3.7949E-02 2.5107E-02 2.4685E-02 2.1601E-022.5122E-02 2.0924E-02 2.7380E-02 2.2261E-02 3.0684E-02THE COSINE BASED PHASOR ANGLES ARE:-9.1272E + 01 -9.0034E + 01 -9.0420E + 01 -9.0036E + 01 -9.0246E + 01 -9.0039E + 01-9.0169E + 01 -9.0040E + 01 -9.0123E + 01 -9.0041E + 01 -9.0092E + 01 -9.0040E + 01-9.0069E + 01 -9.0038E + 01 -9.0051E + 01 -9.0034E + 01 -9.0039E + 01 -9.0031E + 01-9.0032E + 01 -9.0038E + 01 -9.0038E + 01 -9.0025E + 01 -9.0025E + 01 -9.0022E + 01-9.0025E + 01 -9.0021E + 01 -9.0027E + 01 -9.0022E + 01 -9.0031E + 01THE MAGNITUDE OF THE HARMONICS AS A PERCENTAGE OF THE FUNDAMENTAL:*** CHECKING FOURIER ANALYSIS*VALUE(361) AT NC = 11.0244E + 00 1.0244E + 00 1.0244E + 00 1.0244E + 00 1.0244E + 00 1.0244E + 001.0244E + 001.0244E + 00 1.0244E + 00 1.0244E + 00 1.0244E + 00 1.0244E + 00 1.0244E + 001.0244E + 00-9.7561E-01 -9.7561E-01 -9.7561E-01 -9.7561E-01 -9.7561E-01 -9.7561E-01-9.7561E-01-9.7561E-01 -9.7561E-01 -9.7561E-01 -9.7561E-01 -9.7561E-01 -9.7561E-01-9.7561E-01-9.7561E-01 1.0244E + 005.0335E-01 1.1156E + 00 9.8548E-01 1.0407E + 00 9.6388E-01 1.0831E + 009.720 E-019.5716E-01 1.0451E + 00 9.9427E-01 1.0748E + 00 9.3542E-01 1.0637E + 001.022 E + 00-8.8539E-01 -9.6836E-01 -9.0364E-01 -1.0717E + 00 -9.4528E-01 -9.9542E-01-9.457 E-01-9.5329E-01 -1.0331E + 00 -9.2584E-01 -1.0258E + 00 -9.5068E-01 -1.0254E + 00-9.310 E-01-5.0985E-01 5.0335E-01THE MAXIMUM DIFFERENCE IS : 0.521040AT I = 1.00000***********MAIN******************CURVT( 29 1)=:1.2717E + 00 3.4488E-02 4.1969E-01 3.6450E-02 2.4608E-01 3.8742E-021.6870E-01 4.0489E-02 1.2299E-01 4.1067E-02 9.1722E-02 4.0183E-026.8664E-02 3.7860E-02 5.1330E-02 3.4436E-02 3.8801E-02 3.0685E-023.1712E-02 3.8456E-02 3.7949E-02 2.5107E-02 2.4685E-02 2.1601E-022.5122E-02 2.0924E-02 2.7380E-02 2.2261E-02 3.0684E-02THE COSINE BASED PHASOR ANGLES ARE:-9.1272E + 01 -9.0034E + 01 -9.0420E + 01 -9.0036E + 01 -9.0246E + 01 -9.0039E + 01-9.0169E + 01 -9.0040E + 01 -9.0123E + 01 -9.0041E + 01 -9.0092E + 01 -9.0040E + 01-9.0069E + 01 -9.0038E + 01 -9.0051E + 01 -9.0034E + 01 -9.0039E + 01 -9.0031E + 01-9.0032E + 01 -9.0038E + 01 -9.0038E + 01 -9.0025E + 01 -9.0025E + 01 -9.0022E + 01-9.0025E + 01 -9.0021E + 01 -9.0027E + 01 -9.0022E + 01 -9.0031E + 01THE MAGNITUDE OF THE HARMONICS AS A PERCENTAGE OF THE FUNDAMENTAL:1.0000E + 02 2.7118E + 00 3.3001E + 01 2.8661E + 00 1.9350E + 01 3.0463E + 001.3265E + 01 3.1837E + 00 9.6712E + 00 3.2292E + 00 7.2123E + 00 3.1597E + 005.3992E + 00 2.9770E + 00 4.0362E + 00 2.7078E + 00 3.0510E + 00 2.4128E + 002.4936E + 00 3.0239E + 00 2.9840E + 00 1.9742E + 00 1.9410E + 00 1.6985E + 001.9754E + 00 1.6453E + 00 2.1530E + 00 1.7504E + 00 2.4127E + 00 0.444086 0.139151 0.139172