Resampling with and without replacement for estimating significance of spectral components
...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:


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. 


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.