Making own function to use in fourier transform ?

I want to create 10,000 sample (random phase but same power spectra) of my original time series using fourier spectrum analysis. The method I want to opt is given in figure. The actual time series data is attached.

so far I understood that I first need to use discrete power spectrum (note: method is different than default method the matlab command uses).

I am new to fourier spectral analysis and any sort of help would be really helpful to me. I look forward to your reply on the method of generating surrogate data using the method specified in figure.

Thanks in advance.