A test for the absence of aliasing or white noise in two-dimensional locally stationary wavelet processes

Either intentionally or unintentionally, sub-sampling is a common occurrence in image processing and can lead to aliasing if the highest frequency in the underlying process is higher than the Nyquist frequency. Several techniques have already been suggest in order to prevent aliasing from occurring (for example applying anti-aliasing filters), however there is little work describing methods to detect for it. Recently, Eckley and Nason (Biometrika 105(4), 833–848, 2018) developed a test for the absence of aliasing and/or white noise in locally stationary wavelet processes. Following Eckley and Nason (Biometrika 105(4), 833–848, 2018), we derive the corresponding theoretical consequences of sub-sampling a two-dimensional locally stationary wavelet process and develop a procedure to test for the absence of aliasing and/or white noise confounding at a fixed point, demonstrating its effectiveness and use through appropriate simulation studies and an example. In addition, we outline some possibilities for extending these methods further, from images to videos.


Introduction
When sampling a time series at a fixed rate, if the Nyquist frequency is lower than the highest frequency in the process, aliasing occurs, since high frequency components cannot be distinguished from low frequency ones due to the sampling rate not being high enough. This can lead to unreliable analyses and distorted spectral and auto-covariance estimates. For example [see Priestley (1983)], when sampling a real valued stationary process X (t) with corresponding non-normalised spectral density function h X (ω) at intervals t, we obtain a new process Y t = X (t · t) with corresponding nonnormalised spectral density function Hence, if the sampling rate is too low, without prior knowledge, high frequency information will be indistinguishable from low frequency components in the spectrum. Images may be intentionally down-sampled for memory or processing purposes. For example, convolutional neural networks often down-sample images in various layers, the effects of which are explored in Ribeiro and Schön (2021), Vasconcelos et al. (2021), Vasconcelos et al. (2020), with Vasconcelos et al. (2021), Vasconcelos et al. (2020) also describing the benefits of applying anti-aliasing filters at certain locations within the network architecture. Recently, Li et al. (2021) even suggest applying wavelets to convolutional neural networks in order to counteract the effects of aliasing. Therefore, it seems that aliasing can definitely impact the performance of these models. As another example, consider autonomous vehicles. These systems often have to carry out many tasks simultaneously and have a limited amount of memory and computing power available (Liu et al. 2019(Liu et al. , 2021. Due to these limitations, the sampling rate when collecting a video feed might be set too low or one might adjust the sampling rate in order to free computing power for other tasks. As a result, aliasing may be induced, which in turn may have a serious impact if critical high frequency information is confounded. The ability to detect aliasing would allow these systems to adjust their sampling rate accordingly. Aliasing can be prevented if one knows the highest frequency in the process prior to sampling it, but this is rare in practice. Alternatively one could apply anti-aliasing filters (as in the convolutional neural network example above (Vasconcelos et al. 2021(Vasconcelos et al. , 2020) that ensure the sampling rate is high enough in order to capture all the high frequency information in the process. However, this is not always applicable and is not guaranteed to prevent aliasing from occurring. Aliasing and possible methods to prevent it are described in more detail in Priestley (1983), Gonzalez and Woods (2018), Burger and Burge (2010). Since it is not always possible to prevent aliasing, developing a test to detect for it becomes crucial.
Some methods for detecting aliasing in images have already been developed. For example, Hinich and Wolinsky (1988) use the bispectrum (Hinich and Messer 1995) of a signal to formulate a test for aliasing by using a modified version of the Hinich bispectrum test for Gaussianity. More recently, Coulange and Moisan (2010) suggest using the Fourier domain of an image in order to develop an aliasing detection algorithm. Alternatively, Patney and Lefohn (2018) and Douglas (2012) propose using deep neural networks and non-linear factor analysis respectively to detect for aliasing a sequence of images. Our approach differs in that we use locally stationary wavelet processes and are not constrained to using a sequence of images when carrying out our test.
Locally stationary wavelet processes were first introduced in Nason (2000) and were used in order to develop a test for the absence of aliasing and/or white noise confounding in locally stationary wavelet time series in Eckley and Nason (2018). In this paper, we extend these methods to two dimensions with a specific application to images.
PhD thesis (Gott 2012) proposed extending (Eckley and Nason 2018), but only did so for the case of Shannon wavelet processes, assumed that identical dyadic sub-sampling was carried out in both directions and also used an invalid testing regime (which was subsequently fixed in Eckley and Nason (2018)). Use of Shannon wavelets causes considerable simplification of the mathematics underlying the problem, but the assumption is probably far too strong in practice as the wavelets are infinite in extent. Likewise, assuming identical sub-sampling in horizontal and vertical image directions is strong, when one is usually uncertain on how the image has been processed before acquisition. Our work extends the methodology to Daubechies' (Daubechies 1992) compactly supported wavelets, permits different dyadic sub-sampling in horizontal and vertical directions and provides a rigorous testing framework inherited from Eckley and Nason (2018). This paper is organised as follows: Sect. 2 provides a review of the necessary background material on locally stationary wavelet processes, Sect. 3 derives the theoretical consequence of sub-sampling on two-dimensional locally stationary wavelet processes and the raw wavelet periodogram, Sect. 4 outlines the test procedure, Sect. 5 demonstrates the test using simulation studies, Sect. 6 gives an example and finally Sect. 7 outlines some possibilities for extending this procedure from images to videos.

Locally stationary wavelet processes
We review the main concepts from the theory of locally stationary wavelet processes as described by Nason (2000). These were later extended to two-dimensions in Eckley (2001). For more information on locally stationary wavelet process, please see Nason (2000), Eckley and Nason (2018), Eckley (2001) or the book by Nason (2008).
Definition 1 (One-dimensional discrete mother and father wavelets). Let {h n } and {g n } be the low and high-pass quadrature mirror filters used in the construction of the Daubechies (1992) compactly supported orthogonal continuous time wavelets. The one-dimensional discrete mother and father wavelets at scale j, denoted by ψ j = (ψ j,n ) n=0,...,N j and φ j = (φ j,n ) n=0,...,N j respectively, are N j -dimensional vectors obtained using the following formulae where N j = (2 j −1)(N h −1)+1 and N h denotes the number of non-zero elements of {h k }.
Definition 3 (Two-dimensional locally stationary wavelet process). A two-dimensional locally stationary wavelet process is a sequence of stochastic processes {X r ,s }, r = 0, . . . , R, s = 0, . . . , S, with the representation where {ξ j,u,v } are uncorrelated mean zero random variables with unit variance, {ψ j,u,v } are two-dimensional discrete mother wavelets and {w j,u,v } are amplitudes satisfying the following conditions: ∀ j ∈ N, and = h, v, d there exists a Lipschitz continuous function W j : the Lipschitz constants L j , corresponding to the functions W j , are uniformly bounded in j and and ∞ j=1 2 j L j < ∞, (iii) there exist {C j } such that ∞ j=1 C j < ∞ and for each R and S, where, for each j ∈ N and = h, v, d, the supremum is taken over u = 0, . . . , R − 1 and v = 0, . . . , S − 1.
Definition 5 (One-dimensional discrete autocorrelation wavelets). The one-dimensional discrete autocorrelation wavelet and scaling functions are given by respectively, for τ ∈ Z and j ∈ N.
The inner product operators A = (A j,k ) j,k∈N , B = (B j,k ) j,k∈N and C = (C j,k ) j,k∈N are J × J dimensional matrices constructed from the autocorrelation wavelets as follows Although not strictly necessary for defining for defining two-dimensional locally stationary wavelet processes, these inner product operators will be useful in Sect. 3 when studying the effects of dyadically sampling such a process.

Remark 2
We note that, for real valued discrete wavelets, although A = A T and B = B T , it is not generally true that C = C T .

Definition 7 (Two-dimensional raw wavelet periodogram).
The raw wavelet periodogram of Y r ,s is defined as I j,u,v = (d j,u,v u, v ∈ Z and = h, v, d, where d j,u,v are the non-decimated wavelet coefficients of Y r ,s given by Definition 8 (Two-dimensional discrete autocorrelation wavelets). The two-dimensional autocorrelation wavelets are given by where j ∈ N, and = h, v, d.
Definition 9 (Two-dimensional inner product operator). The two-dimensional inner product operator A (2) is a 3J × 3J dimensional matrix constructed from the autocorrelation wavelets as follows where and 1 ≤ j, k ≤ J .
Remark 3 Note that the one-dimensional inner product matrices A, B and C are related to the two-dimensional inner product matrix A (2) as follows where 1 ≤ j, k ≤ J .

Theoretical consequences of dyadic sub-sampling
Following the work in Eckley and Nason (2018), we now examine the effects of dyadic sub-sampling on a twodimensional locally stationary wavelet process. We note that some work on the two-dimensional case is presented in Gott (2012), however, only processes represented using Shannon wavelets and the effects of dyadic sub-sampling that is the same in both directions is considered. Hence, in this work we generalise these concepts by deriving results that are applicable to any choice of wavelets and to dyadic sub-sampling that is not necessarily the same in both dimensions. Analogously to the one-dimensional case (Eckley and Nason 2018), dyadic sub-sampling of a two-dimensional locally stationary wavelet process results in a new process defined on a smaller grid which is the sum of two separate components: a locally stationary wavelet process and a process which is asymptotically white noise. The variance of the latter depends on the S 1 (z), which may not be constant, at the points that remain after the sub-sampling has occurred, in essence, the information at the finest scale that was lost. This is formalised in the following theorem: Theorem 1 Let X r ,s be a two-dimensional stationary wavelet process with evolutionary wavelet spectrum given by ,s or (c) Y r ,s = X 2r ,2 s , the {Y r ,s } r ,s∈Z can be decomposed as Y r ,s = L r ,s + F r ,s where L r ,s is a two-dimensional locally stationary wavelet process with the same underlying wavelets as X r ,s , with raw wavelet periodogram expectation given by (24), (25) and (26) for directions = h, = v and = d respectively, and F r ,s is a zero mean process and autocovariance By repeatedly applying the above theorem to a process, we obtain: Corollary 1 Let X r ,s be a two-dimensional stationary wavelet process with evolutionary wavelet spectrum given by If Y r ,s = X 2 p r ,2 q s , then it can be decomposed as Y r ,s = L r ,s + F r ,s , where L t is a twodimensional locally stationary wavelet process with the same underlying wavelets as X r ,s , with raw wavelet periodogram expectation given by (24), (25) and (26) for directions = h, = v and = d respectively, and F r ,s is a zero mean process and autocovariance cov (F r ,s Our next result shows what happens to the expectation of the raw wavelet periodograms for the vertical, horizontal and diagonal wavelet directions.
Theorem 2 Let X r ,s be a two-dimensional stationary wavelet process with evolutionary wavelet spectrum given by

with Daubechies compactly supported wavelets. Without loss of generality, let p ≤ q (the case p ≥ q is the same but with the expectations of the horizontal and vertical raw wavelet periodograms switched).
Then the expectations of the raw wavelet periodograms of Y r ,s = X 2 p r ,2 q s are given by Note that the first terms in the above expressions are equivalent to the folded spectrum for stationary processes as stated in equation (1). Due to the non-orthogonality of the autocorrelation wavelets, we obtain additional terms consisting of information at finer scales that has been lost due to subsampling and which reappears as confounded information at coarser scales.

Remark 4
Note that if p = q then the expectation of the raw wavelet periodogram is given by for i = h, v, d, which is analogous to the result for onedimensional locally stationary wavelet processes derived in Eckley and Nason (2018).
To aid our comparison, we write the expectation of the raw wavelet periodogram for the original process using inner product operators given in Definition 5: where Comparing these to the expectations of the raw wavelet periodograms for processes that have undergone dyadic subsampling, we note a few key differences. The first term on the right hand side of (24), (25) and (26) appears at all scales and locations and is the same for all directions. This corresponds to power at finer scales from all wavelet directions re-appering at coarser scales due to dyadic sub-sampling in both dimensions at the same rate. The second term in each of (24), (25) and (26) is also surplus, however in this case it differs between directions -the terms in the sum are weighted by the one-dimensional inner product operators from Definition 5. Note also that these are shifted by p. This corresponds to power at finer scales being redistributed to coarser scales due to additional dyadic sub-sampling in only one dimension (in this case the first dimension, since p ≥ q). Finally, similarly to the one-dimensional case (Eckley and Nason 2018), the terms on the right hand sides of (24), (25) and (26) have their inner product matrices shifted by p in one dimension and q in the other and now lie on the grid defined by (2 p u/R, 2 q /S).
As described in Eckley and Nason (2018), a locally stationary wavelet process that is white noise with variance σ 2 has a raw wavelet periodogram expectation of σ 2 at all scales and locations, which means that we cannot distinguish between white noise or aliasing. In essence, if spectral power appears at all scales at a given location than white noise confounding or aliasing has occurred, although we can't tell them apart. We formalise this in the next section.

Testing for the absence of aliasing or white noise confounding
Let's assume that we would like to test whether there is evidence of aliasing or white noise confounding at a point We test whether aliasing occurred in each wavelet direction = h, v, d independently.
In an ideal situation, we would have access to the true underlying evolutionary wavelet spectrum S j at all scales and so we would test null hypothesis against the alternate hypothesis for ∈ {h, v, d}. However, as described by Eckley and Nason (2018), we have to estimate S j from the realisation of the process and, due to boundary effects, only obtain reliable estimates for a finite set of scales j ≤ J † < J . Therefore, we would have to test the null hypothesis against the alternate hypothesis where in this case, S j is the estimated wavelet spectrum of the process. As noted in Eckley and Nason (2018), if H 0 is true then automatically H (I ) 0 is true as well. The same cannot be said for the alternate hypothesis, since there could exist a scale j > J † such that S j (z 0 ) = 0. For this reason, the test becomes a test for the absence of aliasing or white noise confounding in a locally stationary wavelet process.
Following Eckley and Nason (2018), we describe how the test can be carried out in practice. Let b ∈ N >0 , define the set and collect the sample {Ŝ j (z)} z∈ , whereŜ j (z) denotes the estimate of the wavelet evolutionary spectrum at the point z computed as described in Eckley (2001). In order to do so, we first obtain an estimate of the raw wavelet periodogram I j (z), which is a biased estimate of the evolutionary wavelet spectrum. We then use the inner product operator A from Definition 9 to correct for this bias, obtaining the unbiased estimatorŜ j (z). Note that the above set corresponds to the set of spectral estimates lying on a 2b + 1 × 2b + 1 square grid centred at z 0 . From the Lipschitz continuity assumption on S j (z), we have that and so under H 0 we have that ∃ j ∈ {1, . . . , J † } such that which is approximately zero for small u and v and large R and S. Therefore we would like to test whether the mean of the sample {Ŝ j (z)} z∈ is equal to zero for each scale j ∈ {1, . . . , J † }. To do this we can use a Student's t-test on the sample {Ŝ j (z)} z∈ with test statistic where S j,z 0 andσ j,z 0 are the sample mean and standard deviation respectively. Since we are testing multiple hypotheses simultaneously, we use Holm's method (Holm 1979) to ensure that the family-wise error rate is not too high.
As the size of the image increases, we can estimate the evolutionary wavelet spectrum more reliably for a larger number of scales and so it is only natural to increase the number of scales J † to consider in the test. The choice of b also varies with the size of the image. As the latter increases, a constant value of b means that the sample becomes increasingly localised around the point z 0 . In essence, there is a trade-off between size and localisation of the sample-a larger sample (and b) loses local information. Hence, with larger images one can choose a larger value of b if desired.
As previously noted, rejecting the null hypothesis does not guarantee that aliasing or white noise confounding definitely occurred, but rather suggests that this may have happened. Furthermore, the test does not distinguish between the two cases. Hence, further investigation would need to be carried out. For example, one could sample the process at a higher rate and check the finer scales of the corresponding wavelet spectrum, however, it is not always possible to carry this out in practice.

Dealing with correlation in the samples
As noted in Eckley and Nason (2018), due to the Lipschitz continuity assumption, the samples at each individual scale are likely to contain high autocorrelation between observation in close proximity to one another. In this case, the Student's t test will perform poorly as the independence assumption is violated. Some possible solutions to this problem are discussed in Zwiers and Storch (1995) and Thiebaux and Zwiers (1984). A very simple solution is to sub-sample the data; Zwiers and Storch (1995) argue that, in the case of temperature data, observations separated by 5 days from one another could be assumed to be independent. The Student's t test would then be carried out on this smaller sample. Since we are testing data that lies on a regularly spaced square grid, this method can be easily adapted by choosing observations that lie a certain distance from one another on this grid. However, this requires that a significant proportion of the data be discarded, resulting in a loss of information, a concern that was also raised in Zwiers and Storch (1995). Furthermore, one would have to determine an appropriate distance between observations that are kept for the test.
Another approach is to compute the equivalent sample size (ESS) (Zwiers and Storch 1995;Thiebaux and Zwiers 1984). This can be interpreted as the number of effectively independent observations in the sample. For spatial data lying on an n × n grid, the equivalent sample size is defined as where ρ(τ 1 , τ 2 ) is the spatial autocorrelation function with lags τ 1 and τ 2 (a full derivation of this result, following the one-dimensional one in Thiebaux and Zwiers (1984), is available in the appendix). However, in most cases, this cannot be reliably estimated from the data. Some improvements can be obtained by truncating the sum in (40) after a certain number of terms (Zwiers and Storch 1995;Thiebaux and Zwiers 1984). Further improvements to the estimate can also be obtained by constraining the equivalent sample size to values that are deemed to be realistic, as suggested by Zwiers and Storch (1995), and so for an n × n grid, Since we are testing a small region of the spectrum at each scale, we have a few possible ways for estimating the ESS. The natural way of doing so, is to estimate the ESS directly from the sample we are testing. Alternatively we could estimate the ESS for the entire spectrum at a given scale and then re-scale this to the number of observations contained in our sample. Since we are now using a significantly larger number of observations, our estimate of the spatial autocorrelation function (and hence the ESS) will be more robust. The cost of doing so is that regardless of the region we are testing, the ESS will be the same and hence we may lose some local information.

Computational feasibility, parallelisation and speed up
In this subsection, we first address the estimation of the evolutionary wavelet spectrum, focusing on the number of computations required to perform this operation and ways in which this can be parallelised. Consider the two-dimensional process {X m,n }, defined on the square grid m, n = 0, . . . , N .
To compute an estimate of the evolutionary wavelet spec-trum, several steps are involved. Firstly, the non-decimated wavelet transform of the data {X m,n } needs to be computed, as explained in greater detail in Nason (2008) and Eckley (2001). Using the pyramid scheme due to Mallat (1989), this can be performed in O (N 3 log N ) operations. The resulting transform is then squared to obtain the raw wavelet periodogram, requiring a further O(N 2 log(N )) operations. However, due to the inherent bias of the wavelet periodogram, it is necessary to correct it using the inner product operator from Definition 9. This correction yields an unbiased estimate of the evolutionary wavelet spectrum and has a computational complexity of O(N 3 log N ), since we are performing a matrix multiplication at each scale. For further information, please refer to Nason (2000), Nason (2008) or Eckley (2001). Therefore, overall, the computational complexity of this procedure is O (N 3 log N ). Now we turn to the possibility for parallelisation. We can parallelise the non-decimated discrete wavelet transform. As the two-dimensional wavelet transform is composed of sequences of one-dimensional transforms performed on the rows and columns of the data, these sequences can be distributed across multiple processing units, allowing for simultaneous computation. A more comprehensive description of this approach can be found in Marino et al. (1999) and Chaver et al. (2002). Additionally, Nielsen (1998) provides insights into this approach and discusses methods for parallelising the one-dimensional wavelet transform as well.
Following this, parallelising the computation of the wavelet periodogram is quite simple. As it involves squaring the non-decimated wavelet coefficients, this operation can be easily divided among multiple processors.
We now turn our attention to the correction of the raw wavelet periodogram, which involves multiple matrix multiplications, specifically the inverses of the inner product operators with the wavelet periodogram. This is quite straightforward, as matrix multiplications are simply sequences of inner products between the rows of the first matrix and the columns of the second, which can easily be distributed over multiple processors. Numerous existing procedures already address this step, and a detailed discussion is beyond the scope of this section. Comprehensive coverage can be found in Schatz et al. (2016) and Gunnels et al. (1998), for example. By parallelising each step involved in the estimation of the evolutionary wavelet spectrum, notable advancements in computational speed can be achieved.
When testing for aliasing in a two-dimensional locally stationary wavelet process, the natural approach would be to apply the test at each point. However, in the case of images, conducting such tests on every individual pixel can be computationally infeasible, especially for large-scale data sets.
To mitigate this computational burden, an alternative approach is to test regions or specific pixels in the image instead of individually testing every pixel. This approach Fig. 1 Test spectrum for S 1 (z 1 , z 2 ), = h, v, d takes advantage of the fact that testing for the absence of aliasing and/or white noise confounding at a point necessitates analysing the entire region around it. By selecting representative regions or points, the computational workload can be significantly reduced while still providing reliable insights into the presence or absence of aliasing and/or white noise. This sampling strategy enables faster computation, making the test more feasible for practical applications involving large data sets. Additionally, parallelisation techniques can be applied to the testing process itself by distributing the test computations across multiple pixels at the same time.

Simulation study
The test-bed spectrum used in our simulation study is defined in Eq. (42) and Fig. 1.
We begin by conducting the experiment on realisations without sub-sampling and test the points z 0 = (1/2, 1/2) and z 0 = (1/4, 1/4). For the true spectrum defined above, at these points, we expect to accept H 0 and reject H 0 in favour of H A respectively for all directions. We drew 1000 realisations of size 1024 × 1024 from the test-bed spectrum using Daubechies D5 wavelets and test each wavelet direction l = h, v, d independently for scales j = 1, . . . , 4 using the method described above and a nominal size of 5%. The results are displayed in Tables 1 and 2. As we can see, for We also note that in Table 2, our rejection rates for large window widths could be viewed as being too high, and we attribute this to the difficulty in correcting for sample correlations when running the tests. Our next experiment uses the same test-bed spectrum as before, except that in this case the entire spectrum in the horizontal direction is set to zero as well as the spectrum at scale j = 3 in the diagonal direction. We again drew 1000 realisations of size 1024×1024, however in this case we subsample in both dimensions, reducing these to processes of size 512 × 512. As before, we test the points z 0 = (1/2, 1/2) and z 0 = (1/4, 1/4) and for each wavelet direction independently, with the only difference being that we only test j = 1, 2, 3 (since we lost a scale due to sub-sampling). In this case for the true underlying spectrum we would reject H 0 for = v and accept H 0 for = h and d at z 0 = (1/2, 1/2) and reject H 0 for all wavelet directions at z 0 = (1/4, 1/4). The results are summarised in Tables 3 and 4. When testing the point z 0 = (1/2, 1/2), we see that in most cases our test does not reject the null hypothesis in the horizontal and diagonal directions. This is to be expected, as the true spectrum for = h is set to zero for all scales j and for = d it is set to zero for j = 3. Since in all directions there is no spectral power at z 0 = (1/2, 1/2) for j = 1, we expect no power to be redistributed at finer scales due to subsampling. On the other hand, since there is spectral power at z 0 = (1/4, 1/4) for = v and d at the finest scale prior to sub-sampling, power is redistributed at all coarser scales including all scales in the horizontal direction which used to have no spectral power. Although for b = 32 and b = 64 we would reject H 0 , further investigating would need to be carried out in order to conclude that aliasing occurred with certainty. Furthermore, we note that our rejection rates for b = 32 in the horizontal and diagonal directions are significantly lower than in the vertical direction, which is due to power already being present in this direction, as well as difficulties in correcting the test for correlation in the samples.

A real world example
In this section we demonstrate the use of the test on a real image of size 2048×2048. We tested for the absence of aliasing/white noise confounding at each pixel for J † = 5 using Daubechies D6 wavelets. Computational time overall was quite long, taking approximately a week in total. However a significant amount of that computational time was used to compute the evolutionary wavelet spectrum (which we did only once) since the software we used was not parallelised. We note therefore that efficient parallelisation of the software would reduce computational times significantly, for example by following the suggestions outlined in Sect. 4.3. The image and results are displayed in the top row of Fig. 2, corresponding to vertical, horizontal and diagonal directions respectively. Pixels where we would reject the null hypothesis H 0 are displayed in gray. We note that this is rejected more frequently at pixels located in rougher region of the image, for example the trees at the center, compared to the smother regions such as parts of the sky and ground. If aliasing did occur, then this is to be expected, as rougher patches of the image contain more information, making them susceptible to aliasing, whereas smooth regions are constant in value and therefore do not lose information when down-sampled.
As pointed out in previous sections, although we can conclude that aliasing/white noise confounding did not occur in the white regions of the image, we cannot conclude with certainty that it did occur in the gray regions. Further investigation would be required to reach such a conclusion and, if reached, to determine which of the two (aliasing or white noise confounding) occurred.
We now sub-sample the original image twice in both the horizontal and vertical dimensions, reducing it to an image of size 512 × 512. We carry the test out again using the same underlying wavelets and for J † = 3. Due to the smaller size of the image, the computational times are significantly lower, taking only a few hours. The sub-sampled image and the corresponding results are displayed in the lower row of Fig. 2, where again pixels for which the null hypothesis H 0 was rejected are coloured in gray.
The test is rejected for a proportionally larger number of pixels which correspond to rougher patches containing more information. Since we have access to the original version of the image, we can conclude that aliasing occurred in areas that were not previously coloured in gray.

Extension to videos
A natural extension to the work in this paper is to develop a similar procedure for testing for the absence of aliasing/white noise confounding in videos. We can distinguish between two scenarios: the first involves a fixed camera and lighting, while the second assumes variability in both. In the former case, each pixel in the video consistently represents the same object throughout the entire duration, eliminating the need for specific modifications. Conversely, in the latter case, while each pixel may potentially correspond to a different object across different frames, it is worth noting that it may not be strictly necessary to perform adjustments in all cases. This presents a more intricate challenge in determining the pixel-to-object correspondence. One possible avenue for addressing this situation is to explore methods that compensate for camera movements relative to its initial position and adapt to changes in lighting. However, it is important to note that such an approach would require precise information about camera movements and lighting conditions, which may not always be readily available, and therefore alternative solutions will need to be explored. Consequently, it becomes evident that the challenge of handling a moving camera with changing lighting requires careful consideration.
We now briefly outline two possibilities for going from images to videos. The first is to treat videos as threedimensional processes and therefore treat them as threedimensional locally stationary wavelet processes. Testing for the absence of aliasing/white noise confounding would require extending the methods outlined in Eckley and Nason (2018) and this paper to three-dimensional processes.
Unlike time series or images, the dimensions in a video are not all of the same kind and instead consist of one temporal and two spatial dimensions. When dealing with spatio-temporal data, it is common practice to handle the spatial and temporal components as distinct entities, assuming that this separability is valid. For further insights into spatio-temporal statistics, please refer to Cressie and Wikle (2011) for example. This leads naturally to our second suggestion, which is to apply both the one-dimensional and two-dimensional tests to a given pixel at a given time.
More formally, given a point z 0 in the video, this can be viewed as both a point in the time series of a pixel and as a pixel in an image corresponding to that frame. Hence, one could apply the one-dimensional version of the test developed in Eckley and Nason (2018) to the point in time series and the two-dimensional version of the test to the point in the image simultaneously, obtaining results (i.e. whether to accept or reject the null hypothesis) in time and in the tree wavelet directions. To then determine whether there is evidence for aliasing or white noise confounding at that point, one could then apply a multiple hypothesis scheme. This and other possibilities are left for future work.

Closing remarks
In this paper we have extended the absence of aliasing or white noise confounding test developed in Eckley and Nason (2018) to two-dimensional locally stationary wavelet processes and demonstrated it's application to images. Additionally, we gave some suggestions for extending these methods further and, in particular, a procedure for applying the test to videos.
A clear avenue for future work would be to investigate solutions for dealing with correlation in the samples in more detail. Although the methods outline in this work are able correct the Student's t test for correlation, they are still unreliable at times. Further work could also be conducted on choosing an appropriate window width in which to conduct the test. Finally, it would also be interesting to see these methods being applied to videos as well, as this is a natural extension to the work developed here.
Author contributions HAP and GP wrote the main manuscript text. HAP prepared the figures and tables.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your  Koentjoro (2018) intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Proof of Theorem 1
The proof is similar to the one-dimensional case in Eckley and Nason (2018), although it requires some adaptations to two-dimensions.
(a) We begin by substituting s with 2s into Eq. (11) to obtain where L r ,s corresponds to the j > 1 terms and F r ,s to the j = 1 term, i.e. u,v ψ j,u−r ,v−2s ξ j,u,v and We focus first on the term L r ,s , which we can decompose as where (e, e) denotes terms that are evenly-indexed in both components, (e, o) denotes terms that are evenly-indexed in the first component and oddly-indexed in the second, and so on. Consider where we assume, without loss of generality, that r is even. Note that if r was odd, then the expression on the right hand side of (47) would simply equal L (o,e) r ,s . Now, using the definition of the two-dimensional mother wavelet, we have where in (50) we summed over m = r + a and n = s + b instead of a and b and in (51) we summed over p = u − m and q = v − n instead of m and n. Substituting (51) back into (47), we obtain where Summing over a = u − p and b = v − q instead of p and q we have w j,2(a+ p),2(b+q) ψ j−1,a−r ,b−s ξ j,2(a+ p),2 (b+q) . (56) Define ξ * j−1,a+ p,b+q = ξ j,2(a+ p),2(b+q) , W * j−1 (z) = W j (z) and w * j−1,a+ p,b+q = w j,2(a+ p),2(b+q) for j = 2, 3, . . . a, b ∈ Z and z ∈ (0, 1) 2 . Then From assumption (iii) in Definition 3, we have and so U ( p,q) r ,s is itself a two-dimensional locally stationary wavelet process. Since only a finite number of the h 2 p+r and h 2q terms are non-zero, L (e,e) r ,s is equivalent to the sum of a finite number of two-dimensional locally stationary wavelet processes with constant coefficients, which does not depend on r or s, and so is a two-dimensional stationary wavelet process as well. We can proceed in a similar manner for L (e,e) r ,s , L (e,e) r ,s and L (e,e) r ,s to conclude that L r ,s is a two-dimensional locally stationary wavelet process with the same underlying wavelets as X r ,s .
We now compute the covariance of F r ,s . This process has mean zero since ξ j,u,v has mean zero for all j, u, v and . Now, We include the last step, namely (62) where in (64) we used the fact that ψ h 1,u,v = h u g v . We can obtain similar expressions for the vertical = v and diagonal = d directions as well. We now consider two cases: firstly, if w 1,u,v = w 1 = is constant, then, for = h, we have cov(F h r ,s , F h r +τ 1 ,s+τ 2 ) We obtain analogous results for = v and = d, and so cov(F r ,s , F r +τ 1 ,s+τ 2 ) = S h where in (70), we sum over m = u −r and n = v −2s instead of u and v. (71) follows from assumption (iii) in Definition 3 and in (72) we used the Lipschitz continuity of the process. Again, this is analogous for = v and = d, and so cov(F r ,s , F r +τ 1 ,

Proof of Theorem 2
We begin by computing the expectation of the wavelet periodogram of the process Y r ,s = X 2 p r ,2 q s , for i = h, v, d. k, a, b, j, u, v, i, p, q) (80) j, u, v, i, p, q) (81) k, . . . , p, q) j, u, v, i, p, q) where in (77) we substitute in the definition of the twodimensional locally stationary wavelet process, in (80) we set P ( , k, a, b, j, u, v, i, p, q and in (81) we sum over m = a − 2 p u and n = b − 2 q v instead of a and b. In (82) we use assumption (iii) from Definition 3 and in (83) we make use of the Lipschitz continuity assumption. Now, where in (89) we used the definition of the two-dimensional autocorrelation wavelet and in (90), we sum over a = x − r and b = y − s instead of x and y. Therefore Assuming without loss of generality that p < q and considering the case i = h, we have Note that here we used the result derived in Eckley and Nason (2005) that k (2 q b) = k−q (b) for k > q and k (2 q b) = δ 0,b for k ≤ q. An analogous property can be derived for the father autocorrelation wavelet j , which we also used above. Furthermore, we can decompose the expectation as follows    . As we noted previously, the case where p > qsimply results in the expectations for the vertical and horizontal directions being switched.

Derivation of the Spatial Equivalent Sample Size
Let X i, j , i, j = 1, . . . , n be observations on an n × n equally spaced grid, with mean μ and variance σ 2 . Following the derivation in Thiebaux and Zwiers (1984), where C(τ 1 , τ 2 ) and ρ(τ 1 , τ 2 ) are the symmetric spatial covariance and autocovariance functions, respectively, at lags τ 1 and τ 2 . In the case where we have n e independent observations, we would have σ 2 X = σ 2 /n e . Equation the two, we obtain the following expression for the equivalent sample size n −1 e = n −2 n−1 τ 1 ,τ 2 =−(n−1)

Appendix B -Additional tables of results
See Tables 5, 6, 7 and 8.