I have 40 samples that contain information about 6 variables (hence a 40x6 data matrix). Each variable (column) has a cumulative distribution function (marginal distribution) based on the 40 values, which is not Gaussian, whereas correlation matrix exist that relates the variable among each other (I am using Spearsman rho ). What I would like to do, is to generate 1000 samples for each variable (a 1000x6 data matrix) that satisfies the correlation values of the original sampling datasets, and follows the empirical cumulative distribution function of the original 40 samples for each variable. My understanding is that I can use copula for that purpose. So, the steps that I have in mind are:

Calculate the mean and standard deviation for each variable, and generate gaussian marginal distributions for a given correlation (Spearsman rho) matrix. The values of the gaussian variables are not restricted in [0,1].
Next, convert the marginal distribution of each variable into a unit distribution (range of values for each variable lies within [0,1]). This is done without changing the ordering of the data, hence the ranking correlation coefficient remains the same.
Lastly, I am using the cumulative distribution functions of the initial dataset (40x6 data matrix) to generate simulated values for each variable, particularly by quantile-matching (for a specific value of the unit variable, which I assume that equals the cumulative distribution value (vertical axis), I find the simulated value by finding the value at the x-axis that corresponds to this cdf value).
Is the above description accurate? Sorry for the long description, I am new on this and I am not very confident to jump parts of the discussion.I have 40 samples that contain information about 6 variables (hence a 40x6 data matrix). Each variable (column) has a cumulative distribution function (marginal distribution) based on the 40 values, which is not Gaussian, whereas correlation matrix exist that relates the variable among each other (I am using Spearsman rho ). What I would like to do, is to generate 1000 samples for each variable (a 1000x6 data matrix) that satisfies the correlation values of the original sampling datasets, and follows the empirical cumulative distribution function of the original 40 samples for each variable. My understanding is that I can use copula for that purpose. So, the steps that I have in mind are:

Calculate the mean and standard deviation for each variable, and generate gaussian marginal distributions for a given correlation (Spearsman rho) matrix. The values of the gaussian variables are not restricted in [0,1].
Next, convert the marginal distribution of each variable into a unit distribution (range of values for each variable lies within [0,1]). This is done without changing the ordering of the data, hence the ranking correlation coefficient remains the same.
Lastly, I am using the cumulative distribution functions of the initial dataset (40x6 data matrix) to generate simulated values for each variable, particularly by quantile-matching (for a specific value of the unit variable, which I assume that equals the cumulative distribution value (vertical axis), I find the simulated value by finding the value at the x-axis that corresponds to this cdf value).
Is the above description accurate? Sorry for the long description, I am new on this and I am not very confident to jump parts of the discussion.