An algorithm is described that can generate random variants of a time series or image while preserving the probability distribution of original values and the pointwise Holder regularity. Thus, it preserves the multifractal properties of the data. Our algorithm is similar in principle to well-known algorithms based on the preservation of the Fourier amplitude spectrum and original values of a time series. However, it is underpinned by a dual-tree complex wavelet transform rather than a Fourier transform. Our method, which we term the Iterated Amplitude Adjusted Wavelet Transform (IAAWT) method can be used to generate bootstrapped versions of multifractal data and, because it preserves the pointwise Holder regularity but not the local Holder regularity, it can be used to test hypotheses concerning the presence of oscillating singularities in a time series, an important feature of turbulence and econophysics data. Because the locations of the data values are randomized with respect to the multifractal structure, hypotheses about their mutual coupling can be tested, which is important for the velocity-intermittency structure of turbulence and self-regulating processes.