CN105067129A - Power-spectrum-based detection method for analyzing light variability period of quasi-stellar object - Google Patents

Abstract

The invention discloses a power-spectrum-based detection method for analyzing a light variability period of a quasi-stellar object, which relates to the technical field of the non-uniform sampling data periodic signal analyses. The detection method comprises steps of a signal sampling step, a data set completion step, a power spectrum calculation step, and a period detection step. To be specific, at the signal sampling step, equal-interval mean-value sampling is carried out on optical variability original data of a quasi-stellar object; at the data set completion step, data point completing based on linear interpolation for the data after sampling and data point completing based on zero interpolation for the data after sampling are carried out; at the power spectrum calculation step, the optical variability power spectrum amplitude of the quasi-stellar object is calculated; and at the period detection step, a power spectrum difference value is calculated to obtain light variability period frequency-domain signal identification information of the quasi-stellar object. The measurement precision is high; the calculation speed is fast; the analysis result is improved to the great extent; the needed period can be obtained accurately; and the stability of the light variability period of the quasi-stellar object can be disclosed.

Description

A kind of detection method based on power spectrumanalysis quasar period of light variation

Technical field:

The present invention relates to a kind of detection method period of light variation based on power spectrumanalysis quasar 3C273, belong to non-uniform sampling data periodic signal analysis technical field.

Background technology:

The method of traditional time frequency analysis is the various analysis based on Fourier transform, such as: power spectrum, small echo etc.But normal signal analytic transformation can not be suitable for the special circumstances of nonuniform sampling.Astronomical sight is subject to the restriction of a lot of objective factor, and such as weather reason, telescope equipment working condition, observed object are chosen, atmospheric seeing change etc.The disappearance of the raw data obtained, particularly accidental or random shortage of data can to analyzing to obtain the error brought of result.Disappearance degree with raw data sampling promotes, and in analysis result, pseudoperiod becomes the accuracy of branch's impact analysis result.Power spectrum analysis method is one of the basic cycle analytical approach based on Fourier transform, and the data processing relating to time-frequency domain conversion is drawn in this way mostly.Although power spectrum analysis method has unique advantage in process stationary signal, it just seems and scarce capacity is certain to occur deviation when practical application in process nonuniform sampling irregular data.The method of a kind of cycle detection based on power spectrumanalysis quasar of the present invention adopts special power spectrum to mend " 0 " to sampled data linear interpolation polishing and sampled data to missing information composition and asks poor method, solve the cycle frequency-region signal computational analysis of nonuniform sampling preferably.

Summary of the invention:

For the problems referred to above, the technical problem to be solved in the present invention is to provide a kind of detection method based on power spectrumanalysis quasar period of light variation.

A kind of detection method based on power spectrumanalysis quasar period of light variation of the present invention, its detecting step is signal sampling step; The complete step of data group; Spectra calculation step; Cycle detects step; Described signal sampling step is that quasar light change raw data carries out average sampling at equal intervals; The complete step of described data group is data " 0 " interpolation polishing data point after will sample rear data linear interpolation polishing data point and sampling; Described spectra calculation step compute classes celestial body light Variable power spectral amplitude ratio; The described cycle detects step and obtains quasar frequency-region signal period of light variation identification information for asking power spectrum difference.

As preferably, described data component does not obtain cycle frequency-region signal information by linear interpolation, " 0 " polishing and power spectrum difference.

Beneficial effect of the present invention is: its measuring accuracy is high, and computing velocity is fast, improves the result of analysis to a great extent, not only can obtain the required cycle more accurately, and can disclose the stability of quasar period of light variation.

Accompanying drawing illustrates:

For ease of illustrating, the present invention is described in detail by following concrete enforcement and accompanying drawing.

Fig. 1 is principle of the invention figure;

Fig. 2 is the analytical control of power spectrum method described in the present invention signal, and no signal point data is eliminated;

Fig. 3 is the analytical control of power spectrum method described in the present invention signal, and the signaling point of about 25% is eliminated at random;

Fig. 4 is the analytical control of power spectrum method described in the present invention signal, and the signaling point of about 50% is eliminated at random;

Fig. 5 is the analytical control of power spectrum method described in the present invention signal, and the signaling point of 25% is eliminated at equal intervals;

Fig. 6 is the analytical control of power spectrum method described in the present invention signal, and the signaling point of 50% is eliminated at equal intervals;

Fig. 7 is B wave band flow and the spectrum index light curve of quasar 3C273 in the present invention;

Fig. 8 is B wave band flow and the spectrum index light curve of quasar 3C273 in the present invention, by described power spectrum method and Jurkevich methods analyst result.

Embodiment:

For making the object, technical solutions and advantages of the present invention clearly understand, below by the specific embodiment shown in accompanying drawing, the present invention is described.But should be appreciated that, these describe just exemplary, and do not really want to limit the scope of the invention.In addition, in the following description, the description to known features and technology is eliminated, to avoid unnecessarily obscuring concept of the present invention.

As shown in Figure 1, this embodiment by the following technical solutions: its detecting step is signal sampling step; The complete step of data group; Spectra calculation step; Cycle detects step; Described signal sampling step is that quasar light change raw data carries out average sampling at equal intervals; The complete step of described data group is data " 0 " interpolation polishing data point after will sample rear data linear interpolation polishing data point and sampling; Described spectra calculation step compute classes celestial body light Variable power spectral amplitude ratio; The described cycle detects step and obtains quasar frequency-region signal period of light variation identification information for asking power spectrum difference.

Further, described data component does not obtain cycle frequency-region signal information by linear interpolation, " 0 " polishing and power spectrum difference.

The Method And Principle of this embodiment: in random sampled data, there is inevitably long or short data slit in time series, during power spectrum method process data, requires that gap is very short, best without sampled data disappearance.For power spectrum analysis method, as long as artificial supplying missing data, adopt which kind of approximating method all likely to introduce error, cause the out of true of cycle analysis; And know for the possibility cycle debating the pseudoperiod that error is formed.

The application of this embodiment:

If there is signal function f (t) to meet Di Liheli condition, then its Fourier transform is:

$F\left(\mathrm{\ω}\right)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}f\left(t\right)e^{-i2\mathrm{\π}\mathrm{\ω}t}dt$

If fully little at limited region [-T, T] outer f (t) fourier integral, then what can be similar to thinks:

$F\left(\mathrm{\ω}\right)={\∫}_{-T}^{T}f\left(t\right)e^{-i2\mathrm{\π}\mathrm{\ω}t}dt$

[-T, T] is carried out 2N decile, gets [-T, T] length of an interval degree and Along ent is:

$\begin{array}{cc}\mathrm{\Δ}t=\frac{2T}{2N}=\frac{T}{N}& t_k=k\mathrm{\Δ}t,k=-N,-N+1,...,N-1,N\end{array}$

Then

$F\left(\mathrm{\ω}\right)={\∫}_{-N\mathrm{\Δ}t}^{N\mathrm{\Δ}t}f\left(t\right)e^{-i2\mathrm{\π}\mathrm{\ω}t}dt=\underset{k=-N}{\overset{N-1}{\Σ}}\underset{k\mathrm{\Δ}t}{\overset{(k+1)\mathrm{\Δ}t}{\∫}}f\left(t\right)e^{-i2\mathrm{\π}\mathrm{\ω}t}dt$

By ω discretize, get

$\mathrm{\ω}=\frac{1}{N\mathrm{\Δ}t},{\mathrm{\ω}}_j=j\mathrm{\ω},j=0,1,...,N-1$

Have

$\begin{array}{c}F\left({\mathrm{\ω}}_j\right)\≈\underset{k=-N}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_k\right)e^{-\frac{i2\mathrm{\π}jk}{N}}\mathrm{\Δ}t\\ =\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_{k-N}\right)e^{-\frac{i2\mathrm{\π}j\left(k-N\right)}{N}}\mathrm{\Δ}t+\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_k\right)e^{-\frac{i2\mathrm{\π}jk}{N}}\mathrm{\Δ}t\\ =\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}\[f\left(t_{k-N}\right)+f\left(t_k\right)\]e^{-\frac{i2\mathrm{\π}jk}{N}}\mathrm{\Δ}t,j=0,1,\mathrm{...},N-1\end{array}$

Note:

A _k＝[f(t _k-N)+f(t _k)]Δt,k＝0,1,…,N-1

Obtain

$\begin{array}{c}{x}_j=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{A}_k{e}^{-\frac{i2\mathrm{\π}jk}{N}}=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{A}_k{W}_N^{-jk}\\ j=0,1,\mathrm{...},N-1,i=\sqrt{-1},W_N=e^{\frac{2\mathrm{\π}i}{N}}\end{array}$

Note:

$A_k=f\left(t_k\right)g\left(\frac{t_k}{T}\right),k=0,1,...,N-1,$ Then

$F\left({\mathrm{\ω}}_j\right)=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{A}_k{W}_N^{-jk}=\underset{k=0}{\overset{N-1}{\Σ}}f\left(t_k\right)g\left(\frac{t_k}{T}\right)W_N^{-jk}=\frac{1}{\mathrm{\Δ}t}\underset{0}{\overset{T}{\∫}}f\left(t\right)g\left(\frac{t}{T}\right)e^{-i2\mathrm{\π}j\mathrm{\ω}t}dt,j=0,1,...,N-1$

Then seasonal effect in time series power Spectral Estimation is:

$P_j=P\left(j\mathrm{\ω}\right)=\frac{1}{T}{\left|\underset{0}{\overset{T}{\∫}}f\left(t\right)g\left(\frac{t}{T}\right)e^{-i2\mathrm{\π}j\mathrm{\ω}t}dt\right|}^2=\frac{\mathrm{\Δ}t}{N}{\left(x_j\right)}^2,j=0,1,...,N-1$

ω ₀it is the natural frequency of signal f (t).When n is large, P _j=P (j ω) is at ω ₀the neighbouring peak forming a projection, and less at the density value of other frequency place power spectrum.According to this feature, we just can find the periodic component in signal.Equally when surely containing multiple periodic component in signal, P _jthere will be multiple maximum value in=P (j ω), need the cycle of distinguishing the true from the false.

For power spectrum method, select different window functions effectively can suppress spectral aliasing and spectral line leakage phenomenon.The window function of usual selection has BoxcarWindow, TriangleWindow, HammingWindow, HanningWindow and WelchWindow etc.

BartlettWindow

$G\left(t\right)=\frac{1}{2}(1-cos\frac{2\mathrm{\π}j}{N})$ HanningWindow

$G\left(t\right)=1-{\left(\frac{j-\frac{1}{2}N}{\frac{1}{2}N}\right)}^2$ WelchWindow

If suppose that sampling function is f (t), then the function of missing at random data can be defined as:

Actual discrete data sequence then after sampling can be defined as the function of f (t) δ (t), and wherein t is the sequential variable of periodic signal.Reference power spectral method derivation, if function f (t) meets Di Liheli condition, then its Fourier transform is:

$F\left(\mathrm{\ω}\right)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}f\left(t\right)\mathrm{\δ}\left(t\right)e^{-i2\mathrm{\π}\mathrm{\ω}t}dt$

If fully little at limited region [-T, T] outer f (t) fourier integral, then what can be similar to thinks:

$F\left(\mathrm{\ω}\right)={\∫}_{-T}^{T}f\left(t\right)\mathrm{\δ}\left(t\right)e^{-i2\mathrm{\π}\mathrm{\ω}t}dt$

[-T, T] is carried out 2N decile, gets [-T, T] length of an interval degree and Along ent is:

$\begin{array}{cc}\mathrm{\Δ}t=\frac{2T}{2N}=\frac{T}{N}& t_k=k\mathrm{\Δ}t,k=-N,-N+1,...,N-1,N\end{array}$

Then

$F\left(\mathrm{\ω}\right)=\underset{k=-N}{\overset{N-1}{\Σ}}\underset{k\mathrm{\Δ}t}{\overset{(k+1)\mathrm{\Δ}t}{\∫}}f\left(t\right)\mathrm{\δ}\left(t\right)e^{-i2\mathrm{\π}\mathrm{\ω}t}dt$

By ω discretize, get

$\mathrm{\ω}=\frac{1}{N\mathrm{\Δ}t},{\mathrm{\ω}}_j=j\mathrm{\ω},j=0,1,...,N-1$

Have

$\begin{array}{c}F\left({\mathrm{\ω}}_j\right)\≈\underset{k=-N}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_k\right)\mathrm{\δ}\left(t_k\right)e^{-\frac{i2\mathrm{\π}jk}{N}}\mathrm{\Δ}t\\ =\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}\[f\left(t_{k-N}\right)\mathrm{\δ}\left(t_{k-N}\right)+f\left(t_k\right)\mathrm{\δ}\left(t_k\right)\]e^{-\frac{i2\mathrm{\π}jk}{N}}\mathrm{\Δ}t,j=0,1,\mathrm{...},N-1\end{array}$

Note: A _k=[f (t _k-N) δ (t _k-N)+f (t _k) δ (t _k)] Δ t, k=0,1 ..., N-1

From the derivation of power spectrum method, owing to introducing window function to process spectral aliasing and spectral line leakage phenomenon, if therefore sampling is continuous signal, or the discrete sampling of continuous signal, then should there be violent spectral aliasing and spectral line leakage phenomenon in result, but, can notice to be exactly because A here _kbe defined as [f (t _k-N)+f (t _k)] Δ t, if f is (t _k-N) and f (t _k) not desirable periodic signal, then A _kvalue will depart from far away power spectrum method derive in A _k.

Solve A _k=[f (t _k-N) δ (t _k-N)+f (t _k) δ (t _k)] δ (t) function in Δ t, notice the Fourier transform for power spectrum method, have the principle of conservation of energy in Du Hameier integration and Fourier transform, note:

$F\left({\mathrm{\ω}}_j\right)=\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{x}_j\mathrm{\δ}(t-t_k)+\underset{0}{\overset{t}{\∫}}\frac{df\left(t_k\right)}{{\mathrm{dt}}_k}\mathrm{\δ}(t-t_k){\mathrm{dt}}_k$

Wherein $x_j=\underset{K=0}{\overset{N-1}{\Σ}}f\left(t_k\right)g\left(\frac{t_k}{T}\right)W_N^{-jk},$ And $\underset{0}{\overset{t}{\∫}}\frac{df\left(t_k\right)}{{\mathrm{dt}}_k}\mathrm{\δ}(t-t_k){\mathrm{dt}}_k=\underset{0}{\overset{t}{\∫}}\mathrm{\δ}(t-t_k)df\left(t_k\right)$

If order $\stackrel{\‾}{f\left(t\right)}=\frac{1}{N}\underset{t}{\overset{t_k}{\Σ}}f\left(t_k\right)\mathrm{\δ}\left(t_k\right)$

Obtain

$A_k=\[f\left(t_k\right)\mathrm{\δ}\left(t_{k-N}\right)+\stackrel{\‾}{f\left(t\right)}\]\mathrm{\Δ}t$

Then this method seasonal effect in time series power Spectral Estimation is:

$P_j=P\left(j\mathrm{\ω}\right)\≈\frac{\mathrm{\Δ}t}{N}\left|{\left(\underset{K=0}{\overset{N-1}{\mathrm{\Σ}}}\stackrel{\‾}{f\left(t\right)}g\left(\frac{t_k}{T}\right)W_N^{-jk}\right)}^2-{\left(\underset{K=0}{\overset{N-1}{\mathrm{\Σ}}}f\left(t_k\right)g\left(\frac{t_k}{T}\right)\mathrm{\δ}\left(t_k\right)W_N^{-jk}\right)}^2\right|$

When applying described power spectrum method, it should be noted that sampling rate must be greater than the minimum sample rate required by Nyquist Sampling Theorem equally.

Fig. 1 be equal interval sampling power spectrum analysis calculate (filtering) process schematic diagram, wherein circle represents the step to former data manipulation, the step of box indicating analytical calculation.In the diagram, we to f (t) and employ identical window function.

In this method, former history light curve does not significantly introduce function or the data point of hypothesis.And after comparing with primitive period collection of illustrative plates method of estimation, the computation process of whole cycle collection of illustrative plates, except adding window function to process except spectral aliasing and spectral line leakage, does not carry out other interference; And for the interval produced due to historical data data, eliminate again after power Spectral Estimation calculates, the analysis therefore for data information has certain confidence level.

In order to check the reliability of this power spectrum method, we check with the periodic signal of a group.Actual result as astronomical sight all includes error, we used two sin wave functions and adds white Gaussian noise to simulate the long period light parameter certificate of celestial body.Data sequence employs 13000 data points.

$x\left(t\right)=14+\underset{i=1}{\overset{2}{\Σ}}{A}_{i}sin2{\mathrm{\πf}}_{i}t+\mathrm{\ϵ}$

Wherein, f ₁and f ₂natural period be T ₁a cycle is repeated 65 data points; T ₂a cycle is repeated 65 data points.The amplitude of two Sin wave functions is A ₁=4.3 magnitudes; A ₂=6.5 magnitudes, and average white Gaussian noise is 0.01 magnitude (random variation, unit is magnitude).In order to the impact of matching observed result and shortage of data, we first stochastic censored except partial test signal, then calculate its cycle collection of illustrative plates, two the sin wave function abilities supposing that light becomes can be judged according to cycle collection of illustrative plates, judge the reliability of this power spectrum method.Fig. 2 to Fig. 6 is the result of inspection equal interval sampling power spectrum analysis method reliability, Fig. 2 is the analysis result of periodic inspection signal when not deleting any data point, obviously can see in Fig. 2 having two obvious peak values from result, the frequency of these two peak values correspond to two simulation cycles in inspection signal.Fig. 3 and Fig. 4 is respectively periodic inspection signal in the result of stochastic censored except the equal interval sampling power spectrum analysis in about 25% and 50% situation, can see that from Fig. 3 and Fig. 4 the crest frequency that two simulation cycles are corresponding does not obviously change, even if when having illustrated for random loss in astronomical observation below 50% data point, described analytical approach still effectively.Fig. 5 and Fig. 6 is another kind of situation, namely 25% and 50% periodic inspection signal is deleted at equal intervals, carry out the shortage of data of the time period at equal intervals in simulated light varied curve, Fig. 5 and Fig. 6 is respectively periodic inspection signal in the result deleting the sampling period atlas analysis in 25% and 50% situation at equal intervals.Can see from Fig. 5 and Fig. 6, equally when having lacked below 50% data point, described analytical approach still effectively.

The practical application of described invention: radio source 3C273 is first celestial body being recognized as quasar.3C273 is the quasar of the brightest (magnitude of average B wave band is 12.83) and nearest (red shift is 0.158), the multiband observation of 3C273 has been had to the history of 40 years, therefore 3C273 is one of celestial body that research is the most thorough from radio wavelength band to gamma-rays wave band.In history, this celestial body shows significantly light and becomes on different wave bands, and a lot of scientist shows keen interest to the potential cycle finding this celestial body.Particularly at optical region, Smith has searched out the possible cycle of 13 years; Use Fourier analysis method, the result of Kunkel shows the faint light variable period of about 17.9.After employing Jurkevich method, Jurkevich have found the possible cycle of 15 years.Searched out the possible cycle of 9 years and 16 years at Ozernoi in 1997, and the cycle of 16 years by the VLBI of radio observe project confirm.It is extremely important that the theoretical model of research quasar radiation mechanism analyzes light curve, and therefore for the history light curve research of 3C273, cycle analysis is very meaningful.

In order to study the power spectrum change of quasar 3C273, have collected celestial body in 35 years up to now optics B, V, R wave band can data to build optical narrow-band spectrum change light curve, and life cycle atlas analysis method calculate and use Jurkevich method to confirm spectrum change periodicity.The magnitude value of 3C273 is scaled data on flows, with logF _bform represent B wave band light-metering flow.In the history light curve of B wave band, the peak value demonstrating outburst appears at 1971,1978,1982,1992 and 2002, and calculates optical region spectrum index.Fig. 7 is B wave band flow and the spectrum index light curve of quasar 3C273.Light curve data point in figure comes from the paper of 1968,1969,1970,1975,1980 of Burkhead, the paper of 1998 of Courvoisier, the paper of 1989 of Impey, the paper of 2008 of Soldi, the paper of 1999 of T ü rler and network data base Http: //isdc.unige.ch/3c373.

The result of calculation that have employed described power spectrum method and Jurkevich method provides in fig. 8.In fig. 8, the consistent cycle is obtained from described power spectrum method and Jurkevich method, namely 15.7 and 10.5, these cycles are consistent in B wave band history light curve with spectrum index, and the cycle of 15.7 and Jurkevich analyze B wave band history light curve and the long-term light curve of spectrum index are consistent, in addition, these two results with in history to the discussion result of the possible period of light variation of 3C273 also closely.

Shown by described power spectrum research, the narrow-band spectrum change of the optical region of quasar 3C273 has periodically, and its cycle is in a disguised form consistent with the multiband optical-flow light of celestial body.In the history light curve of optical region, there is not the situation of obvious time delay, the radiation therefore between different-waveband because of this in same relatively large region out, can not the detailed information in resolving radiation region in history light curve.Because spectrum change has periodically, and identical with the light curve cycle, and therefore the radiation of different-waveband has identical origin, and a nearly step describes the Binary Black Hole model that may exist.

More than show and describe ultimate principle of the present invention and principal character and advantage of the present invention.The technician of the industry should understand; the present invention is not restricted to the described embodiments; what describe in above-described embodiment and instructions just illustrates principle of the present invention; without departing from the spirit and scope of the present invention; the present invention also has various changes and modifications, and these changes and improvements all fall in the claimed scope of the invention.Application claims protection domain is defined by appending claims and equivalent thereof.

Claims (2)

1. based on the detection method of power spectrumanalysis quasar period of light variation, it is characterized in that: its detecting step is signal sampling step; The complete step of data group; Spectra calculation step; Cycle detects step; Described signal sampling step is that quasar light change raw data carries out average sampling at equal intervals; The complete step of described data group is data " 0 " interpolation polishing data point after will sample rear data linear interpolation polishing data point and sampling; Described spectra calculation step compute classes celestial body light Variable power spectral amplitude ratio; The described cycle detects step and obtains quasar frequency-region signal period of light variation identification information for asking power spectrum difference.

2. a kind of detection method based on power spectrumanalysis quasar period of light variation according to claim 1, is characterized in that: described data component does not obtain cycle frequency-region signal information by linear interpolation, " 0 " polishing and power spectrum difference.