Stability of cycle in samplogram and spurious cycles in solar activity
aa r X i v : . [ a s t r o - ph . S R ] S e p Stability of cycle in samplogram and spurious cycles in solar activity
Kim Chol-jun
Faculty of Physics, Kim Il Sung University, DPR Korea email address: [email protected]
September 22, 2020
Abstract
The spectral analysis with stochastic significance cannot distinguish the spurious cycles effectively, be-cause a non-stationary signal can make a significant peak in spectrum. I show that the random separationbetween the grand extremes such as grand maxima and minima in solar activity can make a spurious butsignificant peak in spectrum. And it is possible to pick even the weak stable cycle by applying an averagingdown-sampling and the differencing to a random signal. This is because both operations variously changethe height of peak in spectrum and the spurious cycle is in turn unstable in both operations. I introduce asamplogram showing these operations and stability intuitively.
The cycles in solar activity have long been discussed. In mid-19th century Schwabe and Wolf had discoveredthe well-known 11-yr solar cycle. Later, the Gleissberg cycle ( ∼
80 yr) was proposed, which, however, has dimly(60 ∼
120 yr) identified in the spectrum.The reconstruction of the past solar activity makes it possible to study the long-term solar cycles. The recon-struction from cosmogenic radionuclides such as C and Be (Beer, Tobias & Weiss, 2018; Chol-jun & Jik-su,2020; Eddy, 1976; Solanki et al., 2004; Steinhilber et al., 2012; Usoskin, 2008, 2016; Wu et al., 2018) have im-plied the various solar cycles: The Gleissberg (88-yr), Suess/de Vries (208-yr), Hallstatt ( ∼ ∼ ∼ ∼
500 and ∼ In this section we see the variation of peak for the sinusoidal signal when the sampling interval changes.
We start with a multi-cyclic sinusoidal signal: u ( t ) = X j U j cos (cid:18) πtT j + φ j (cid:19) , (1)where U j , T j and φ j are the amplitude, period and phase of the j -th mode of the signal, respectively. A timeseries sampled from this signal with a time-step ∆ t is expressed as u ( r ) = X j U j cos (cid:18) πr ∆ tT j + φ j (cid:19) = 12 X j U j (cid:20) exp( iφ j ) exp (cid:18) πir ∆ tT j (cid:19) + exp( − iφ j ) exp (cid:18) − πir ∆ tT j (cid:19)(cid:21) , (2)where r is the index of point in the time series and runs from 1 to N that is the length of the time series. Also, t ≈ r ∆ t .Apply a down-sampling to the time series, which means the sampling interval getting greater or a coarsersampling. The down sampling can have two approaches: hopping and averaging. A hopping down-samplingtakes points of the series after every several steps (e.g. m steps) to make a new series, which causes the loss of2nformation because some data points are missed in the down-sampled series. However, the data points of thedown-sampled series are already included in the original series. When the sampling interval increasing in thisdown-sampling, the total number of points is reduced by the same factor as the sampling interval increases.This approach is used in signal acceptance of machines with cares of low acceptance frequency. The hoppingdown-sampled series is expressed as u m,k ( r ) = 12 X j U j (cid:20) exp( iφ j ) exp (cid:18) πi ∆ t ( k + m ( r − T j (cid:19) + exp( − iφ j ) exp (cid:18) − πi ∆ t ( k + m ( r − T j (cid:19) (cid:21) , (3)where m is the sampling step and k is an integer between 1 and m , that is, we can select one series among m down-sampled series. r = 1 ∼ n is the index of point in the down-sampled series. The length of thedown-sampled series is approximated to be n = IntegerPart ( N/m ) , (4)that is, if a perfect factorization N = mn is impossible, then integer part of N/m is taken as the length of thedown-sampled series, which means that in down-sampling the series is truncated. However, this approach willsimplify analysis.On the other hand, an averaging down-sampling is to average every m points and take it as a point of thenew down-sampled series. The averaging down-sampled series is expressed as u ¯ m ( r ) = 1 m m X k =1 u m,k ( r ) (5)= 12 m X j U j (cid:20) exp( iφ j ) exp (cid:18) πi ( r − m ∆ tT j (cid:19) exp (cid:18) πi ∆ tT j (cid:19) − exp (cid:16) πim ∆ t T j (cid:17) − exp (cid:16) πi ∆ t T j (cid:17) + exp( − iφ j ) exp (cid:18) − πi ( r − m ∆ tT j (cid:19) exp (cid:18) − πi ∆ tT j (cid:19) − exp (cid:16) − πim ∆ t T j (cid:17) − exp (cid:16) − πi ∆ t T j (cid:17) (cid:21) , which can be derived from Eq. (3) by using formula for the geometric series: m X k =1 exp (cid:18) πik ∆ tT (cid:19) = exp (cid:18) πi ∆ tT (cid:19) − exp (cid:16) πim ∆ t T (cid:17) − exp (cid:16) πi ∆ t T (cid:17) . (6)In the hopping down-sampling the data points are inherited intact to the down-sampled series. This gives arisk to inherit the spurious cycle. However, in the averaging down-sampling, the new data points are forged byaveraging and, in that way, the noise is suppressed and the risk that the spurious cycle is inherited intact to thedown-sampled series lowers. From now on, in this paper, we will mainly consider the averaging down-sampling.The power spectrum is obtained by the Fourier transform of autocorrelation of the signal. If two sampledsignals x ( r ) and y ( r ) are given, the correlation between them is defined as a normalized covariance: Cor ( x, y ) = Cov ( x, y ) p Var ( x ) Var ( y ) . (7) Cov ( x, y ) stands for the covariance between x ( r ) and y ( r ) of length N : Cov ( x, y ) = 1 N − N X r =1 ( x ( r ) − ¯ x )( y ( r ) − ¯ y ) ∗ , (8)where the asterisk means complex conjugate. If we deal only with real signals, the conjugate symbol disappears.¯ x and ¯ y stand for the mean values of x ( r ) and y ( r ): ¯ x = N P Nr =1 x ( r ) and ¯ y = N P Nr =1 y ( r ). The variances of x ( r ) and y ( r ) are defined in a similar way as the covariance: Var ( x ) = 1 N − N X r =1 ( x ( r ) − ¯ x )( x ( r ) − ¯ x ) ∗ , Var ( y ) = 1 N − N X r =1 ( y ( r ) − ¯ y )( y ( r ) − ¯ y ) ∗ . (9)3he autocorrelation of a sampled signal x ( r ) is defined in the same way as the correlation (Eq. 7). AC ( x, g ) = Cov ( x ( r ) , x ( r + g )) p Var ( x ( r )) Var ( x ( r + g )) , (10)where g is shifting steps (or lag), that is, the autocorrelation of signal can be the correlation between a subseries x ( r ) starting at the first point of the series and some (e.g. g ) steps shifted (or lagged) subseries x ( r + g ) endingat the last point of the series. If the length of the time series is fixed, the lengths of both shifted series must betruncated by the shifting steps.To evaluate the power spectrum or autocorrelation (AC) periodogram we use the discrete Fourier transform(DFT). The DFT v ( s ) of a time series u ( r ) with time-step ∆ t and length N is defined as v ( s ) = 1 √ N N X r =1 u ( r ) exp (cid:18) πi ( r − s − N (cid:19) , (11)where s is the index of point of the series v ( s ) and stands for the frequency index. The frequency and periodcorresponding to the s -th point of the series v ( s ) are evaluated as f = s ∆ f and T = N ∆ ts − . (12) s runs from 2 to N and the frequency step ∆ f has relation with the time step ∆ t as ∆ f ∆ t = N . For simplicity,take following replacements: s − q, (13)∆ tT j = q j N ≈ q j mn . (14)The frequency index s is replaced with q = s − T j of the j -th mode is replaced with frequency q j . Remember that in the averaging down-sampling the whole length of the series N is truncated into areducible product mn where m is averaging step and n = IntegerPart ( N/m ) is the length of the down-sampledseries(Eq. 4).Now, it is possible to evaluate the power or AC periodogram of the down-sampled time series u ¯ m ( r ) of length n by Eq. 11: P Sm ( s ) = 1 √ n n X g =1 ACm ( g ) exp (cid:18) πi ( g − s − n (cid:19) , (15)where ACm ( g ) is the autocorrelation series of the averaging down-sampled time series u ¯ m ( r ) and the shiftingstep g is used as a time variable for the Fourier transformation. Do not confuse g with the frequency index q .Taking care of Eq. (13), a long algebraic calculation with numerical comparisons gives us the expression for ACperiodgram at the leading order P Sm ( q ) = 1 √ n n X g =1 ACm ( g ) exp (cid:18) πi ( g − qn (cid:19) ≈ √ n n X g =1 P j U j A j P k U k A k (cid:20) exp (cid:18) πigq j n (cid:19) + exp (cid:18) − πigq j n (cid:19)(cid:21) exp (cid:18) πi ( g − qn (cid:19) = 12 √ n P k U k A k X j U j A j exp (cid:18) πiq j n (cid:19) − exp (2 πi ( q + q j ))1 − exp (cid:16) πi ( q + q j ) n (cid:17) + exp (cid:18) − πiq j n (cid:19) − exp (2 πi ( q − q j ))1 − exp (cid:16) πi ( q − q j ) n (cid:17) , (16)where A j = 1 − exp (cid:18) πiq j n (cid:19) − exp (cid:18) πiq j mn (cid:19) − exp (cid:18) − πiq j n (cid:19) − exp (cid:18) − πiq j mn (cid:19) . (17)4
00 400 600 800 1000Period1234Power (a) (b)
Figure 1:
The power spectrum and variation of the power in averaging down-sampling of a composite signal of U =1 , U = 2 , U = 3 , T = 95 , T = 210 , T = 330 , φ = 1 . , φ = 1 , φ = 2 . , ∆ t = 10 and N = 200. (a) The powerspectrum at sampling step m = 1. (b) Variation of power of the j = 3-rd mode when increasing the sampling step m . The solid line shows the numerical evaluation and dashed line does the approximate one. The dotted line in (b)means the degradation proportional to √ m . The red vertical bar means the Nyquist sampling step for the 3-rd mode: m = T t = 16 .
5. At the Nyquist sampling the sampling peak appears.
Using Eq. (16), we can simulate the power spectrum of the signal and variation of the power in the averagingdown-sampling. Fig.1 shows the numerical and approximate evaluations of power spectrum and the variation ofthe power in the averaging down-sampling (the numerical one is evaluated by the default functions in language
Mathematica and the approximate one is by Eq. 16). The numerical and approximate evaluations show a goodagreement.In Eq. (16), the quantity in square brackets reflects the frequency dependency and gives some “peaks”:
P Smcore ( q ) = exp (cid:18) πiq j n (cid:19) − exp (2 πi ( q + q j ))1 − exp (cid:16) πi ( q + q j ) n (cid:17) + exp (cid:18) − πiq j n (cid:19) − exp (2 πi ( q − q j ))1 − exp (cid:16) πi ( q − q j ) n (cid:17) , (18)where there appear two ”singularities” at q = q j and q + q j = n . In fact, they are not the real singularities.In approaching q → q j , we find that lim q → q j − exp (2 πi ( q − q j ))1 − exp (cid:16) πi ( q − q j ) n (cid:17) = n. (19)This corresponds to the maximum of the left hand side. So the both peaks appear: the former “singularity” inEq. (18) makes a true peak at frequency q = q j in spectrum and the latter does an aliasing peak at frequency q = n − q j . The frequency domain of the Fourier spectrum is given in 0 ≤ f ≤ f s , where f s = t is a samplingfrequency. A region of 0 ≤ f ≤ f s is a true frequency region and a region of f s ≤ f ≤ f s is an aliasing frequencyregion. According to the sampling theorem, only modes of frequency f < f s , i.e. frequency less than theNyquist sampling frequency f N = f s , are effectively evaluated. The aliasing peak at q = n − q j in the aliasingregion f s ≤ f ≤ f s is no more than a reflection of the true peak at q = q j in the true region, which is obviousfrom a symmetry of the above expressions. Ordinarily, the aliasing region (and aliasing peak) is neglected inthe power or Fourier spectra. Actually, the frequency q j at which the true peak occurs does not depend on thesampling interval m ∆ t or the length n of the down-sampled series. This says that when increasing the samplinginterval m the true peak keeps at the same frequency q = q j . This is just a down-sampling (or sampling)stability of cycle . Such endured peaks when increasing sampling interval form a track of peaks like a ridge,which we will see later.In fact, in the DFT, the true peak should happen not exactly at q = q j but at q = Round ( q j j ) and thealiasing peak at q = Round ( n − q j ) because q is not always integer, where Round () means the rounding value.If q = n , the true peak around q = q j and the aliasing peak around q = n − q j are overlapped, which makesanother type of peak. We call this peak a sampling peak of the j-th mode because the overlapping frequency q = n corresponds to the Nyquist sampling m ∆ t = T (Eq. 14). In comparison with the sampling peak, thetraditional true peak at m = 1 could be called a resonance peak , which represent the cycle in the Fourier andpower spectra. In fact, if we consider a case of the Nyquist sampling of the j -th mode, it holds m Nj ∆ t = T j n Nj = 2 q j , (20)5
10 15 20 m0.51.01.5Power (a) (b)
Figure 2:
The variations of the power of the j = 2-nd mode in Fig. 1. The length of the Nyquist down-sampled seriesis (a) n = 20 and (b) n = 19. The solid line stands for the numerical evaluation, the dashed line for the approximateevaluation according to Eq. (16) and the dotted line for the degradation proportional to √ m . In (a) the approximateevaluation (dashed line) shows a valley at the Nyquist sampling. In (b) the valley disappears. The red vertical bar meansthe Nyquist sampling m = 10 .
5. The blue dotted vertical bar stands for the nodal peak with the j =1-st mode and theorange dashed vertical bar for the nodal peak with the j =3-rd mode (see § from which P Sm ( q j , n Nj ) ≈ √ n Nj P k U k A k U j A j (cid:20) exp (cid:18) πiq j n Nj (cid:19) + exp (cid:18) − πiq j n Nj (cid:19)(cid:21) ≈ − √ n Nj P k U k A k U j A j , (21)where m Nj and n Nj are the averaging step and the length of the down-sampled time series u ¯ m ( r ) in Eq. 5.Thus, it is obvious that both terms in square brackets in Eq. (21) have the same phase so they seem not cancelout each other and form the peak at the Nyquist sampling.However, instead of the sampling peak it is possible to appear a valley. An example is shown in numericalsimulation of the power degradation in the down-sampling for the j = 2-nd mode in Fig. 1. A valley appears atthe Nyquist sampling! (Fig. 2(a)) In this case, the power at the Nyquist sampling is evaluated in Eq. (16) with m = 10 , n = 20 , q = 10 and q = q j = 9 . T = 210. Here we can find 2 q = n . It seems that thereal value of power at 2 q = n should be related to the appearance of valley through a possibility that the powercan become zero. Investigate a value range of the power in Eq. (18) for the above case. The value of Eq. (18)approaches to 0 when q j comes to 9.5238 and 10.4762. In DFT a resonance peak happens at q = Round ( q j )instead of q j , so in the case of q − q j ≈ . n is odd. For example, we can set n = 2 IntegerPart ( n ) −
1. Then Eq. (18)cannot to be real and, simultaneously, we can avoid 0, i.e. the appearance of valley. However, for the odd n ,the peak at q j = q becomes lower. Actually, the maximum power lowers by ∼ n , which means there disappear a significant sampling peak over the ridge. But the peak becomes broader infrequency domain for the odd n than for the even n . In fact, zero points of q j are 8.0 and 11.0, so the peakbroadens 3 times. Thus, we see that the valley at the Nyquist sampling can be avoided by that the lengthof every down-sampled series should be odd. However, then the peak as well as the valley might disappear(Fig. 2b).Instead of the normalized covariance (Eq. 10), we can define the autocorrelation as the absolute covariance AC ( x, g ) = Cov ( x ( r ) , x ( r + g )) , (22)which is popular in many books (cf. Oppenheim, Schafer & Buck (1999)). Then the power of the averagingdown-sampled time series u ¯ m ( r ) in Eq. 5 can be evaluated as P Sm ( q ) ≈ m √ n X j U j A j exp (cid:18) πiq j n (cid:19) − exp (2 πi ( q + q j ))1 − exp (cid:16) πi ( q + q j ) n (cid:17) + exp (cid:18) − πiq j n (cid:19) − exp (2 πi ( q − q j ))1 − exp (cid:16) πi ( q − q j ) n (cid:17) . (23)We can also induce the Fourier transform of the averaging down-sampled time series u ¯ m ( r ) more easily: v ¯ m ( q ) ≈ X j U j m √ n (cid:26) exp( iφ j ) exp (cid:18) πimn q j (cid:19) − exp (cid:0) πin q j (cid:1) − exp (cid:0) πimn q j (cid:1) − exp [2 πi ( q + q j )]1 − exp (cid:2) πin ( q + q j ) (cid:3) + exp( − iφ j ) exp (cid:18) − πimn q j (cid:19) − exp (cid:0) − πin q j (cid:1) − exp (cid:0) − πimn q j (cid:1) − exp [2 πi ( q − q j )]1 − exp (cid:2) πin ( q − q j ) (cid:3) (cid:27) . (24)6t is obvious that the Fourier amplitude and (normalized and absolute) powers have the sampling peak at q = q j = n Nj . Furthermore, they are real. So their complex argument could be 0 or ± π and the power ofsampling peak should be negative real. However, the DFT violates such an expectation. The Fourier transformis easier than the power to extract the information of cycle such as the amplitude and phase. On the otherhand, the power is based on several averaging (in down-sampling and autocorrelation) so that this can depressthe random part of the signal more effectively than the Fourier transform. We investigate how the power of a resonance peak changes when increasing sampling interval (i.e. the averagingstep) m . The variation trend of the power is not simple. At small m and around the true peak, the firstterm in Eq. (18) that corresponds to the aliasing peak can be neglected. Also at small m the other modes butthe j -th mode can not make the resonance peaks around q ≈ q j . Thus the power of the resonance peak at q = Round ( q j ) can be approximately evaluated: | P Sm ( q j , m ) | ≈ √ n U j A j P k U k A k ≈ r Nm U j A j P k U k A k ∝ √ m U j A j P k U k A k , (25)which reminds us of the Wiener-Khinchin’s theorem that the power is proportional to square of the Fourieramplitude of the mode. In fact, in the case of m = 1 it holds A j = 1 so that the power is proportional to thesquare of the amplitude U j . In a mono-cyclic signal, i.e. where all k = j , the power of the resonance peak willdecrease approximately in inverse proportion to square root of the sampling interval √ m .However, in the multi-cyclic signal, the sampling (i.e. m ) dependency will be complicated due to that A j depends also on sampling interval m . Eq. (17) for A j can be rewritten: A j = 1 − exp (cid:18) πiq j n (cid:19) − exp (cid:18) πiq j mn (cid:19) − exp (cid:18) − πiq j n (cid:19) − exp (cid:18) − πiq j mn (cid:19) ≈ − exp (cid:18) πim ∆ tT j (cid:19) − exp (cid:18) πi ∆ tT j (cid:19) − exp (cid:18) − πim ∆ tT j (cid:19) − exp (cid:18) − πi ∆ tT j (cid:19) = sin (cid:18) πm ∆ tT j (cid:19) sin (cid:18) π ∆ tT j (cid:19) , (26)where Eq. (14) is used. A relation of A j vs. m is shown in Fig. 3(a). A j reaches maximum at m = T j t (i.e.the Nyquist sampling of the j -th mode) and the maximum value is A j max = sin − (cid:18) π ∆ tT j (cid:19) . (27)In denominator of Eq. (25) the sum P k U k A k can be replaced with a representative mode U k A k . A weightedmean mode or predominant mode can be the representative mode. Then Eq. (25) is reduced to a simple ratio A j A k , which is evaluated with approximation: A j A k ∝ sin (cid:16) πm ∆ tT j (cid:17) sin (cid:16) πm ∆ tT k (cid:17) ≈ T k T j (cid:16) πm ∆ tT j (cid:17) over m ∈ [1 , T j ∆ t ], if T j < T k T j T k (cid:16) πm ∆ tT k (cid:17) − over m ∈ [1 , T k ∆ t ], if T j > T k m ∈ [1 , T j ∆ t ], if T j ≈ T k (28)This expression is a rough estimation that, however, shows that the power of longer cycles than the representativemode will vary exceeding the proportion to √ m and the power of shorter cycle – exceeded by the proportion to √ m (Fig. 3). This can be understood by noting that a decay of short mode at its Nyquist sampling will let thepower of longer modes to increase due to the energy conservation of the signal, which belongs to the normalizedautocorrelation Eq. (10) where the power in Eq. (25) represents a relative fraction of every mode among thewhole energy of the signal. If the case is for a mono-cyclic signal or the predominant mode, the power decreasesin proportion to √ m . In Fig. 4 we can see the different sets of modes. If a short-period mode ( j = 1-st mode inFig. 4) is predominant, then the variation of the power for the longer modes(the j = 2-st mode in Fig. 4) willhave an elevation around the Nyquist sampling above the proportionality to √ m . If a long-period mode(the j = 3-rd mode in Fig. 4) is predominant, we see the opposite variation for the shorter modes. And for thepredominant cycle itself, its degradation when increasing the sampling interval is proportional to √ m .7 A j (a) A j A k (b) A j A k (c) A j A k (d) Figure 3:
The variation of A j in T j = 88 (a) and of the ratio A j A k in (b) T j = 88 , T k = 210, (c) T j = 210 , T k = 88and (d) T j = 92 , T k = 88 vs. sampling interval m . In all cases ∆ t = 1 and N = 200. In (a) the vertical bar stands forthe Nyquist sampling, m = T j t . In (b), (c) and (d) the solid lines stand for the exact evaluation (l.h.s. of ≈ ) and thedashed lines for the approximate evaluation (r.h.s. of ≈ ) for the ratio A j A k in Eq. (28). (a) (b) (c) Figure 4:
The variation of the power for the j = 2-nd mode in a composite signal consisted of three modes: T =95 , T = 210 , T = 330 , φ = 1 . , φ = 1 , φ = 2 . , ∆ t = 10 and N = 200 in the cases of (a) U = 10 , U = 2 , U = 3,(b) U = 1 , U = 10 , U = 3, (c) U = 1 , U = 2 , U = 10.. The variations of its amplitudes are different. The solid lineindicates numerical evaluation, the dotted line – for the approximate evaluation by Eq. (16), the dash-dotted line for thevariation proportional to √ m and the dashed line for the variation according to Eq. (25). The red vertical bar standsfor the Nyquist sampling m ∆ t = T . The peaks after the Nyquist sampling are nodal peaks. The samplogram of a multi-cyclic signal of U = 1 , U = 2 , U = 3 , T = 95 , T = 210 , T = 330 , φ = 1 . , φ =1 , φ = 2 . , ∆ t = 10 and N = 200. The solid blue line means the Nyquist sampling where the sampling interval is half ofthe period. The diagram extends to the Fourier limit where the period equals the sampling interval. The dashed blackhorizontal lines stand for the true ridges and the dashed red curved lines for aliasing ridges. A sampling peak appearsclearly at the Nyquist sampling. In this section we show the general configuration of the samplogram. We can also the D-samplogram and theeffect of differencing.
A samplogram is the diagram that shows how the power changes when the sampling interval increases in theaveraging down-sampling. The samplogram has 2 dimensions: resonance dimension (period or frequency) and sampling dimension (sampling interval). We can consider a third dimension, which is just the power. Asupright on the paper, the power can be represented by density or contours. Fig. 5 shows a samplogram for amulti-cyclic signal Eq. (1) and its sampled time series Eq. (2), where U = 1 , U = 2 , U = 3 , T = 95 , T =210 , T = 330 , φ = 1 . , φ = 1 , φ = 2 . , ∆ t = 10 and N = 200.In samplogram Fig. 5 the Nyquist sampling (solid blue) and the Fourier sampling limit (oblique rightboundary of the diagram) are visible. The samplogram is a version of the power spectrum in two dimensions.Also, the power spectrum is the Fourier transform of autocorrelation. In the Fourier transformation of a band-limited signal, the period that is considered is limited down to the sampling interval including the aliasingregion. In other words, once a sampling interval is given, cycles of period less than the sampling interval arenot considered. This sets the Fourier sampling limit to which the samplogram extends. The Nyquist samplingcorresponds to that the sampling interval is half of the period of cycle. The sampling theorem says that a cyclicsignal can be fully reconstructed only if it is sampled by an interval less than half of the period. As seen atthe end of § m = 1, that means no down-sampling and the ordinate, and theNyquist sampling is called a true periodicity region while the region between the Nyquist sampling and theFourier sampling limit is called an aliasing periodicity region . Ordinary power and Fourier spectra should dealonly with the true periodicity region while the samplogram shows both regions.In Fig. 5 two types of ridges are shown: true ridges and aliasing ridges. Both ridges correspond to extensionsof the two types of peaks: the true peak and the aliasing peak. When increasing the sampling interval, tracksof peaks form the ridges . Both ridges can be located in the frequency domain as follows: q = q j , for the true ridge (29) q = n − q j , for the aliasing ridge (30)where n is the length of the down-sampled time series u ¯ m and q j is the frequency number for the j -th mode.9n the period domain the above equations can be rewritten: T = mn ∆ tq j , for the true ridge (31) T = mn ∆ tn − q j , for the aliasing ridge (32)where we recall n = IntegerPart (cid:18) Nm (cid:19) ≈ Nm , (33) q j = Round (cid:18) mn ∆ tT j (cid:19) ≈ N ∆ tT j . (34)The same kind of ridges of the different modes do not cross each other: the true ridges of the j -th and l -thmodes as well as their aliasing ridges are parallel. However, the true ridge of a mode and the aliasing ridgeof another mode intersect each other. At their crossing a new kind of peak appears, which is called a nodalpeak . (See Fig. 6. Don’t confuse nodal peaks with discrete peaks in the aliasing ridges, which are related tothe sampling discretization.) The sampling peak can be considered a kind of nodal peak where the true andaliasing ridges of the same mode are overlapped. From Eq. (16) the power of nodal peak of the j -th and l -thmodes can be expressed P Sm ( q j , q l ) ≈ √ n P k U k A k (cid:20) U j A j exp (cid:18) πiq j n (cid:19) + U l A l exp (cid:18) − πiq l n (cid:19)(cid:21) . (35)The nodal peaks are seen in Fig. 2 as well as Fig. 4, which corroborates that the nodal peaks appears on theridge. From Eqs. (31) and (33) the nodal peak of the true ridge of the j -th mode and the aliasing ridge of the k -th mode happens at the sampling step m jl ≈ Nq j + q l . (36)We can see that this location is symmetric for q j and q l . Thus it is clear that the nodal peak of the true ridgeof the j -th mode and the aliasing ridge of the l -th mode and another nodal peak of the j -th aliasing ridge andthe l -th true ridge happen at the same sampling interval. We can also find the both nodal peaks have the samepower. The true ridge and aliasing ridge for the same cycle have the same degradation in terms of the samplinginterval. So, the both nodal peaks that appear at the same sampling will have the same power.The samplogram can explain the aliasing phenomenon that at lower sampling frequency than the Nyquistfrequency an aliasing peak appears in the power and Fourier spectra (e.g. Poluianov, 2014). In Fig. 5 it is visiblethat the aliasing ridge of the j = 3-rd mode invades into the true periodicity region from the aliasing region.This means that if the sampling rate is lower than the Nyquist rate for the j = 3-rd mode, then an aliasingpeak appears in spectrum while the true peak disappears so that the frequency of the j = 3-rd mode couldbe estimated faulty. In fact, the ordinary spectra deal with only the true periodicity region so this aliasingphenomenon was difficult to show. The samplogram shows this intuitively. l -samplogram In most physical processes, the difference series is considered to clarify the underlying dynamics. The Brownianmotion is an typical instance where the displacement but the position itself is the basic dynamic quantity. Thedifference series consists of differences between adjacent points of the time series:∆ u ( r ) = u ( r + 1) − u ( r ) for a time series u ( r ) , (37)where r is the index of point. If the original time series is a multi-cyclic one u ( r ) in Eq. (2), then its differenceseries can be written ∆ u ( r ) = X j π ∆ tT j U j cos (cid:18) πr ∆ tT j + φ j + π (cid:19) . (38)Here we can be sure that the periodicity in time series conserves by the differencing, which is just a differencingstability of cycle . And we can find that the amplitude of a mode decreases in inverse proportion to its period.10 (a) (b)
50 100 200 500 Period0123456Power (c) (d)
Figure 6:
The power spectra and D l (more precisely, ¯D l , i.e. the samplograms are evaluated by first averaging and laterdifferencing)-samplograms of the l -th order difference series of the multi-cyclic signal in Fig. 5. (a) and (b) show thepower spectrum and D1 (more precisely, ¯D1)-samplogram of the first-order difference series while (c) and (d) – powerspectrum and D3 (more precisely, ¯D3)-samplogram of the third-order difference series. In power spectra the longer cyclesare the more depressed. Dashed lines indicate the ridges. At crossing of ridges, the nodal peaks emerge. So the differencing can be considered a kind of long-period (low- frequency) filter. If the series is l th-orderdifferentiated, a stronger filter will be given:∆ l u ( r ) = X j (cid:18) π ∆ tT j (cid:19) l U j cos (cid:18) πr ∆ tT j + φ j + lπ (cid:19) . (39)According to the Wiener-Khinchin’s theorem, the power of a mode should decrease in inverse proportion to the2 l -th power of its period in the l -order differencing (see Fig. 6). If the period of a mode is long enough, itspower might become neglected at all. Furthermore, according to the regular degradation in Eq. (25), the ridgeswould be degenerated when increasing sampling interval. So, at the Nyquist sampling the power could be morenegligible than a resonance peak at m = 1.However, Fig. 6 (in particular (d)) shows that the ridges are elevating just before the Nyquist sampling forthe differences. Why are the ridges of negligible modes in power spectra elevated when approaching its Nyquistsampling? It can be explained from that the differencing is acting like a long-period (or low-frequency) filterwhile the Nyquist sampling acting like a short-period (or high-frequency) filter. Consider a mode, e.g. the j = 2-th mode in Fig. 6. At the Nyquist sampling of this mode, the shorter j = 1-st mode has already passedwell beyond its Nyquist sampling so the mode would decay at all. On the other hand, through the differencing,the longer j = 3-rd mode will be depressed much enough. Thus at the Nyquist sampling of the j = 2-th mode,only the j = 2-th mode itself will remain, in spite of that its power is negligible in power spectrum of m = 1 andthe 1-st and 3-rd modes are depressed. So, the power of the j = 2-th mode will be increased up at its Nyquistsampling. This fact can be explained in another way. As said in § √ m ). The effect of differencing is so strong that we can neglect theamplitude of modes in original time series and can reorder the power of modes merely by their period in somehigh-order differences. From this fact, the diagram can show the successive upheavals of ridges in the order oftheir period around their Nyquist sampling.The nodal peaks swell this elevation more around the Nyquist sampling. In fact, nodal peaks are distributedmuch more around the Nyquist sampling. Thus, the nodal peaks can exceed the sampling peaks (Fig. 7). Thenodal peaks distributed around the Nyquist sampling look like the Laue’s X-ray diffraction pattern that is apowerful tool to analyze the structures of the material so we could expect an information of signal from thosenodal peaks. 11 (a) (b) (c) (d) (e) (f) Figure 7:
The variations of power for modes in Fig. 5 in D1- and D2-samplogram along with sampling interval m .(a), (c), (e) correspond to the modes of j =1, 2 and 3 in D1-samplogram and (b), (d), (f) to the same modes in D2-samplogram. The vertical bars mean the nodal and sampling peaks. The peaks look to be shifted due to truncation ofseries, which is said in Eq. (4). The shortest j =1-st mode has shown a ridge degrading below the proportionality to √ m but longer modes have shown some elevation, in particular, just before the Nyquist sampling. Nodal peaks strengthensuch elevations. We can find that the higher the differencing order would be the shorter the period of the predominantcycle should be so that the more cycles should have the elevation around the Nyquist sampling. The samplogram accents the behavior of the power at the Nyquist sampling. This behavior looks rathermore obvious in the differences. We’d like to call the samplogram of l -th order difference series a D l -samplogram .According to experiences, the D2- or D3- samplogram will make the ridges elevating successively in the orderof their period around the Nyquist sampling.The D l -samplogram needs both differencing and averaging. Normal procedure consists of first differencingand later averaging. We can take the reverse order of operations: first averaging and later differencing. Thereverse procedure shows similar properties as the normal one. We denote this kind of samplogram as ¯D l -samplogram. Fig. 6 shows the ¯D1- and ¯D3-samplogram depicted by the reverse order of operations. Thedifference between D l - and ¯D l -samplograms appears for a stochastic signal while they appear to be the samefor a deterministic signal. We can enhance the visuality of the samplogram such as the resolution, contrast and so on by some sophisticatedtricks. We can also extend the samplogram analysis to the well-known Lomb-Scargle periodogram, which willgive us more advantages.
The resolution and contrast in the samplogram can be improved by using some sophisticated methods.The length of the down-sampled time series decreases inversely proportional to the sampling step m (Eq. 4): n = IntegerPart ( N/m ), which leads to lowering the resolution of period, in particular, at low frequencies. Indeed, using the expression for period in DFT (Eq. 12), the resolution of period can be defined∆ T ( q ) = T ( q + 1) − T ( q ) = N ∆ tq + 1 − N ∆ tq = N ∆ tq ( q + 1) ≈ N ∆ tq , (40)where q is the frequency index. If the length N is changed into N ′ , the frequency number q for the given period12 is changed: q ′ = N ′ ∆ tT = q N ′ N , (41)from which ∆ T ′ ( q ) ≈ N ′ ∆ tq ′ = NN ′ N ∆ tq = NN ′ ∆ T ( q ) (42)So the resolution of period is proportional to the length of the time series. To avoid lowering the resolution inthe down-sampling, the zero padding can be applied at the end of the autocorrelation series ACm ( g ) of everydown-sampled time series, which is called the post-zero padding. Padding any number of zeros at the end of asignal has no effect on its DFT (Orfanidis, 2010), while the length of the series increases and the resolution isimproved.Here we consider a question in detail: may the resonance peaks be shifted in spectrum by the zero padding?Evaluate the location of the resonance peak in the frequency and period domain. Recall the principal termsof the power of a sinusoidal signal as in Eq. (16). If the post-zero padding is applied to the series of theautocorrelation ACm ( g ) so the length of series, n , grows to n ′ , then the power is reevaluated as P Sm ′ ( q ) = 1 √ n ′ n ′ X g =1 ACm ( g ) exp (cid:18) πi ( g − qn ′ (cid:19) = 1 √ n ′ n X g =1 ACm ( g ) exp (cid:18) πi ( g − qn ′ (cid:19) ≈ √ n ′ P k U k A k X j U j A j " exp (cid:18) πiq j n (cid:19) − exp (cid:16) πi (cid:16) qn ′ + q j n (cid:17) n (cid:17) − exp (cid:16) πi (cid:16) qn ′ + q j n (cid:17)(cid:17) + exp (cid:18) − πiq j n (cid:19) − exp (cid:16) πi (cid:16) qn ′ − q j n (cid:17) n (cid:17) − exp (cid:16) πi (cid:16) qn ′ − q j n (cid:17)(cid:17) . (43)Then the resonance peak of m = 1 locates at qn ′ = q j n or q = q j n ′ n , (44)which can be rewritten in the period domain T ′ j = mn ′ ∆ tq = mn ∆ tq j = T j , and T ′ aj = mn ′ ∆ tn ′ − q = mn ∆ tn − q j = T aj , (45)where T j , T aj are the true and aliasing periods of the j -th mode. Thus we can conclude that the resonancepeaks or the ridges should not be shifted by the post-zero padding, while the resolution of period is improved.However, lengthening the time series will accompany lowering Fourier amplitude or power in inverse proportionto the length. It will be considered in next subsection.The zero-padding implies us an important thing: a non-stationary signal can make a peak in the spectrum.For example, if only one period or short segment of cycle is existing for a much longer time, this segment ofthe cycle can make a peak in spectrum, which means spurious cycle because the peak in spectrum means thestationary cycle while the segment means non-stationary. For example, if in some period of solar activity thereare only two grand minima, then the separation between them pretends to be a stationary cycle via Fourieranalysis, because this approach decompose any signal into the series of stationary sinusoidal signals. As seen in Eq. (25), the power of a cycle decreases inversely proportional to square root of the sampling step m , which makes it difficult to observe the behavior at the Nyquist sampling where the long cycles are muchlowered. Besides that, the post-zero padding also lowers the power inversely proportional to square root ofextended length n ′ . In order to visualize the power at the Nyquist sampling more clearly, it would be neededto lift the ridge for the great sampling interval. To do this, for example, a lifting proportional to m could beapplied to the power of every down-sampled series.Figure 8 shows the raw samplogram and visuality-enhanced samplogram with a extending (by zero padding)factor n ′ n = √ m and a lifting factor √ m . Manipulation begins with that the series of autocorrelation for every13 a) (b) (c) Figure 8:
The raw samplogram (a), visuality-enhanced samplogram (b) with an extending factor √ m and a lifting factor √ m and the peak-barred samplogram(c) with an extending factor √ m and a unit lifting factor. down-sampled series of length n = Nm is so zero-padded that its length becomes n ′ = N √ m , that is, n ′ − n zeros arepadded at the end of time series of autocorrelation. The length of the time series is extended by a factor √ m thatmay be user-defined, where N is the length of the initial time series before the down-sampling and m is averagingstep. Lowering power of every down-sampled series due to extending is compensated through multiplication bya factor p √ m = √ m to the down-sampled time series. Then, the power is multiplied additionally by the liftingfactor √ m that may also be user-defined. In Fig. 8, it is visible that the contrast of ridges against backgroundand the resolution of period for long cycle are enhanced, in particular, at the Nyquist sampling.As mentioned earlier, in DFT the peak appears at Round ( q j ) but not at the decimal frequency number q j ,where Round () means the rounding value. The period determined by Eq. (31) may vary when the samplinginterval varies. Thus, if we evaluate the period of cycle by a point-like peak, we might find a biased period byrounding, especially in a great sampling interval. In fact, if we found the peak at
Round ( q j ), the real frequencynumber could be in a range of [ Round ( q j ) − . , Round ( q j ) + 0 . Round ( q j ) − . , Round ( q j ) + 0 .
5) has all the samepower as the peak at frequency
Round ( q j ) (Fig. 8c). This will help to find the exact frequency for a cycle,especially for a low-frequency or long-period cycle that has the low resolution of period. The Lomb-Scargle (LS) periodogram is somewhat different from the aforementioned AC periodogram (Eq. 15).This can be applied to the unevenly sampled time series. Suppose that there are N data points u ( r ) = u ( t r ) , r =1 , , . . . , N . Then the LS periodogram of the time series can be defined by (Press et al., 2007) P N ( ω ) ≡ Var ( u ) ( [Σ r ( u ( r ) − ¯ u ) cos ω ( t r − τ )] Σ r cos ω ( t r − τ ) + [Σ r ( u ( r ) − ¯ u ) sin ω ( t r − τ )] Σ r sin ω ( t r − τ ) ) , (46)where ω is the frequency while ¯ u and Var ( u ) are the mean value and variance of u ( r ), respectively (Eq. ?? andEq. 9). And τ is defined by tan(2 ωτ ) = Σ r sin 2 ωt r Σ r cos 2 ωt r , (47)which implies the mean value of the time steps t r (Scargle, 1989). The LS periodogram was induced throughthe least-square method, so this has been often called the least-square (LS) periodogram.Eq. (46) appears to give the power at an arbitrary frequency ω so to be able to provide an high resolution.However, as (Scargle, 1989) indicated, the number of independent frequencies are limited to the length oftime series. So we can not expect higher resolution in LS periodogram in comparison with the classical ACperiodogram. However, LS periodogram has a peculiarity that this can be applied to the unevenly sampleddataset, so if we combine the samplogram and LS periodogram, that is if we depict the samplogram in termsof the peaks in LS periodogram, we can expand the samplogram to any unevenly sampled dataset. Figure 9shows the peaks in the samplogram of the AC periodogram and LS periodogram for the reconstructed sunspotnumber (RSSN) and the reconstruction of total solar irradiance (RTSI) . The dataset is available through the MPS sun-climate web-page at .See also Wu et al. (2018). The dataset is available through the NOAA web-page at .See also Steinhilber et al. (2012). yyyyyyyyyyyyyyyyyyyyyyyyy yyyyyyyyyyyyyyyyyyyyyyyyy yyyyyyyyyyyyyyyyyyyyyyyyy yyyyyyyyyyyyyyyyyyyy yyyyyyyyyyyyyyyy yyyyyyyyy yyyyyyy yyyyyyy yy yy y H yr L P e r i od H y r L (a) yyyyyyy yyyyyyy yyyyyyy yyyyyyyy yyyyyyy yyyyyyy yyyyyyy yyyyyy yyyyyyy yyyyyyy yyyy yyyyy yyy yyyy yy yy yyyy y y H yr L P e r i od H y r L (b) Figure 9: The peaks in the AC samplogram of RSSN within the period of 80 to 250 yr (a) and the peaks inLS samplogram for RTSI within the period of 450 to 850 yr (b) including D-samplograms upto the 3rd order.The marks “0”, “y”, “2” and “3” stand for the peak points in spectra of the 0-th to 3-rd order differences,respectively. The error bars mean the period resolution related to the rounding of frequency index (cf. the peakbarring in § So far we have seen only the deterministic periodic signal. The random signal is discussed in this section. Andthe effectivity of sampling and differencing stability analysis is discussed as well.
According to the idea of Fourier analysis, an arbitrary signal can be linearly decomposed into the series ofsinusoidal signals. So it is first important to analyze the sinusoidal signal that we saw in §
2. The feature of astationary cycle in the sinusoidal signal in the traditional spectrum is almost unique: • The resonance peak
The sinusoidal stationary cycle makes a resonance peak in the Fourier or power spectra.
The resonancepeak is a peak that appears in the Fourier or power spectra. The resonance peaks also appear in ordinateof the samplogram, i.e. at averaging step m = 1. As we know, however, many noisy and random peakspollute the spectra for a random signal and it is difficult to distinguish them.The samplogram can show the additional features of the sinusoidal stationary cycle, which we have seenbefore: • The ridge
The sinusoidal stationary cycle makes a ridge in samplogram.
The ridge is an extension of resonancepeaks along the sampling interval increasing. The ridge expresses the sampling stability of cycle.The sinusoidal signal has a regular degradation of the ridge in the averaging down-sampling.
Whenincreasing the sampling interval or sampling step m , the degradation of the longer cycle than the weightedmean or predominant cycle has an elevation in comparison with the proportionality to √ m , while for theshorter cycle a sinking appears. The predominant mode itself degrades in proportion to √ m . For a randomsignal, the randomness is added to this normal degradation, which threatens the existence of peak andridge.The elevation of the ridge for the longer cycle before the Nyquist sampling is very special. Peaks that arenegligible in the power spectrum become significant before the Nyquist sampling. This property appearsespecially in the D l -samplogram, where the predominant or weighted mean cycle is shifted to higherfrequency by differencing and successive elevation of ridges before the Nyquist sampling appears. We canadd a virtual “predominant and very short” cycle to the signal. Then the elevation of ridges for almostall of the cycles are observable including the probable shortest cycle. • The sampling peak 15 he sinusoidal stationary cycle makes the sampling peak at its Nyquist sampling.
This peak is formed bysuperposition of the true and aliasing ridges at the Nyquist sampling. However, in DFT, the appearanceof sampling peak is conditional: in the place of this peak a valley might appear.Incidentally, the sampling peak has another characteristic which is concerned with the complex argumentof power.
At the sampling peak, the power should be a negative real and must have the argument of ± π . In DFT, however, the argument could vary because the frequency number is given by an integer but adecimal number. The argument depends also on the parity of length of time series. • The nodal peak
In a multi-cyclic signal, different sinusoidal stationary cycles make the nodal peak in the samplogram.
Thenodal peak is formed by superposition of a true ridge for one cycle and an aliasing ridge for another cycle.The sampling peak is a kind of nodal peak, which is formed by the true and aliasing ridges for the samecycle.
The two nodal peaks for two cycles appear at the same sampling interval in both true and aliasingperiodicity regions and, furthermore, they have the same power.
Although the ridge for cycle may benegligible, the nodal peak on the ridge can appear significant. Like the sampling peak, the appearance ofthe nodal peak might be conditional.Here the sampling and nodal peaks are conditional, so the ridge can be regarded as a necessary propertyof the cycle. And the resonance peak is engulfed in the adjacent peak or background, it woulbe better for theexistence of cycle to be determined by ridge in samplogram rather than resonance peak in ordinary spectrum.In fact, a negligible cycle in sampling of m = 1 can become significant at its Nyquist sampling.The above features have been obtained for a sinusoidal stationary cycle. However, in practice we meet hardlywith such a signal. In fact, the solar activity is considered to be a non-stationary and stochastic process. If arandom signal is mixed with the deterministic cycles, then what will happen to the above features? Will theybe hardly detectable? We consider the random signal and effectivity of the samplogram to detect a weak cyclein it. Consider a multi-cyclic signal with noise: u ( r ) = U + U cos (cid:20) πr ∆ tT + φ (cid:21) + U cos (cid:20) πr ∆ tT + φ (cid:21) + + U cos (cid:20) πr ∆ tT + φ (cid:21) + w ( r ) , (48)where U = 10 , U = 5 , U = 1 , U = 1 , T = 420 , T = 50 , T = 12 , φ = 0 . , φ = 1 , φ = 0 . , ∆ t = 1 , and N = 800. The index r runs from N to N to avoid the initial transition. The first mode corresponds to the“whole time trend”, the second – a mid-term cycle and the third – a short-term cycle. A random signal w hasa distribution of N (0 ,
10) that is much stronger than the cycles. This signal and its Fourier and power spectraare shown in Fig. 10. We can see that the second and third modes of T = 50 , T = 12 seem to be almostspurious in spectra. Especially, in the Lomb-Scargle periodogram (Fig. 10(d)) the false-alarm probabilities forthose cycles are near to 0.999, which means that those cycle might be almost noisy and can not be considereddeterministic.We test the stability of cycles in the samplogram and D-samplograms. Figure 11 shows the peaks of thesamplogram and D-samplograms up to the 3rd order. We consider that the real decimal frequency is roundedto integer in DFT so we show the rounding region of a given peak within ± . § ∼
45 t.s. The weak cycle of 137 t.s. could be supposed. The 9.5 and 137-t.s. cycles are obviously thespurious cycles. In spite of appearance of the spurious cycles, the example shows that the stochastic stabilitycan filter many other spurious cycles effectively.The shorter region than 8 t.s. that we neglect here have only 3 or less sampling steps to check the stability.The shorter cycles have the higher resolution of the period. And these cycles almost correspond to “the noisecycles”. So we would meet the stronger requirement for the shorter cycles that the multiplet of the peaksof various orders should be consistent within the rounding region. However, the longer cycles have the widevariance of the peaks due to low resolution and long range of sampling interval to check the stability. Even theridge bending appears for the long cycle. The consistency of the peaks of various orders are expected only near16
00 600 700 800time - (a) A m p lit ud e (b) P o w e r (c) Pfa = = = = L SP o w e r (d) Figure 10: (a) A noisy multi-cyclic signal (Eq. 48), (b) its Fourier spectrum, (c) the power spectrum and (d) theLomb-Scargle periodgram. In (d) the “Pfa” indicates the false-alarm probability. yyyyy yyyyyyyy yyyyyyyyyy yyyyyyyy yyyy yyy y
222 222222222 22222222 222222 222 222 23333333 33333333 33333333 3333 33 3 P e r i od (a) yy yyyy yyyy yyy yyy yyy yyyyy yyy yyy y yy y y
222 222 2222 22222 222 2222 222 22 22 22 2 23333 333 33333 33333 3333 33333 333 333 33 33 3 P e r i od (b)
000 000 000 000 0000 0000 000 000 000 00 00 000 00 00 00 0 0 0 0 yy yy yyyy yy yy yyy yyy yy yyyy yyyy yyy yy yyy yy yyy yy y y y y y y
222 2222 222 2222 222 2 222 222 22 222 222 2222 2 22 22 2 2 2 2 233 333 33 333 333 3 33 33 333 333 333 33 3 33 33 3 3 3 3 P e r i od (c) yy y y y y y y y yy y y yy yy yy y y y y yy y yy y y y y y y y y P e r i od (d) Figure 11:
The peaks of (a) the samplogram, (b) the D1-samplogram, (c) the D2-samplogram and (d) the D3-samplogramof the noisy multi-cyclic signal (Eq. 48). The error bar stands for the rounding region of frequency index. The bluediagonal means the Nyquist sampling. The marks of “0”, “y”, “2”, “3” mean the peaks of the samplogram and D1 toD3-samplograms. a) (b) (c) (d) Figure 12: The histogram for the separations between grand extremes of reconstructed solar activity. (a) standsfor the grand minima and (b) for the maxima from Be reconstruction in Inceoglu et al. (2015), while (c) forthe grand minima and (d) maxima in Usoskin (2016).the Nyquist sampling. So we can require “the consistency” for the short cycle and “the persistency” for thelong cycle.If the cycle is significant, its stability is obvious because in spite of the degradation in the down-sampling theridge will remain significant. However, the example proves that the non-significant cycles have some stabilityagainst the down-sampling and differencing.
Beside manifesting the effectivity of the samplogram analysis, the above example also shows that the featuresof the sinusoidal stationary cycle would appear in the samplogram even for a random signal. This is becausethe random signal is decomposed into the sinusoidal signals in the Fourier analysis. However, among thosesinusoidal stationary cycles there are the spurious cycles. What is the origin of those spurious cycles? Ofcourse, they come from the random noise. But what we are interested in is what process the significant spuriouscycles are made through. If we could know this process, we could find the stable cycles more precisely. We canconsider several factors:First, the spurious cycles can come from the non-sinusoidal waveform. For example, the rectangular ortriangular signals cause many spurious cycles. The current SSN dataset accompanies the 11-yr cycle as wellas 5.5-yr cycle, which is related to its non-sinusoidal waveform (Petrovay, 2010). However, such harmonicssometimes become the essential argument to detect any cycle in astronomical phenomena (Koll´ath & Ol´ah,2009; Ol´ah et al., 2009) and are not bothering in our case of solar activity.Secondly, the spurious cycles can come from the random occurrences of the grand extremes such as grandminima and grand maxima. As we have seen in § Be) and Usoskin (2016). We plot the histogram for the separations with bin of 20-yr width. The low count(equal to 1) probably means the spurious cycle because this is non-iterative. However, if the count is no lessthan 2, we can expect the iterative appearances of the separation so this might mean a stable cycle.We can find the Suess/de Vriece cycle appear more frequently between the grand maxima rather than thegrand minima. This is consistent with the result obtained in Chol-jun & Jik-su (2020) where we verified thisresult with the superposed epoch analysis (SEA) method. Beside this cycle, we can see the 130, 170, 190, 250-yrseparations to be more iterative.
In fact, the power spectrum is based on the Fourier approach where an arbitrary signal is decomposed into thesinusoidal stationary cycles. So the random signal is also transformed into the series of sinusoidal cycles andis often regarded to have short period and weak power. This transformation should be, of course, arbitraryand non-stable in any operation on the signal such as the down-sampling, because the down-sampling includesthe averaging in itself. In deed, we can observe that the ridge of spurious cycle degrades much randomly anddoes not follow the normal degradation, which affects the maintenance of ridge in the averaging down sampling.18o we can expect that the spurious cycles representing the random signal should vary in operations such asdown-sampling and differencing.However, for the deterministic cycle we can assume two properties:- the (down-) sampling stability.and- the differencing stability.The (down-) sampling stability implies that the peak for the cycle must be kept when the sampling intervalchanges (see § § § §
1, the drawback of the LS periodgram comes from that peak reflects the frequency consistencyand the cycle amplitude. If we want to remove the significant spurious cycle, we should divide both factors.So we introduced the cycle-keeping operations: the down-sampling and differencing. What we want from theseoperations is the power variation of the cycles. As we have seen in § § So far we have seen the analysis of the stable cycles in the samplogram. However, there are some problems inthe analysis.First, the samplogram is effective to remove the high-frequency spurious cycles but not for the low-frequencyones. This is related to that the averaging suppresses only the high-frequency or short-period noise. We haveseen that long spurious cycles remain significant in the samplogram.Secondly, no quantity is developed still yet to evaluate the stochastic stability of the cycle. We want thequantity like the false-alarm probability to evaluate the stochastic significance. If there is no such a quantity,we have to perform some tedious and ineffective work in terms of examination with the naked eye to match theD-samplograms and analyze the sampling stability.Thirdly, we can see a ridge bending in the samplogram for some signal. The ridge bending means thatthe true ridge bends up (to the longer period) or down (to the short period) when the sampling interval isincreasing, in particular near the Nyquist sampling. This may be related to the averaging sampling itself or thesignal peculiarity. The former may come from that we have truncated the time series. The latter may comefrom the modulation in signal. The ridge bending makes a trouble to the precise estimation of the period.We can consider another way to evaluate the stability of cycle. And we can add a signal artificially. Ifwe add a shortest and significant signal, the all existing cycles including even short cycle will have the ridgeselevating around the Nyquist sampling.We see samplogram analysis for the solar activity in Chol-jun (2020).19
Conclusion
We introduced a samplogram analysis that can infer stable cycles by applying cycle-keeping operations to anysignal. The spurious cycle is unstable in those operations. In this paper, we considered a averaging down-sampling and differencing as such cycle-keeping operations. In a averaging down-sampling the long cycle isamplified and the short cycle is depressed while in differencing vice versa.The spectral analysis with stochastic significance cannot distinguish the spurious cycles effectively, becausethe peak in spectrum involves not only the frequency consistence but also the amplitude of the signal. Even non-frequency-consistent or non-stationary random rise and fall of signal can make a significant peak in spectrum. Weshowed that the random separation between the grand extremes in long-term solar activity can be misunderstoodas a significant stationary (spurious) cycle in traditional spectral analysis.The sampling and differencing stability can separate the amplitude from the frequency stability in peakheight in spectrum, at least weakly, through the amplitude variation in various operations. We saw that thismethod is more effective than a spectral analysis with stochastic significance in distinguishing the spurious andstable cycles.
Acknowledgement
K. Chol-jun has been supported by Kim Il Sung University during investigation.
References
Attolini, M., R., Galli, M., Nanni, T, 1988, Long and short cycles in solar activity during the last millennia. InStephenson, F.R. and Wolfendale, A.W., editors, Secular solar and geomagnetic variations in the last 10,000years, Dordrecht: Kluwer, 49-68Beer, J., Tobias, S., M., Weiss, N., O., On long-term modulation of the Sun’s magnetic cycle, Monthly Noticesof the Royal Astronomical Society, 473, 1596, http:/doi:10.1093/mnras/stx2337
Beutler, F., J., Leneman, O., A., Z., Random sampling of random process, Information and Control, 9, 325-346Cameron, R., H., Sch¨ussler, M., 2019, Solar activity: periodicities beyond 11 years are consistent with randomforcing, Astronomy and Astrophysics, 625, A28, https://doi.org/10.1051/0004-6361/201935290
Chol-jun, K., Jik-su, K., 2020, Solar activity cycle of ∼
200 yr from mediaeval Korean records and recon-structions of cosmogenic radionuclides, Monthly Notices of the Royal Astronomical Society, 492, 384-393, http://doi:10.1093/mnras/stz3452
Chol-jun, K., 2020, Stable solar cycles, submittedEddy, J., A., 1976, The Maunder minimum, Science, 192, 1189-1202Inceoglu, F., Simoniello, R., Knudsen, M., F., Karoff, C., Olsen, J., Turck-Chi´eze, S., Jacobsen, B., H., 2015,Grand solar minima and maxima deduced from Be and C: magnetic dynamo configuration and polarityreversal, Astronomy and Astrophysics, 577, A20, https://doi.org/10.1051/0004-6361/201424212
Koll´ath, Z. and Ol´ah, K., 2009, Multiple and changing cycles of active stars. I. Methodsof analysis and application to the solar cycles, Astronomy and Astrophysics, 501, 695-702, http://dx.doi.org/10.1051/0004-6361/200811303
Lomb, N., R., 1976, Least-Squares Frequency Analysis of Unevenly Spaced Data, Astrophysical and SpaceScience, 39, 447-462.Ma, L. H., Vaquero, J. M., 2009, Is the Suess cycle present in historical naked-eye observations of sunspots?,New Astronomy, 14, 307-310, https://doi.org/10.1016/j.newast.2008.04.001
Ogurtsov, M., G., Nagovitsyn, Yu., A., Kocharov, G., E., Jungner, H., 2002, Long-Period Cycles of the Sun’sActivity Recorded in Direct Solar Data and Proxies, Solar Physics, 211: 371-394,Ol´ah, K. et al. 2009, Multiple and changing cycles of active stars. II. Results, Astronomy and Astrophysics,501, 703-713, http://dx.doi.org/10.1051/0004-6361/200811304
Petrovay, K., 2010, Solar Cycle Prediction, Liv. Rev. Sol. Phys. 7, 6,
Poluianov, S., Usoskin, I., 2014, Critical Analysis of a Hypothesis of the Planetary Tidal Influence on SolarActivity, Solar Physics, 289: 2333, http://arxiv.org/abs/1401.3547v1
Press, W., H., Teukolsky, S., A., Vetterling, W., T., Flannery., B., P., 2007, Numerical Recipes, The Art ofScientific Computing, 3rd ed., Cambridge Univ. Press, Cambridge, pp. 685-692.Scargle, J., D., 1989, Studies in Astronomical Time Series Analysis. III. Fourier Transforms, Autocorrelationand Cross-correlation Functions of Unevenly Spaced Data, Astrophysical Journal, 343, 874-887.Solanki, S. K., Usoskin, I. G., Kromer, B., Sch¨ussler, M., Beer, J., 2004, Unusual activity of the sun duringrecent decades compared to the previous 11,000 years, Nature, 431, 1084Steinhilber, F., Abreu, J. A., Beer, J., Brunner, I., Christl, M., Fischer, H., Heikkil¨a, U., Kubik,P. W., Mann, M., McCracken, K. G., Miller, H., Miyahara, H., Oerter, H., Wilhelms, F., 2012,9,400 years of cosmic radiation and solar activity from ice cores and tree rings, PNAS, 109, 5967,
Usoskin, I. G., 2008, A History of Solar Activity over Millennia, Living Review of Solar Physics, 5, 3,
Usoskin, I. G., Gallet, Y., Lopes, F., Kovaltsov, G. A., Hulot, G., 2016, Solar activity during the Holocene: theHallstatt cycle and its consequence for grand minima and maxima, Astronomy and Astrophysics, 587, A150, http://arxiv.org/abs/1602.02483v1
Wu, C. J., Usoskin, I. G., Krivova, N., Kovaltsov, G. A., Baroni, M., Bard, E., Solanki, S. K., 2018, Solaractivity over nine millennia: A consistent multi-proxy reconstruction?, Astronomy and Astrophysics, 615,A93, https://doi.org/10.1051/0004-6361/201731892https://doi.org/10.1051/0004-6361/201731892