7.0 PERIOD CREATION : Kurt D. Schwehr - August, 1991 LAST CHANGED: Kurt D. Schwehr - December 15, 1991 7.1 PROGRAM PERIOD.EXE 7.2 SYNOPSIS Generate a Periodogram/Power Spectrum, find the peaks, and generate a graph and a report. 7.3 PURPOSE To analyze robotic telescope data for periodic events (or any time series data). 7.4 DESCRIPTION PERIOD uses an algorithm developed by Scargle which allows the signal power to computed at any frequency. From this power it is possible to determine what the probability that the power represents a real frequency in the data. 7.5 COMPUTER REQUIREMENTS - IBM PC/compatible or greater - 256 Kb or more - VGA graphics adapter - HP DeskJet or LaserJet printer, any model (for hardcopies) 7.6.0 INSTRUCTIONS/USER MANUAL 7.6.1 Preparing a data file and the data. The data file is a file containing whatever data you wish to analyze. It can contain up to 5000 data points and must be in ASCII text format. You may have columns containing other data for each point. Now write down the column number for the time index and magnitude of the data. Write down the number of lines that precede your data. They are usually notes and column titles. Now check that your data columns are separated only Spaces, Commas or TABS. If there are other non numeric characters in your data, the program will only analyze a small portion of your data. It is important that as much data be present as possible. The more periods of the cycle that you are looking for that exist in your data the more likely the program will find it. Unfortunately, as the number of samples increases, the time increases. Large data sets will take a long time to run. 7.6.2 Files required to run. - PERIOD.EXE -> Must be in the current directory or path - TMSRB.FON -> Must be in the current directory - your data file -> Must be in the current directory 7.6.3 Description of a Sample Run Here is the data file, "131156.NIT" used for the sample: - Night Table - Star : 131156 Date : 1991-11-06 Group : 148 Ref Star : 131265 Day # Mag 2 % Err # Mag 3 % Err # Clr 2-3 % Err 346.680 10 5.1374 .44 11 4.6359 .31 11 .7004 .35 347.633 3 5.1284 .11 3 4.6269 .28 3 .7015 .38 348.696 11 5.1162 .22 12 4.6222 .56 11 .6937 .35 350.682 3 5.1135 .11 3 4.6226 .22 3 .6909 .27 351.698 11 5.1302 .16 11 4.6309 .37 11 .6988 .37 353.699 9 5.1360 .31 9 4.6313 .43 8 .7060 .32 354.653 3 5.1322 .86 3 4.6371 .68 3 .6951 .18 355.762 15 5.1188 .57 14 4.6231 .35 14 .6959 .33 356.653 3 5.1194 .18 3 4.6234 .01 3 .6960 .18 357.654 6 5.1302 .22 6 4.6318 .22 6 .6984 .20 359.655 6 5.1390 .23 6 4.6375 .18 6 .7014 .31 360.702 9 5.1306 .35 8 4.6267 .49 9 .7016 .65 361.706 10 5.1216 .35 9 4.6234 .25 10 .6969 .45 364.658 3 5.1416 .57 3 4.6484 1.19 3 .6932 .67 367.745 11 5.1249 .44 11 4.6280 .39 12 .6958 .70 371.662 3 5.1398 .15 3 4.6379 .13 3 .7019 .15 372.665 3 5.1339 .24 3 4.6345 .11 3 .6993 .31 373.663 3 5.1263 .06 3 4.6230 .17 3 .7033 .22 374.677 6 5.1174 .43 6 4.6184 .80 6 .6989 .89 377.695 9 5.1372 .16 9 4.6355 .27 10 .7020 .46 378.737 8 5.1394 .66 8 4.6389 .92 8 .7004 .75 379.667 8 5.1267 .35 7 4.6246 .18 7 .7025 .20 380.714 15 5.1222 .29 15 4.6242 .40 15 .6973 .31 381.728 11 5.1137 .19 11 4.6187 .31 11 .6959 .19 382.663 7 5.1163 .14 7 4.6166 .36 7 .6987 .29 383.681 8 5.1288 .46 7 4.6280 .26 7 .6974 .40 386.666 8 5.1269 .39 7 4.6256 .23 8 .7034 .62 387.710 5 5.1171 .24 5 4.6181 .24 5 .6990 .23 388.684 7 5.1127 .26 8 4.6175 .24 7 .6948 .21 390.689 12 5.1272 .55 11 4.6303 .82 12 .6990 .81 391.737 9 5.1265 .32 10 4.6285 .62 10 .6993 .67 392.707 1 5.1115 .00 1 4.6232 .00 1 .6882 .00 393.753 12 5.1202 .23 12 4.6226 .21 13 .6988 .38 394.696 5 5.1162 .81 5 4.6209 .49 5 .6953 .47 In the sample run: the "!" is followed by comments, computer output is in a different font, and user input is underlined. D:\DATA>PERIOD ! Begin the program Enter filename: 131156.NIT ! Enter your data file Enter Unit Time (EX: Sec,Min,Hr,Day) up to 3 characters: Day Number of text lines to skip in the beginning of your data file: 6 Enter the two column numbers to be read: 1,6 Reading Data . . . Number of Data Points = 34 Maximum Frequency(Nyquist Freq.) = 7.011875E-01 (Cycles/Day) Minimum Frequency = 1.041319E-02 (Cycles/Day) Change Minimum(y/N)? N Change Maximum(y/N) ? N Mean for 2nd col # of data = 4.628 Variation = .000 Enter the subject (up to 60 characters): Star: HD 131156 Filter: Mag3-VIS GRP:148 File:131156.NIT ! Enter whatever title you like Do you want to send the results to the printer(y/N)? N ! Enter 'y' if you want a hardcopy Computing Initial Survey . . . xx% Done ! This is how far through calculation it is Analyzing 3 peaks . . . ! How many peaks it is looking for Here is the graph of the periodogram/power spectrum: Here is the report: Frequency range: .0000 TO .7012 (Cycles/Day) Nyquist Freq.: .7012 Frequency Omega Period Power Probability --------- ----- ------ ----- ----------- .15636 .98 6.40 10.4 .999 .01823 .11 54.85 3.3 .251 .61775 3.88 1.62 2.6 .077 D:\DATA> ! Back to DOS 7.6.4 The report Omega(w) is defined in f(t) = sin w(t) + x. Period is the length of time to complete one cycle. Power is the relative power contained in the signal. Probability is the chance that a signal is real. 7.6.5.0 Important Definitions 7.6.5.1 Minimum Frequency This is the lowest frequency that the power spectrum can distinguish in the data. In order for the equation to detect the frequency it must see at least one full cycle or period of the wave. The more periods in the data, the more likely that the power spectrum will show the wave energy. 7.6.5.2 Nyquist Frequency The Nyquist frequency is the point a which it is no longer possible to detect frequencies above. It is given by: fn = 1 / (Average Time Step). If you need more information on this please refer to Presses' Numerical Recipes or another text on power spectrum analysis. 7.7 COMPILING This program was compiled with MS FORTRAN 5.0 under the large memory model. The command line for compiling was: "FL /FPi PERIOD.FOR POWER.FOR GRAPHICS.LIB" 7.8 ALGORITHM The algorithm used was taken from Horne(1986). Please refer to the paper for more information as this section is completely explained in it. Horne, page 757: For a time series X(tj), where i = 1,2,...,No Px(w) = (1/2) [(Sum12 / Sum2) + (Sum32 / Sum4)] where: tan (2wt) = (sum 2wtj)/(sum cos 2wtj) => sum is from j=1 to j=No sum1 = X(tj) cos w(tj - t) => " sum2 = cos2 w(tj - t) => " sum3 = X(tj) sin w(tj - t) => " sum4 = sin2 w(tj - t) => " For the following equations to produce the correct probability, they must be properly normalized: 1) The data must have the mean value subtracted from is. 2) PN(w) = PX(w)/s2 Ni = Number of Selected Frequencies Ni = -6.362 + 1.193*NumPoints + 0.00098 * NumPoints2 F = Probability that frequency is not real F = 1 - (1 - e (-PN) )Ni Probability = 1 - F 7.9 PROGRAM FLOW This is a description of the program. FILE: PERIOD.FOR - Initialize: Question User Read Data File Calculate constants Prepare Data - Calculate Survey: Get an idea what it looks like - Find Peaks: Survey each peak Find the maximum value for peak - Start Graphics: Find the range and domain Set graphics mode Draw graph and ticks Put on the title Label the axises - Plot Survey: Draw the points of the survey - Plot Peak Surveys: Plot all the points Draw a line at each peak - Print: Send to printer if selected Reset screen to text mode - Report: Show each peak and it's probability Send to printer if selected - Exit: Return to DOS 7.10 BUGS The one serious bug is in reading data files. If the program hits a character between numbers that is not a separator (CR, LF, Space, or TAB) then it thinks that it has hit the end of the data. This is a problem to our work, as we mark certain data points with an asterisk (*). Therefor we must strip them from the data file before running it through PERIOD. The way to see if this is happening to a particular data set is to check the answer given after "Number of data points:." If this value is far lower than you expect, this is probably what is happening. 7.11 SUGGESTIONS FOR FUTURE WORK - Incorporate this program to work automatically in other programs - Set it up for other video displays - I attempted to set up the program to use the LFIT routine from Numerical Recipes. This would find a quadratic approximation (c + cX + cX^2) to a number of samples of the peak. It would then use the derivative of this to find a better approximation of the peak. I was unable to get it working. 7.12 REFERENCES HORNE, J.H. and BALIUNAS S., "A Prescription for Period Analysis of Unevenly Sampled Time Series," The Astrophysical Journal, 302:757-763, 1986 March 15. PRESS, H.P. et. al.,Numerical Recipes, The Art of Scientific Computing (FORTRAN Version), Cambridge Press, 1989. SCARGLE, J. D., "Studies In Astronomical Time Series Analysis. II. Statistical Aspects of Spectral Analysis of Unevenly Space Data," The Astrophysical Journal, 263:835-853, 1982 December 15.