how to replace N(0,1) innovations with fractional gaussian noise ?

hi everyone,

my question is, how can i use fractional gaussian noise in place of standard gaussian noise to describe the noise in a mean reverting process of the form X(t+1)=phi_i*X(t)+c_i+sigma_i*X(t)*N(0,1). My problem is that i can replicate N(0,1) with the function normpdf(Data,mean,sigma) but the only matlab function that gives gaussian noise is "wfbm" ,that does not take my data as input. Please , anyone who is familiar with these notions to help me!!!!