In python: missing randomness (fractional Brownian motion)

I am new to Python. I have two scripts for generating and plotting a 2D lattice which values represent a spatially distributed attribute in the form of a fractal image. The first script contains the function generating the fractal (orthodox name: Spectral Synthesis with Inverse Fourier Transform), while the second simply calls the function, plots it, ad writes it to a .csv file. The function relies on multiple random.gauss(0.0, sigma) calls for obtaining normally distributed random numbers. It also calls random.random() since it needs other random numbers. Below are the codes. Thanks to anyone who will help!

The problem: the plot always shows the same pattern, hence making me think that there is something NOT random in the code. My intention is to obtain a different fractal image every time I run the code. What makes it NOT random? Is that seed=122 in the second script? PS: please note that since the procedure introduces a given degree of periodicity, and since I do not want that periodicity to be shown, the plotted lattice is smaller than the generated lattice. In particular, I am using size(plotted)=101x101, while size(generated)=256x256.

First script: the function itself

from __future__ import division #Avoids the floor of the mathematical result of division if the args are ints or longs
import numpy
import random
import math
import pylab

def SpectralSynthesisFM2D(max_level, sigma, H, seed=0, normalise=True, bounds=[0,1]):
  """
  ________________________________________________________________________
  Args:
      max_level : Maximum number of recursions (N = 2^max_level)
      sigma     : Initial standard deviation
      H         : Roughness constant (Hurst exponent), varies form 0.0 to 1.0
      H=0.8 is a good representation of many natural phenomena (Voss, 1985)
      seed      : seed value for random number generator
      normalise : normalizes the data using bound
      bounds    : used for normalization of the grid data
  Result:     
      Output is given in the form of an array (grid) which holds surface
      values for a square region.  
  _________________________________________________________________________
  """   

  N = int(2**max_level)
  A = numpy.zeros((N,N), dtype = complex)
  random.seed() #Seeds the random number generator
  PI = 3.141592
  for i in range(0,int((N-1)/2)):
    for j in range(0,int((N-1)/2)):
      phase = 2*PI*random.random() #/random.randrange(1,Arand)
      if i != 0 or j != 0:
        rad = pow((i*i + j*j),(-(H+1)/2) )*random.gauss(0.0, sigma)
      else:
        rad = 0.0

      A[i][j] = rad*math.cos(phase) + rad*math.sin(phase)*j 

      if i ==0: 
        i0 = 0
      else:
        i0 = N - i

      if j==0:
        j0 = 0
      else:
        j0 = N - j

      A[i0][j0] = rad * math.cos(phase) - rad*math.sin(phase)*j

  for i in range(1,int((N-1)/2)):
    for j in range(1,int((N-1)/2)):
      phase = 2*PI*random.random() #/random.randrange(1,Arand)
      rad = pow((i*i + j*j),(-(H+1)/2) )*random.gauss(0.0, sigma)
      A[i][N-j] = rad * math.cos(phase) + rad* math.sin(phase)*j
      A[N-i][j] = rad * math.cos(phase) - rad* math.sin(phase)*j

  Grid = numpy.real(pylab.ifft2(( A ))) #Implements the Discrete Inverse Fast Fourier Transform on a 2D lattice (grid)
  if(normalise):
        Grid += numpy.amin(Grid)*-1 + bounds[0]
        Grid = (Grid/numpy.amax(Grid)) * bounds[1]
  return Grid

Second script: plotting the function and obtaining the fractal image:

from __future__ import division #Avoids the floor of the mathematical result of division if the args are ints or longs
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mc
import SSFM2D 


#Parameter assignments
max_level=8               #This is the exponent controlling the grid size. In this case N=2^8=256. Use only integers.
N=2**max_level
sigma=1                   #Variation for random Gauss generation (standardised normal distribution)
H=0.8                     #Hurst exponent (0.8 is a recommended value for natural phenomena)
seed=122                    #Setting the seed for random Gauss generation
print ('The lattice size is '+str(N)+'x'+str(N))


#Lattice initialization
Lattice=np.zeros((256,256))


#Calling Spectral fBm function
Lattice=SSFM2D.SpectralSynthesisFM2D(max_level, sigma, H, seed, normalise=True, bounds=[0,1])


#Plotting the lattice
M=np.zeros((101,101))     #Lattice to be plotted. Size(M) != size(Lattice) to minimize visual fBm periodic features. 
for i in range(0,101):
    for j in range(0,101):
        M[i][j]=Lattice[i][j]


#Lattice statistics              
print M[-101:,-101:].sum()
print M[-101:,-101:].max()
print M[-101:,-101:].min()
print M[-101:,-101:].mean()        


#Writing X, Y and values to a .csv file from scratch
import numpy
import csv
with open('C:\\Users\\Francesco\\Desktop\\Python_files\\csv\\fBm_256x256.csv', 'w') as f: # Change directory if necessary
    writer = csv.writer(f)
    writer.writerow(['X', 'Y', 'Value'])
    for (x, y), val in numpy.ndenumerate(M):
        writer.writerow([x, y, val])

# Plotting - Showing interpolation of randomization
plt.imshow(M[-101:,-101:].T, origin='lower',interpolation='nearest',cmap='Blues', norm=mc.Normalize(vmin=0,vmax=M.max()))
title_string=('fBm: Inverse FFT on Spectral Synthesis')
subtitle_string=('Lattice size: 101x101')
plt.suptitle(title_string, y=0.99, fontsize=17)
plt.title(subtitle_string, fontsize=9)
plt.show()

# Makes a custom list of tick mark intervals for colour bar (assumes minimum is always zero)
numberOfTicks = 5
ticksListIncrement = M.max()/(numberOfTicks)
ticksList = []
for i in range((numberOfTicks+1)):
    ticksList.append(ticksListIncrement * i) 

cb=plt.colorbar(orientation='horizontal', format='%0.2f', ticks=ticksList) 
cb.set_label('Values') 
plt.show()
plt.xlim(0, 100)
plt.xlabel('Easting (Cells)') 
plt.ylim(100, 0)
plt.ylabel('Northing (Cells)')