# Resampling with and without replacement for estimating significance of spectral components

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).

@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.

• Thank you very much for the reference. Do you think you can expand a little bit on what the OP is asking to see how bootstrapping is relevant to this question with the data it provides? We are talking about what is mentioned here specifically, that is, shuffling $x[n]$ and using the shuffled sequence for statistical testing. Not bringing in surrogate data. I don't see where is the "population" here to apply the bootstrapping to. Where are you going to randomly sample from?
– A_A
Mar 8 '18 at 16:41
• Hi: Not a lot room here. Take an AR(1). you construct an. AR(1) model and can't test for significance of coefficient ( say $\phi$ ) because underlying distribution is not known. So, take the original observations and sample with replacement (so that you have the same # of observations ) and estimate the same AR(1). do that again, and again, then you have many estimates of $\phi$. to estimate it's distribution . then, find where the actual estimate falls in the estimated distribution. if it's in the tails, that supports significance. google "intro to boostrapping" for clearer explanation. Mar 8 '18 at 20:25
• Note that although above was just in time domain, the same "concept" appliues any time you want to test for significance. Mar 8 '18 at 20:26
• Thank you for the explanation. I am not convinced that it applies entirely here. What do you train the AR(1) on? Random values of $x[n]$? And then, how does it work for spectral estimation and how is this different from a straightforward test using the average power of a noise sequence (?). These are the points I do not get from this answer.
– A_A
Mar 9 '18 at 12:06
• Hi: You estimate the AR(1) over and over on the shuffled (with replacement) observations. This will give you the distribution of the phi estimate under the null that there is no relationship. Then you estimate phi using the actual ( not shuffled ) data and compare that estimate to the estimated distribution. I'm pretty confident that it applies but the details definitely matter. When time allows, I'll read his question again and see if I can be more specific. I may not be able to though because I may not be able to understand what he's doing exactly because I'm new to DSP. Mar 9 '18 at 16:57

...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.

1. 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.

• This one will be more applied but you need to have access to Jstor. jstor.org/stable/2336278?seq=1#page_scan_tab_contents Mar 8 '18 at 16:24
• My apologies, but I am slightly confused as to what your suggestion for the correct way is? You refer to a 'deterioration', but each resort is entirely random (uniform pdf). I am only resorting multiple times in order to acquire a reliable averaged random spectrum against which to compare my spectral components. Are you a) saying that the easiest way to do this would be to compare it against a random process of the same energy? and b) that this is still incorrect? If so, I agree that my method is more time consuming, but why does this make it incorrect? Mar 10 '18 at 15:43
• P.S. My time series data are also very short - only 30 time points. A single white noise process of this length does not look very random! Mar 10 '18 at 15:52
• @MattVowels 1)Additional explanation clearer than original.2)Shuffling will take any $x[n]$ and convert it to white noise.White noise has equal magnitude across all frequencies.Why do you have to go through all of this shuffling when you can simply divide the power of the signal to the number of components that you have and then compare their amplitude? 3) For Welch's method, you can put confidence intervals in this way...
– A_A
Mar 11 '18 at 23:28
• @MattVowels ...4) Testing multiple $X[k]$ has to be adjusted. I would not trust it. It's invalid, the components are not independent. Can I please ask what is your data source (?) (in the time domain). 30 points (in time domain (?)) are very few anyway.
– A_A
Mar 11 '18 at 23:30