Simulating Gaussian Fractional Noise

I am trying to generate Gaussian Fractional Noise using the following code :

#autocovariance of Gaussian Noise
import numpy as np
from scipy.linalg import circulant,eigvals,eig,inv
import matplotlib.pyplot as plt

def autocovG(H,n):
    empty_mt = np.zeros(n)
    for m in range(n):
        empty_mt[m] = (np.abs((m-1)/252)**(2*H) - 2 * np.abs(m/252)**(2*H) + np.abs((m+1)/252)**(2*H))*0.5
    return empty_mt

H = 0.10 #Hurst Parameter
first_columns = autocovG(H,1024)
C = circulant(first_columns)
eigen = np.diag(np.fft.fft(first_columns))
eigenv = np.array([[(np.exp((-2*np.pi*1j/eigen.shape[0])*j*k))/eigen.shape[0] for j in range(0,eigen.shape[0])] for k in range(0,eigen.shape[0])])
M = np.dot(eigenv,np.dot(np.sqrt(eigen),inv(eigenv)))
G = np.dot(M,np.random.default_rng().normal(0,scale = 1,size=(1024,1))).real
plt.plot(np.cumsum(G))
plt.show()

However, after the eigendecomposition I find the cumulative sum of my fractional Gaussian noise yielding this:

Which is obviously not a fractional Brownian motion, which should be stationary with a hurst parameter of 0.10. What is wrong with the code?

Update
When I recompose my circulant matrix through CVC*, with C* being the reverse eigenvector's matrix, I find that it gives the transpose of my original circulant matrix.