Period detection device, period detection method, and program

Abstract

To improve the accuracy of detecting the correct period while suppressing the increase in the number of candidate periods.SOLUTION: A period detection device comprises a power spectrum calculation part that calculates the power spectrum of astronomical data that is time-series data representing changes in the brightness of celestial bodies, for each given candidate period, a candidate period reset part that resets the candidate period by adding the periods surrounding each candidate period with a high power spectrum as a new candidate period, on the basis of a predetermined number of candidate periods with high power spectrums among power spectra for each candidate period calculated by the power spectrum calculation part, and a period calculation part that detects the period of the celestial body from all candidate periods reset by the candidate period reset part.SELECTED DRAWING: Figure 5

Description

[Technical Field]

[0001] The present invention relates to a period detection device, a period detection method, and a program. [Background technology]

[0002] Variable objects, which are celestial objects whose brightness changes over time, include flare stars, novae, supernovae, Cepheid stars, and RR Lyrae stars. Many variable objects exhibit periodic changes in brightness, and this phenomenon is useful for elucidating the structure and evolution of the universe. Fourier transforms are a typical method for extracting periodicity in data, but Fourier transforms cannot be applied to periodic analysis of celestial objects. This is because observational data of celestial objects contains many missing values ​​and is very sparse.

[0003] Figure 11 shows an example of time-series data for a certain celestial object. The horizontal axis in Figure 11 represents the observation date and time, and the vertical axis represents the brightness (magnitude) of the celestial object. Observation data for celestial objects is often limited by weather conditions or the observation of other celestial objects, so it is inevitable that the data will be non-uniformly sampled, as shown in Figure 11. Furthermore, the observed values ​​are affected by various noises, such as noise from the atmosphere, noise from the observation equipment, and noise from other stars. For these reasons, it is difficult to identify the period of a celestial object using a Fourier transform.

[0004] Conventionally, methods such as the string length method (see, for example, Non-Patent Document 1) and the Lomb-scargle periodogram method (see, for example, Non-Patent Documents 2, 3, and 4) have been proposed for analyzing the periods of other celestial bodies. [Prior art documents] [Non-patent literature]

[0005] [Non-Patent Document 1] D. Clarke, “String / Rope length methods using the Lafler-Kinman statistic”, Astronomy & Astrophysics. [Non-patent document 2] NRLOMB, “LEAST - SQUARES FREQUENCY ANALYSIS OF UNEQUALLY SPACED DATA”, Astrophysics and space science, 1976. [Non-patent document 3] Jacob T. VanderPlas, “Understanding the Lomb-Scargle Periodogram”, The Astrophysical Journal Supplement Series, 236:16 (28pp), May 2018. [Non-patent document 4] JEFFRREY D. SCARGLE, “Studies in astronomical time series analysis. ii. statistical aspects of spectral analysis of unevenly spaced data”, Astrophysical Journal, 1982. Summary of the Invention [Problem to be solved by the invention]

[0006] The stringlength method has the advantage of being able to perform period analysis regardless of the number of missing values. However, the stringlength method has the problem that in order to achieve relatively high accuracy, the number of candidate periods becomes enormous, and the longest period among the candidate periods tends to be judged as the best period. The Lomb-scargle perodogram method is an extension of the Fourier transform to accommodate non-uniform sampling, and is a signal processing technique used in various fields such as the automotive industry and medicine. However, the Lomb-scargle perodogram method is not very accurate for very sparse data.

[0007] All of the above-mentioned methods require the specification of candidate periods (hereinafter referred to as "candidate periods"). However, there is a problem in that the number of candidate periods becomes very large when applying these methods to astronomical object data in which the periods differ for each celestial object.

[0008] In view of the above circumstances, an object of the present invention is to provide a technique that can improve the accuracy of detecting a correct period while suppressing an increase in the number of candidate periods. [Means for solving the problem]

[0009] One aspect of the present invention is a period detection device that includes a power spectrum calculation unit that calculates the power spectrum of celestial body data, which is time-series data that represents changes in the brightness of a celestial body, for each given candidate period; a candidate period resetting unit that resets the candidate periods by adding, as new candidate periods, periods surrounding each of the predetermined number of candidate periods with high power spectra, based on the power spectra for each candidate period calculated by the power spectrum calculation unit.; and a period calculation unit that detects the period of the celestial body from all of the candidate periods reset by the candidate period resetting unit.

[0010] One aspect of the present invention is a period detection method that calculates the power spectrum of celestial body data, which is time-series data that represents changes in the brightness of a celestial body, for each given candidate period, and based on a predetermined number of candidate periods with high power spectra from the calculated power spectra for each candidate period, adds periods surrounding each of the predetermined number of candidate periods with high power spectra as new candidate periods to reset the candidate periods, and detects the period of the celestial body from all of the reset candidate periods.

[0011] One aspect of the present invention is a program for causing a computer to execute the following steps: a power spectrum calculation step of calculating, for each given candidate period, the power spectrum of celestial body data, which is time-series data that represents changes in the brightness of a celestial body; a candidate period resetting step of resetting candidate periods by adding, as new candidate periods, periods surrounding each of a predetermined number of candidate periods with high power spectra, based on the power spectra of the candidate periods calculated by the power spectrum calculation step; and a period calculation step of detecting the period of the celestial body from all of the candidate periods reset in the candidate period resetting step. [Effects of the Invention]

[0012] According to the present invention, it is possible to improve the accuracy of detecting the correct period while suppressing an increase in the number of candidate periods. [Brief explanation of the drawings]

[0013] [Figure 1] FIG. 1 shows the results of the folding method. [Figure 2] This is a diagram (part 1) to explain the problems with the Stringlength method. [Figure 3] This is a diagram (part 2) to explain the problems with the Stringlength method. [Figure 4] FIG. 1 is a diagram for explaining an overview of the present invention. [Figure 5] FIG. 2 is a diagram illustrating a specific example of a period detection device according to the first embodiment. [Figure 6] 4 is a flowchart showing the flow of processing of a period detection method performed by the period detection device in the first embodiment. [Figure 7] 10 is a flowchart showing the flow (part 1) of candidate period resetting processing performed by a candidate period resetting unit in the first embodiment. [Figure 8] 10 is a flowchart showing the flow (part 2) of the candidate period resetting process performed by the candidate period resetting unit in the first embodiment. [Figure 9] FIG. 10 is a diagram showing an algorithm illustrating the procedure of processing performed by a power spectrum calculation unit 110 in the second embodiment. [Figure 10] FIG. 10 is a flowchart showing the processing of the algorithm shown in FIG. 9. [Figure 11] FIG. 1 is a diagram showing an example of time-series data of a certain celestial body. DETAILED DESCRIPTION OF THE INVENTION

[0014] Hereinafter, an embodiment of the present invention will be described with reference to the drawings. In the following explanation, symbols placed above letters in mathematical expressions or functions (hereinafter referred to as "mathematical expressions, etc.") are written before the letters. For example, the symbol "-" placed above letters in mathematical expressions, etc. will be written as "-x" before the letter "x" below.

[0015] Before describing the present invention in detail, a conventional method will be described. [Folding method] The folding method is a technique that is often used to confirm the periodicity of celestial bodies. Details of the folding method are described in Non-Patent Document 1. The folding method is a technique that performs a process called folding, which converts and rearranges the time information of observation data so as to extract phase information for a certain period (hereinafter referred to as the "trial period"). Here, the procedure of the folding method will be explained in detail. First, in the folding method, the time series data of the celestial body is divided into equal intervals at the candidate period. Next, in the folding method, all the divided sections are superimposed onto one another. Specifically, the vertical axis is the brightness, the horizontal axis is the phase, and the candidate period is set to 1, and they are superimposed. The observation date and time of the data of each celestial body is set to t, and the candidate period is set to P c Then, the phase φ of the candidate period after folding is expressed as the following equation (1).

[0016]

number

[0017] An example of the folding method is shown in Figure 1. Figure 1(A) shows the results when the graph is divided by a period different from the correct period (e.g., 2π), while Figure 1(B) shows the results when the graph is divided by the correct period. As shown in Figures 1(A) and 1(B), by dividing the graph by the correct period of 2π and overlaying them, a clean wave of one period is obtained, as shown in Figure 1(B). On the other hand, if the candidate period is far from the correct period, the result will be scattered, like the original graph with non-uniform sampling, as shown in Figure 1(A). As can be seen, with the folding method, if the candidate period is very close to the correct period of the celestial body, the graph will have a shape similar to a wave of one wavelength, resulting in a neat shape. However, if the candidate period is not close to the correct period, the shape will be very scattered.

[0018] [String length method] Next, we will explain the Stringlength method. The Stringlength method is a widely used method for analyzing astronomical periods. The Stringlength method first performs folding on each candidate period. Next, the Stringlength method calculates the sum of the Euclidean distances between adjacent observation points for each observation point on the folded graph. The Stringlength method performs the same calculation on all candidate periods, and determines the period with the smallest sum of the Euclidean distances between adjacent observation points as the best period. The Stringlength method can perform period analysis without being affected by the sparseness of the data observation date and time, making it suitable for handling very sparse astronomical data, as shown in Figure 9. However, the Stringlength method has the problem of determining the longest candidate period as the best.

[0019] Figure 2 is a diagram (part 1) to explain the problems of the string length method. (A) in Figure 2 shows the correct period P c = 0.33306, and (B) in Figure 2 shows a graph with a period P close to a multiple of the correct period for a certain celestial object. c = 0.99904. As shown in Figure 2, the actual correct period (P c =0.33306) than the correct period (period P c = 0.99904) has a smaller sum of Euclidean distances, so the stringlength method may judge a period different from the correct period to be the best.

[0020] In addition, the Stringlength method requires that the candidate period be very close to the correct period in order to obtain a clean waveform, so the number of candidate periods must be increased, which can result in high computational costs. Figure 3 is a diagram (part 2) to explain the problems with the Stringlength method. Figure 3 (A) shows the correct period P for a certain celestial body. c = 0.501086, and (B) in Figure 3 shows the period Pc = 0.501042. As shown in Figure 3, although the difference in period between the graphs shown in Figure 3(A) and Figure 3(B) is only about 0.00004, it can be seen that there are parts in the graph shown in Figure 3(B) that do not form a waveform-like shape. As such, even when there is a minute difference, such as five decimal places, between the correct period and the candidate period, it can be seen that the folding results do not form a neatly organized shape. As a result, the sum of the Euclidean distances for each observation point becomes long, and the Stringlength method may determine that a period different from the correct period is the best.

[0021] [Lomb-Scargle periodogram method] The Lomb-Scargle periodogram is an algorithm used to detect and characterize periodic signals in non-uniformly sampled data, and is used in a variety of fields, including astronomy, automotive, communications, and medicine. k N (N is an integer greater than or equal to 1) observations x k For (k=1,...,N), the Lomb-Scargle periodogram method calculates the power spectrum for frequency (period) using the following equation (2).

[0022]

number

[0023] -x and σ in equation (2) 2 are the mean and variance of the observed values, respectively, and τ is a value that satisfies the following equation (3).

[0024]

number

[0025] The Lomb-Scargle periodogram method is a method that can find the correct period without requiring a huge number of candidate periods and while keeping computational costs low. However, it has the problem that accuracy is significantly reduced due to the influence of data sparseness.

[0026] [Method of the present invention] Next, the method used in the present invention to detect the periods of variable celestial bodies will be described. First, an overview of the present invention will be described using Figure 4. The present invention uses the Lomb-scargle periodogram method and folding. First, the present invention performs period analysis on all candidate periods using the Lomb-scargle periodogram method. Next, the present invention extracts the top K (K is an integer equal to or greater than 2) periods whose power spectra are obtained as a result of period analysis using the Lomb-scargle periodogram method. Note that while Figure 4 shows the top 20 candidate periods (K = 20), this is just an example, and the number of the top K periods to be extracted is not particularly limited. Next, the present invention further analyzes each of the top K extracted periods using folding.

[0027] Specifically, folding is performed on each of the cycles obtained by further dividing the periphery of the top K cycles into smaller parts, thereby folding a cycle closer to the correct solution. Then, in this invention, an approximate curve is drawn for each graph obtained by folding using Polyfit or Group Lasso, and the mean least squares error between the approximate curve and the observed value is calculated. The mean least squares error represents the magnitude of deviation from the approximate curve. Therefore, in this invention, the cycle with the smallest mean least squares error can be determined to be the best cycle.

[0028] The Lomb-scargle periodogram method is used in this invention because it can detect power spectra of candidate periods close to the correct period even when the number of candidate periods is small. However, because astronomical observation data is very sparse, the accuracy is not very high using the Lomb-scargle periodogram method alone. Furthermore, if folding is performed on all candidate periods, the correct period cannot be found unless the number of candidate periods is enormous. Therefore, this invention solves this problem by using the Lomb-scargle periodogram method to narrow down the candidate periods and then searching the surrounding areas of the candidate periods using folding. Furthermore, the use of the mean least squares error is expected to solve the problem expected in the stringlength method (the problem of a certain long period being easily determined to be the best period).

[0029] (First embodiment) [Device configuration] Next, a device that performs the above-mentioned processing will be described. FIG. 5 is a diagram showing a specific example of a period detection device according to the first embodiment. The period detection device 1 is a device that detects the period of a celestial body using celestial body data. The period detection device 1 includes a control unit 10. The control unit 10 includes a central processing unit (CPU), memory, auxiliary storage device, and the like, all connected via a bus, and executes a program. By executing the program, the control unit 10 functions as a device including a power spectrum calculation unit 110, a candidate period resetting unit 120, and a period calculation unit 130. Note that all or part of the functions of the control unit 10 may be realized using hardware such as an application-specific integrated circuit (ASIC), a programmable logic device (PLD), or a field-programmable gate array (FPGA). The program may also be recorded on a computer-readable recording medium. Examples of computer-readable recording media include portable media such as flexible disks, magneto-optical disks, read-only memories (ROMs), and CD-ROMs, as well as storage devices such as hard disks built into computer systems. The program may also be transmitted and received via a telecommunications line.

[0030] The power spectrum calculation unit 110 calculates the power spectrum of each of all candidate periods given for the input celestial body data using the Lomb-scargle periodogram method, and extracts K candidate periods from all the candidate periods in descending order of power spectrum. For example, if there are 100 candidate periods, the power spectrum calculation unit 110 applies the Lomb-scargle periodogram method to each of the 100 candidate periods and extracts K (e.g., 20) candidate periods with the largest power spectrum.

[0031] While the Stringlength method uses a huge number of candidate periods and misses the correct period unless the candidate period is very close to the correct period, the Lomb-scargle periodogram method can find the correct period with a smaller number of candidate periods.Furthermore, since it does not have the problem of determining a long period as the best period, which is present in the Stringlength method, the present invention uses the Lomb-scargle periodogram method as a method for narrowing down the candidate periods.

[0032] The candidate period resetting unit 120 adds neighboring periods of each of the top K candidate periods extracted by the power spectrum calculation unit 110 as new candidate periods. For example, if four neighboring periods are added as new candidate periods for each candidate period, the number of candidate periods is 20 + 20 × 4 = 100. In this way, the candidate period resetting unit 120 resets the candidate periods by adding new candidate periods to the extracted top K candidate periods. The candidate period resetting unit 120 performs folding on each of the top K candidate periods and the added candidate periods. For example, if the number of candidate periods is 100 as described above, the candidate period resetting unit 120 performs folding on each of the 100 candidate periods.

[0033] There are two methods by which the candidate period resetting unit 120 searches for periods surrounding each of the top K periods. The first method divides the area before and after each of the top K periods into equal intervals using a predetermined number of divisions, and performs additional period analysis on each divided period. The second method uses a spline function to create an approximate curve and search for a more optimal period. Because the processing by the candidate period resetting unit 120 can create candidate periods that are very close to the correct period, it can solve the problem with the Lomb-scargle periodogram method, in which periods that are not very close to the correct period for folding cannot be organized into a wavelength shape.

[0034] The period calculation unit 130 calculates the period of the celestial body based on the folding results for each of the candidate periods (the top K candidate periods and the added candidate period) reset by the candidate period resetting unit 120. Specifically, the period calculation unit 130 first draws an approximate curve for each of the multiple graphs obtained as a folding result by the candidate period resetting unit 120. For example, if the candidate period resetting unit 120 performs folding for each of 100 candidate periods, the period calculation unit 130 draws an approximate curve for each of the 100 graphs obtained as a folding result.

[0035] There are two ways to draw an approximation curve. The first method uses a simple moving average, and the second method uses the Group Lasso. The simple moving average is a method of calculating the average value of data within a specific interval. The method using the Group Lasso solves the regression problem defined in the following equation (4) to create an approximation curve using the sum of a cosine function and a sinus function.

[0036]

number

[0037] In equation (4), N represents the number of observation points, and t kwhere m represents the observation date and time, and m represents the number of candidate periods. The method using a simple moving average has significantly lower computational costs than the method using Group Lasso, but the method using Group Lasso is characterized by higher fitting accuracy. After drawing an approximate curve for each of the multiple graphs, the period calculation unit 130 calculates the mean least squares error (mse) between the drawn approximate curve and the graph. The magnitude of the mean least squares error indicates how much the folded graph deviates from the approximate curve. Therefore, a graph folded with a period close to the correct period will have a small mean least squares error because it resembles a wavelength. Conversely, a graph folded with a period not close to the correct period will have a large mean least squares error. The period calculation unit 130 outputs the period with the smallest mean least squares error as the correct period. Note that the period calculation unit 130 may select, as the correct period, a period with a mean least squares error less than a threshold, rather than the period with the smallest mean least squares error. In this case, for example, the period calculation unit 130 may select the period with the second smallest mean least squared error, or the third smallest. As mentioned above, the larger the mean least squared error, the further away from the correct period it is, so the threshold value for comparison with the mean least squared error should not be too large. In the following explanation, we will use an example configuration in which the period calculation unit 130 selects the period with the smallest mean least squared error as the correct period. Like the present invention, the string length method also performs period analysis using a folded graph, but there is a problem in that the longer period is determined to be the best period. However, this problem can be solved in the present invention by using the mean least squared error.

[0038] (Period detection method using period detection device 1) FIG. 6 is a flowchart showing the flow of processing of the period detection method performed by the period detection device 1 in the first embodiment. Assume that astronomical data is input to the period detection device 1 from outside. Note that the astronomical data may not be input from outside but may be stored inside the period detection device 1. The power spectrum calculation unit 110 calculates the power spectrum for each of all given candidate periods for the input celestial data using the Lomb-scargle periodogram method (step S101). Then, based on the calculated power spectrum, the power spectrum calculation unit 110 extracts K candidate periods from all the candidate periods in descending order of power spectrum (step S102). The power spectrum calculation unit 110 outputs information indicating the extracted K candidate periods to the candidate period resetting unit 120.

[0039] The candidate period resetting unit 120 performs a candidate period resetting process based on the information indicating the K candidate periods output from the power spectrum calculation unit 110 (step S103). The candidate period resetting process is a process for resetting candidate periods, and represents a process of adding a new candidate period in addition to the top K candidate periods as described above. The candidate period resetting unit 120 resets the candidate periods by executing the candidate period resetting process. As a result, the number of candidate periods becomes K periods + the number of added candidate periods. The candidate period resetting unit 120 performs folding on each of the reset candidate periods (step S104). This results in a graph for each reset candidate period (a graph like that shown in FIG. 2 or FIG. 3).

[0040] The candidate period resetting unit 120 outputs information indicating the multiple graphs obtained for each of the reset candidate periods to the period calculation unit 130. The period calculation unit 130 detects the correct period based on the information indicating the multiple graphs output from the candidate period resetting unit 120 (step S105). Specifically, the period calculation unit 130 first draws an approximation curve for each of the multiple graphs. Next, the period calculation unit 130 calculates the mean least square error (mse) between the approximation curve and the graph for each of the multiple graphs. Then, the period calculation unit 130 detects the period with the smallest mean least square error as the correct period.

[0041] Fig. 7 is a flowchart showing the flow (part 1) of the candidate period resetting process performed by the candidate period resetting unit 120 in the first embodiment. Fig. 7 describes a case where the candidate period resetting unit 120 resets the candidate period by dividing it into equal intervals (first method). The candidate period resetting unit 120 receives as input information indicating K candidate periods output from the power spectrum calculation unit 110 (step S201). The candidate period resetting unit 120 repeatedly executes the processes from loop Ls1 to loop Le1. In loop Ls1 to loop Le1, periods obtained when i is changed from 1 to K are added as candidate periods. The candidate period resetting unit 120 executes the processes from loop Ls1 to loop Le1 until it has added K candidate periods. The candidate period resetting unit 120 assigns 1 to i and extracts the period immediately before and the period immediately after the highest i-th candidate period (step S202).

[0042] When i=1, the candidate period resetting unit 120 extracts the period immediately preceding and the period immediately following the first highest candidate period. The candidate period resetting unit 120 further divides the space between the period immediately preceding and the period immediately following the extracted i-th highest candidate period into l (l is an integer equal to or greater than 1) equal intervals (step S203). For example, when l=4, the candidate period resetting unit 120 divides the space between the period immediately preceding and the period immediately following the extracted i-th highest candidate period into four equal intervals. As a result, four new candidate periods are set for each candidate period. The candidate period resetting unit 120 changes i from 1 to K and adds the newly set candidate periods for each of the first to K-th highest candidate periods to the set C of candidate periods. Thereafter, the candidate period resetting unit 120 executes the processes from loop Ls1 to loop Le1 to obtain the final set C of candidate periods (step S204). The set C of candidate periods is the number of candidate periods described above (K+the number of added candidate periods).

[0043] Fig. 8 is a flowchart showing the flow (part 2) of the candidate period resetting process performed by the candidate period resetting unit 120 in the first embodiment. Fig. 8 describes a case where the candidate period resetting unit 120 resets the candidate period using a spline function (second method). The candidate period resetting unit 120 receives as input information indicating the K candidate periods and information indicating the power spectrum output from the power spectrum calculation unit 110 (step S301). The candidate period resetting unit 120 adds the K candidate periods to a set C of candidate periods (step S302). The candidate period resetting unit 120 calculates an approximate curve for the relationship between the K candidate periods and the power spectrum using a spline function (step S303). The candidate period resetting unit 120 adds periods at which the spline function takes a maximum value to the set C of candidate periods (step S304). The candidate period resetting unit 120 then performs this process on the K candidate periods to obtain the final set C of candidate periods (step S305).

[0044] The period detection device 1 configured as described above includes a power spectrum calculation unit 110 that calculates the power spectrum of celestial body data, which is time-series data that represents changes in the brightness of a celestial body, for each given candidate period; a candidate period resetting unit 120 that resets the candidate periods by adding, as new candidate periods, periods surrounding each of the predetermined number of candidate periods with high power spectra, based on the predetermined number of candidate periods with high power spectra, out of the power spectra for each candidate period calculated by the power spectrum calculation unit 110; and a period calculation unit 130 that detects the period of the celestial body from all of the candidate periods reset by the candidate period resetting unit 120.

[0045] In this way, the period detection device 1 narrows down the candidate periods based on the power spectrum, and then further subdivides the area around the narrowed-down candidate periods and adds these as candidate periods. The period of the celestial body is then detected from the narrowed-down candidate periods and the newly added candidate periods. This allows the candidate period to be detected from periods that are closer to the correct period. This makes it possible to improve the accuracy of detecting the correct period while suppressing an increase in the number of candidate periods.

[0046] Furthermore, the period detection device 1 detects, as the period of the celestial body, a candidate period with an error less than a threshold value (for example, the smallest error) from among all the candidate periods reset by the candidate period resetting unit 120. As described above, the larger the error, the further away from the correct period. Therefore, by detecting, as the period of the celestial body, a candidate period with an error less than a threshold value (for example, the smallest error) from among all the candidate periods, the correct period can be detected with high accuracy. As a result, it is possible to improve the accuracy of detecting the correct period.

[0047] (Second embodiment) In the first embodiment described above, a configuration was described that can improve the accuracy of detecting a correct period while suppressing an increase in the number of candidate periods. In contrast, the Lomb-Scargle periodogram method used in the first embodiment has room for further speedup by improving the efficiency of calculations. Therefore, in the second embodiment, a configuration will be described that can speed up the Lomb-Scargle periodogram method, thereby suppressing an increase in the number of candidate periods, and detecting a correct period quickly and with high accuracy.

[0048] As a method for speeding up the Lomb-Scargle periodogram method, in the second embodiment, the process of detecting candidate periods is speeded up by applying the double angle formula to the calculations of the Lomb-Scargle periodogram method used to narrow down candidate periods. The device configuration in the second embodiment is the same as in the first embodiment, but the process performed by the power spectrum calculation unit 110 differs from that in the first embodiment. The following description will focus on the differences from the first embodiment.

[0049] [Lomb-Scargle periodogram method] Before describing the processing performed by power spectrum calculation unit 110, differences between the Lomb-Scargle periodogram method and the first embodiment will be described. Equation (2) shown in the first embodiment is an equation for calculating the power spectrum with respect to frequency (period). In contrast, the Lomb-Scargle periodogram method in the second embodiment uses the following equation (5) to calculate the power spectrum with respect to angular frequency (period).

[0050]

number

[0051] -x and σ in equation (5) 2 are the mean and variance of the observed values, respectively, and τ is a value that satisfies the following equation (6).

[0052]

number

[0053] Next, the processing performed by the power spectrum calculation unit 110 in the second embodiment will be described. The power spectrum calculation unit 110 uses the Lomb-Scargle periodogram method for each of all given candidate periods in the input celestial object data. When performing the Lomb-Scargle periodogram method, the power spectrum calculation unit 110 uses a double angle formula, taking into account the correlation between the calculated values ​​between the candidate periods. This speeds up the Lomb-Scargle periodogram method. The speed-up of the Lomb-Scargle periodogram method will be described in detail below. Equation (5), which calculates the power spectrum corresponding to a certain ω in the Lomb-Scargle periodogram method, can be transformed into the following equation (7) using the addition theorem.

[0054]

number

[0055] A, B, C, and D in equation (7) are as shown in equation (8) below.

[0056]

number

[0057] As mentioned above, by using the addition theorem, we can find the trigonometric functions of ωτ and ωt k where t k Since is common among the candidate periods, ωt k The trigonometric functions of are related. For example, ωt k Trigonometric functions and 2ωt k The trigonometric function of holds true from the double angle formula:

[0058]

number

[0059] Equation (9) is sin(ωt k ) and cos(ωt k ) value, we can calculate sin(2ωt) using only addition, subtraction, and multiplication without using trigonometric functions. k ) and cos(2ωt k ) can be obtained. Trigonometric functions are usually calculated using methods such as Taylor series expansion. This results in high calculation costs. In contrast, by using the double angle formula described above, the number of times trigonometric functions are used can be reduced, and the calculation costs can be reduced. As a result, calculation efficiency can be improved. As a specific method for speeding up the calculation, the power spectrum calculation unit 110 uses the already calculated sinωt k and cosωt k Using the double angle formula based on sin2ωt k and cos2ωt k Furthermore, we calculate the sin2ωt k and the value of cos2ωt k By storing the value of , this speed-up method can be recursively applied to the entire candidate period.

[0060] (Processing of the power spectrum calculation unit 110) The processing flow of the power spectrum calculation unit 110 in the second embodiment will be described. Fig. 9 is an algorithm showing the procedure of the processing performed by the power spectrum calculation unit 110 in the second embodiment, and Fig. 10 is a flowchart showing the processing of the algorithm shown in Fig. 9. The processing flow of the power spectrum calculation unit 110 in the second embodiment will be described below with reference to Fig. 9. Here, a given set of candidate periods is converted into angular frequencies, and a power spectrum is calculated based on these.

[0061] The power spectrum calculation unit 110 receives as input a set P (size m) of candidate periods, an array T (size N) for time data, an array y (size N) for observed values, and the number k of periods to be extracted (step S401). Next, the power spectrum calculation unit 110 inputs, into a variable W, period p and a hash map of size m for the power spectrum of period p (step S402). The power spectrum calculation unit 110 determines whether the set P of candidate periods is empty (step S403). If the power spectrum calculation unit 110 determines that the set P of candidate periods is not empty (step S403—YES), the power spectrum calculation unit 110 selects the longest period p included in the set P (step S404). The power spectrum calculation unit 110 calculates the angular frequency ω=2π / p for the selected longest period p (step S405).

[0062] The power spectrum calculation unit 110 executes the processes from step S406 to step S409 to calculate the cosine term C and the sine term S for each angular frequency ω based on the time data (corresponding to lines 7 to 12 of the algorithm in FIG. 9). Specifically, the power spectrum calculation unit 110 updates the cosine term C with an empty array of size N (step S406), and updates the sine term S with an empty array of size N (step S407). Thereafter, the power spectrum calculation unit 110 executes loops Ls2 to Le2. In loops Ls2 to Le2, cos(ωt i ) and sin(ωt i The power spectrum calculation unit 110 calculates cos(ωt i ) is the cosine term C[i], and sin(ωt i ) is set as the sine term S[i] (step S409, step S409). The power spectrum calculation unit 110 calculates the power spectrum by the Lomb-scargle periodogram method using these terms (step S410).

[0063] Thereafter, the power spectrum calculation unit 110 executes the processes from step S411 to step S416 to update the cosine term C and the sine term S using the double angle formula (corresponding to lines 16 to 22 of the algorithm in FIG. 9) in order to perform power spectrum calculations for frequencies obtained by doubling the angular frequency ω. Specifically, the power spectrum calculation unit 110 updates the set P of candidate periods with P\{p} (step S411). The subsequent process is a process for speeding up the Lomb-scargle periodogram method using the double angle formula. The power spectrum calculation unit 110 determines whether the period p / 2 belongs to the set P (step S412). If the power spectrum calculation unit 110 determines that the period p / 2 does not belong to the set P (step S412-NO), the power spectrum calculation unit 110 executes the process of step S403.

[0064] On the other hand, if the power spectrum calculation unit 110 determines that the period p / 2 belongs to the set P (step S412—YES), the power spectrum calculation unit 110 sets the period p to period p / 2 (step S413). Thereafter, the power spectrum calculation unit 110 calculates the angular frequency ω=2π / p for the updated period p (for example, period p / 2) (step S414). The power spectrum calculation unit 110 executes loops Ls3 to Le3. In loops Ls3 to Le3, 2C[i]S[i] and 2C[i] are calculated when i is changed from 1 to N. 2 The power spectrum calculation unit 110 calculates 2C[i]S[i] when i is changed from 1 to N as the sine term S[i], and calculates 2C[i] 2−1 is set as the cosine term C[i] (steps S415 and S416). The power spectrum calculation unit 110 uses these terms to calculate the power spectrum using the Lomb-scargle periodogram method (step S417). In this way, the power spectrum calculation unit 110 calculates the power spectrum for the doubled frequency (corresponding to line 23 of the algorithm in FIG. 9). The power spectrum calculation unit 110 updates the set P of candidate periods with P\{p} (step S418). Thereafter, the power spectrum calculation unit 110 executes the process of step S403.

[0065] If the power spectrum calculation unit 110 determines in the process of step S403 that the set P of candidate periods is empty (step S403-NO), the power spectrum calculation unit 110 extracts K candidate periods in descending order of the power spectrum obtained by the above process (step S419).Then, the power spectrum calculation unit 110 outputs information indicating the extracted K candidate periods to the candidate period resetting unit 120 (step S420).

[0066] The period detection device 1 of the second embodiment configured as described above applies the double angle formula to the process of extracting candidate periods, thereby speeding up the process of detecting candidate periods. This speeds up the process of detecting candidate periods more than in the first embodiment. This makes it possible to detect correct periods quickly and with high accuracy while suppressing an increase in the number of candidate periods.

[0067] As described above, the period detection device 1 according to the second embodiment uses sinωt k and cosωt k sin2ωt calculated using the double angle formula based on k and the value of cos2ωt k By storing the value of , it is used recursively. This makes it possible to further speed up the calculation by eliminating the need for recalculation.

[0068] Although an embodiment of the present invention has been described in detail above with reference to the drawings, the specific configuration is not limited to this embodiment, and includes designs within the scope of the gist of the present invention. [Explanation of symbols]

[0069] 1... Period detection device, 10... Control unit, 110... Power spectrum calculation unit, 120... Candidate period resetting unit, 130... Period calculation unit

Claims

1. a power spectrum calculation unit that calculates the power spectrum of celestial body data, which is time-series data that represents changes in the brightness of a celestial body, for each given candidate period; a candidate period resetting unit that resets candidate periods by adding, as new candidate periods, periods surrounding each of the predetermined number of candidate periods having high power spectra, based on the power spectra of the respective candidate periods calculated by the power spectrum calculation unit; and a period calculation unit that detects the period of the celestial object from among all the candidate periods reset by the candidate period resetting unit; A period detection device comprising:

2. The candidate period resetting unit By applying the folding method to each of all the reset candidate periods, a graph showing phase information for each candidate period is created, and The period calculation unit Based on the graph created for each candidate period, a candidate period with an error less than a threshold is detected as the period of the celestial body from among all the reset candidate periods. The period detection device according to claim 1 .

3. The period calculation unit An approximate curve is drawn for the graph created for each candidate period, and the candidate period for which the mean least square error between the approximate curve and the graph is the smallest is detected as the period of the celestial body.

3. The period detection device according to claim 2.

4. The candidate period resetting unit and adding the new candidate periods by dividing the interval between the preceding and succeeding periods of each of the predetermined number of candidate periods with high power spectra by a predetermined number of divisions at equal intervals. The period detection device according to any one of claims 1 to 3.

5. The candidate period resetting unit adding a new candidate period by creating an approximation curve using a spline function based on each of the predetermined number of candidate periods with high power spectra; The period detection device according to any one of claims 1 to 3.

6. The power spectrum calculation unit Calculating a power spectrum for the celestial object data using the Lomb-scargle periodogram method, and extracting a predetermined number of candidate periods from the given candidate periods in descending order of power spectrum; The period detection device according to any one of claims 1 to 3.

7. The power spectrum calculation unit Calculate the power spectrum by using the double angle formula in the calculation of the Lomb-scargle periodogram method.

7. The period detection device according to claim 6.

8. The power spectrum calculation unit a part of the values ​​calculated using the double angle formula is stored, and a power spectrum is calculated using the part of the values ​​stored; 8. The period detection device according to claim 7.

9. Calculate the power spectrum of the celestial data, which is time-series data showing the change in brightness of the celestial object, for each given candidate period. resetting the candidate periods by adding, as new candidate periods, periods around each of the predetermined number of candidate periods with high power spectra, among the calculated power spectra for each candidate period; Detecting the period of the celestial object from all the reset candidate periods; Cycle detection method.

10. On the computer, a power spectrum calculation step of calculating the power spectrum of celestial body data, which is time series data representing changes in the brightness of a celestial body, for each given candidate period; a candidate period resetting step of resetting candidate periods by adding, as new candidate periods, periods around each of the predetermined number of candidate periods having high power spectra, based on the power spectra of the respective candidate periods calculated in the power spectrum calculation step; a period calculation step of detecting the period of the celestial object from all of the candidate periods reset in the candidate period resetting step; A program to execute.

Images