If the power spectrum is always of the form $1/k^p$ there are some optimisations you can make. The first is to compile a function to generate the reciprocal of the spectrum: $k^p$. This can be very fast since you don't need to check for the $k=0$ point (the reciprocal spectrum is well behaved at $k=0\ $).

We can then use ReplacePart to temporarily set the value at $k=0$ to something other than 0 (I've used 1), then take the reciprocal to get the spectrum we really wanted, and finally use ReplacePart again to set the value at $k=0$ back to 0.

The other optimisation is to create the random noise directly in frequency space, by adding a random phase to the normally distributed amplitudes. Doing this allows you to get 2 fields at once, from the real and imaginary parts of the complex field which is generated.

The code below is only for 2D, but I think it should generalise to 1D or 3D quite straighforwardly.

The timing to generate 2 fields at 1024 x 1024 is about 280 ms on my PC.

fastpower=Compile[{{n,_Integer},{p,_Real}}, Block[{temp}, temp=Range[-n,n-2,2]^2; Outer[Plus,temp,temp]^(p/2)]]; spectrum[n_,gamma_]:=ReplacePart[1/ReplacePart[ RotateRight[fastpower[n,-0.5gamma],{n/2,n/2}], {1,1}->1.],{1,1}->0.]; complexrandoms[n_]:=RandomReal[NormalDistribution[0,1],{n,n}] (Cos[#]+I Sin[#])&[RandomReal[{0,2\[Pi]},{n,n}]]; fields[n_,gamma_]:={Re[#],Im[#]}&@Fourier[complexrandoms[n]spectrum[n,gamma]]; 

Example:

{f1,f2}=fields[1024,-3]; GraphicsRow[Image@Rescale@#&/@{f1,f2}] 

If the power spectrum is always of the form $1/k^p$ there are some optimisations you can make. The first is to compile a function to generate the reciprocal of the spectrum: $k^p$. This can be very fast since you don't need to check for the $k=0$ point (the reciprocal spectrum is well behaved at $k=0\ $).

We can then use ReplacePart to temporarily set the value at $k=0$ to something other than 0 (I've used 1), then take the reciprocal to get the spectrum we really wanted, and finally use ReplacePart again to set the value at $k=0$ back to 0.

The other optimisation is to create the random noise directly in frequency space, by having independent normally distributed real and imaginary components. Doing this allows you to get 2 fields at once, from the real and imaginary parts of the complex field which is generated.

The code below is only for 2D, but I think it should generalise to 1D or 3D quite straighforwardly.

The timing to generate 2 fields at 1024 x 1024 is about 280 ms on my PC.

fastpower=Compile[{{n,_Integer},{p,_Real}}, Block[{temp}, temp=Range[-n,n-2,2]^2; Outer[Plus,temp,temp]^(p/2)]]; spectrum[n_,gamma_]:=ReplacePart[1/ReplacePart[ RotateRight[fastpower[n,-0.5gamma],{n/2,n/2}], {1,1}->1.],{1,1}->0.]; complexrandoms[n_]:=RandomReal[NormalDistribution[0,1],{n,n}]  + I RandomReal[NormalDistribution[0,1],{n,n}] fields[n_,gamma_]:={Re[#],Im[#]}&@Fourier[complexrandoms[n]spectrum[n,gamma]]; 

Example:

{f1,f2}=fields[1024,-3]; GraphicsRow[Image@Rescale@#&/@{f1,f2}] 

If the power spectrum is always of the form 1/k^p there are some optimisations you can make. The first is to compile a function to generate the reciprocal of the spectrum: k^p. This can be very fast since you don't need to check for the k=0 point (the reciprocal spectrum is well behaved at k=0).

We can then use ReplacePart to temporarily set the value at k=0 to something other than 0 (I've used 1), then take the reciprocal to get the spectrum we really wanted, and finally use ReplacePart again to set the value at k=0 back to 0.

The other optimisation is to create the random noise directly in frequency space, by adding a random phase to the normally distributed amplitudes. Doing this allows you to get 2 fields at once, from the real and imaginary parts of the complex field which is generated.

The code below is only for 2D, but I think it should generalise to 1D or 3D quite straighforwardly.

The timing to generate 2 fields at 1024 x 1024 is about 280 ms on my PC.

fastpower=Compile[{{n,_Integer},{p,_Real}}, Block[{temp}, temp=Range[-n,n-2,2]^2; Outer[Plus,temp,temp]^(p/2)]]; spectrum[n_,gamma_]:=ReplacePart[1/ReplacePart[ RotateRight[fastpower[n,-0.5gamma],{n/2,n/2}], {1,1}->1.],{1,1}->0.]; complexrandoms[n_]:=RandomReal[NormalDistribution[0,1],{n,n}] (Cos[#]+I Sin[#])&[RandomReal[{0,2\[Pi]},{n,n}]]; fields[n_,gamma_]:={Re[#],Im[#]}&@Fourier[complexrandoms[n]spectrum[n,gamma]]; 

Example:

{f1,f2}=fields[1024,-3]; GraphicsRow[Image@Rescale@#&/@{f1,f2}] 

If the power spectrum is always of the form $1/k^p$ there are some optimisations you can make. The first is to compile a function to generate the reciprocal of the spectrum: $k^p$. This can be very fast since you don't need to check for the $k=0$ point (the reciprocal spectrum is well behaved at $k=0\ $).

We can then use ReplacePart to temporarily set the value at $k=0$ to something other than 0 (I've used 1), then take the reciprocal to get the spectrum we really wanted, and finally use ReplacePart again to set the value at $k=0$ back to 0.

The other optimisation is to create the random noise directly in frequency space, by adding a random phase to the normally distributed amplitudes. Doing this allows you to get 2 fields at once, from the real and imaginary parts of the complex field which is generated.

The code below is only for 2D, but I think it should generalise to 1D or 3D quite straighforwardly.

The timing to generate 2 fields at 1024 x 1024 is about 280 ms on my PC.

fastpower=Compile[{{n,_Integer},{p,_Real}}, Block[{temp}, temp=Range[-n,n-2,2]^2; Outer[Plus,temp,temp]^(p/2)]]; spectrum[n_,gamma_]:=ReplacePart[1/ReplacePart[ RotateRight[fastpower[n,-0.5gamma],{n/2,n/2}], {1,1}->1.],{1,1}->0.]; complexrandoms[n_]:=RandomReal[NormalDistribution[0,1],{n,n}] (Cos[#]+I Sin[#])&[RandomReal[{0,2\[Pi]},{n,n}]]; fields[n_,gamma_]:={Re[#],Im[#]}&@Fourier[complexrandoms[n]spectrum[n,gamma]]; 

Example:

{f1,f2}=fields[1024,-3]; GraphicsRow[Image@Rescale@#&/@{f1,f2}]