Resampling with and without replacement for estimating significance of spectral components

Question

In order to test the significance of spectral components, it seems reasonable to randomly sort the data in order to destroy all the serial correlations / spectral order e.g. 100,000 times, and then estimate the significance or C.I.s by comparing your sample spectrum against the spectra of all the resorted data. Essentially H0 = there is no serial correlation / spectral order, and therefore the signal is white noise.

Here the resorted data would be resampling without replacement. My question is - is this correct or incorrect for spectral estimation? Would it be better to use resampling WITH replacement? And, if so, then why?

---

Thanks all for your responses! Indeed I am aware it is called bootstrapping - and just to clarify, my intention is not only to resort the time series 10,000 times, but to also create 10,000 spectra from these resorted versions of the time series in order to establish how many times my spectral component of interest is above the resorted time series spectra.

So H0 is that, over 10,000 comparisons, my spectral component is not higher in amplitude than the corresponding component in the resorted spectra more than 5% of the time. My H1 being that the spectral component is greater in amplitude (and therefore not due to chance) in at least 5% or more of the comparisons.

Indeed, I think you can go further than this and also calculate how many times it is BELOW the resorted time series spectra - this would, presumably, give you an indication of its persistence (i.e. it is more persistent than white noise and is therefore unlikely to change in a given period).

The code below calculates: the spectrum for a sine wave X[k] (length 30 time points), an average spectrum for the sine wave resampled 20000 times without replacement, an average spectrum for the sine wave resampled 20000 times with replacement. The original signal spectrum could as well be an average spectrum of, say, 200 signals of interest.

 Fs=1; %Sampling Frequency
 T=30; %Total signal duration in seconds
 t=0:(1./Fs):(T-(1./Fs)); %Time vector
 f1=1; % A frequency figure (in Hz)
 x=cos(f1*t);
 X=fft(x);
 X = abs(X);

 for i = 1:20000
     ps = round(rand(1,30)*29)+1;
     x_p_all(i,:)=(x(ps))';
     X_hat(i,:) = abs(fft(x_p_all(i,:)));
 end

 X_hat_mean = mean(X_hat);

 for i = 1: 20000
     ps2 = randperm(length(x));
     x_p_all2(i,:)=(x(ps2))';
     X_hat2(i,:) = abs(fft(x_p_all2(i,:)));;  
 end

 X_hat_mean2 = mean(X_hat2)

 plot(X_hat_mean)
 hold on
 plot(X)
 plot(X_hat_mean2)

It can be seen that the resampling with replacement averages to be the same as the resampling without replacement average (as expected). Obviously, the energy of an individual resampled signal will be different for with and without replacement, but bootstrapping undertakes the process thousands of times. Now, if we don't permute and take the average, we have the option to compare our spectrum against each individual permutation, where each comparison would be a comparison like this:

Clearly, in this instance, the spectrum of the resampled signal (in this case, without replacement) looks far from white noise because it contains spurious periodicities as a result of the short time duration and limited number of frequency bins. My method compares these resampled signal spectra against the original signal spectrum e.g. 10,000 times to establish significance - a similar method to bootstrapping (which does not, as far as I am aware, require Bonferroni correction).

Answer 1

...is this correct or incorrect for spectral estimation?

From what is provided in this description, I am more inclined to say that it is incorrect.

Here is my understanding of what you are describing:

You have some discrete time series $x[n]$ and its spectrum $X[k]$ which I will assume here that you obtained via the Discrete Fourier Transform (DFT) as $X = \mathcal{F}(x)$. $n$ denotes discrete time, $k$ denotes discrete frequency. Now, you pick some spectral component $X[k]$ and ask yourself: "Could this component have emerged by chance or is it really part of the signal?".

To answer this question, you shuffle $x[n]$ through a series of permutations that swap around some $n+u$ sample to $n+v$ sample. After shuffling a sufficiently long time, what you end up with is some $\hat{x}(n)$, accompanied by some $\hat{X}(k)$, which is essentially a white noise sequence using the same values of $x(n)$.

And then what? Where are the data for your $H_0$?

I am assuming that you are interested on the significance of some $X[k]$ (?).

If you create some distribution of spectral components from $\hat{X}(k)$, you could answer the question "Does $X[k]$ come from this distribution?". That is, could it be because of chance?

BUT, that would be equivalent to considering the average amplitude of a white noise sequence with equal power as your $x[n]$. Then, if $X[k]$ is above that average you can be sure that it cannot be due to noise.

Ideally, what you would like to have had is a distribution for each X[k]. That is, many realisations of the process that generates your signal, many realisations of their spectra and from all of these, one distribution for each $k$. Then, you could answer the question "Given the value of this $X[k]$, does it really come from my signal or not?"

If you are worried about the variance of the components, you might want to use Welch's Method.

Would it be better to use resampling WITH replacement?

That would modify the permutation rules to allow or disallow $u=v$. I don't see where is the re-sampling here. You had one spectrum before, you have one spectrum after the re-shuffling...Except if you mean that you would be storing each and every permutation result (?) in which case, yes, that would be wrong because it would represent the statistics of the deterioration process....Or maybe I am getting the whole thing wrong (?).

By the way, there is another thing you need to be careful of. If you did a number of statistical tests to test each and every spectral component then you would be assuming that spectral components are independent of each other and also you would have to correct your scores for repeated tests.

Hope this helps.

EDIT:

A few notes with illustrations to understand better the problem and how to deal with it:

Let's build a signal to work with:

Fs=8000; %Sampling Frequency
T=4; %Total signal duration in seconds
t=0:(1./Fs):(T-(1./Fs)); %Time vector
p=2.*pi.*t; % Phase vector (These last two vectors save us typing the same things again and again).
f1=440; % A frequency figure (in Hz)
%
%
x=cos(f1.*p); %The signal
ps=randperm(length(x)); %A random permutation
x_p=x(ps); %The permuted signal
X = fft(x); % The spectrum of the signal
X_hat = fft(x(ps)); %The spectrum of the permuted signal

The key "objects of interest" here are x (the signal) and ps which is a permutation sequence. The permutation sequence will "loop" through the samples of x and re-map them to different time instances.

Randperm samples WITHOUT replacement. If you sample WITH replacement, then the power of the permuted sequence will be different to the power of the non-permuted sequence (see below). Even if you pass x through ps once, you will get a totally randomised x but you can try doing it a few times anyway.

Here are spectra of x and x[ps[n]] (i.e. x permuted once):

Here is the code that produces that graph:

subplot(221);plot(t,x);xlabel("Discrete time (s)");title("x[n]");grid;
subplot(222);plot(t,x_p);xlabel("Discrete Time (s)");title("Permuted x (x[ps[n]])");grid;
subplot(223);semilogy(Fs.*(0:length(x)-1)./length(x),(X.*conj(X))./(length(x).^2));grid;title("X[k]");xlabel("Discrete Frequency (Hz)");
subplot(224);semilogy(Fs.*(0:length(x)-1)./length(x),(X_hat.*conj(X_hat))./(length(x).^2));xlabel("Discrete Frequency (Hz)");title("Spectrum of permuted x[n]");grid

Notice the normalisation by length(x). Also notice how x is permuted once by x(ps) which effects the "re-mapping".

Is the power of x different than the power of x(ps). No, but only if you sample without replacement, in which case, each instantaneous value of x is re-used exactly once.

If you sample with replacement, some values of x might be re-used (re-map the same time instance twice), in which case the sum will work out to be different. You can try it with some other permutation sequence like ps = round(rand(1,length(x)).*length(x)); which as you can see uses rand which doesn't care about its "history" of random values produced.

Here is the first thing that I am trying to say / clarify from what you are writing:

If you pick one X[k] and compare it against $\hat{X}$ which is the whole spectrum of the permuted sequence (fig above), then:

1. What you would really be asking would be "Can this component have risen by chance from a random sequence with equal power as my x?". Irrespectively of whether that bottom right spectrum is the result of one permutation or more it will be the spectrum of white noise of a sequence of power equal to that of x. Then, why go to the trouble of permuting x a few times? What is the added value? What is the power of x? The power of x is sum(x.^2). What is the average power? The average power is sum(x.^2)./length(x). How many points are there in the FFT? There are length(x) points. What is the sort of power we expect each bin to take? (sum(x.^2)./length(x))/length(x) = sum(x.^2)./(length(x).^2). This is what the bottom two graphs depict (adjusted for one / two sided spectra, so, basically, each line divided by two).

This is why I am writing above, if you want to proceed this way, then you can simply say, the closer my spectral line gets to that "theoretical" average white noise value, the less I can be sure that it is due to signal.

2. You mention confidence intervals. If you proceed with the above frame of mind, you would be putting confidence intervals on spectral lines obtained by a distribution across the spectrum. In other words, you would put confidence intervals on $\pm X[k]$ from $\hat{X}$ (?). If you want to put confidence intervals on $X[k]$, you need to look within $X[k]$. In other words, put confidence intervals on the variation of $X[k]$ itself. How can you do this? Instead of deriving your $X$ from $x$, you derive it from successive overlapping frames obtained by $x$. In this way, you form $X_0[k], X_1[k], X_2[k], X_m[k]$ where $m$ is the $m^{th}$ frame and $k$ is the discrete frequency component. So, now, combined with the assumption that $x$ is Gaussian and that the power is squaring a Gaussian therefore the distribution will be Chi squared (the DFT is linear so no further changes), you can put confidence intervals from a Chi-square distribution with the degrees of freedom equal to the number of frames minus one. (Which, by the way, is the way pwelch does it and why I am hinting at Welch's method above).

That is all I am trying to say.

With 30 time instances, you might be able to average "across subjects" (?). So, if your study has 30 participants, you can derive the "normal" variation of each $X[k]$ from those. Then, for an unknown spectrum, you can use those distributions to judge if a specific harmonic is "expectable" there or not.

Hope this helps.

EDIT 2:

The above code can be used to experiment further, especially with re-shuffling. In Octave (and possibly MATLAB too), you can apply the reshuffling more than once in two ways, either via x(ps)(ps)(ps)(ps)(ps) and so on, or x(ps(ps(ps(ps(ps(ps(ps))))))) which permutes the permutation of the permutation of the permutation [...] of x. I have a sneaking suspicion that the first way doesn't work in MATLAB (that is the x(ps)(ps)...(ps). Of course, you can even do it with different permutation sequences every time too.

Answer 2

@Matt: In the statistical world, what you are doing is called bootstrapping. You need to sample with replacement so that the distribution converges to the true distribution (under the null of independence ). There is a ton of literature on bootstrapping ( Efron, 79, I think ) so I won't try to explain it but what you are doing seems reasonable to me. I would google for applied related papers. Below is relevant but it's from Annals so probably pretty theoretical-non applied.

https://projecteuclid.org/download/pdf_1/euclid.aos/1176348515