Noise Shape Filter to obtain a given PSD

I have experimental PSD data of disturbances in terms of acceleration on a mass. My first problem is that these data are presented as $ (m/s^2)/\sqrt(Hz)$ which is not the PSD unit. Those unit should be adressed to an amplitude spectral density (ASD) which is directly $ ASD = \sqrt{PSD}$. (I am not sure of that either, beacuase in some books the relations takes into account $\sqrt{2}$. However, my idea is to reconvert into a transfer function (TF) the PSD and then use that TF as a filter. Infact, if the TF is fed by a white noise ($ W(f) = 1$) the output of the filter operation is $ H(f)*1$ i.e. the colored noise in frequency domain $C(f)$. Thus $$ PSD =|H(f)|^2 $$ So I can interpolate the given points and use a least square method to obtain the chosen number of coefficients of the TF $$ min(a_z, b_p) \sum_{i} (|H(2\pi j f)|- \sqrt{PSD}|)^2 $$ ( I don't really know how to do this). My other option is to hit and trial whit some tf Below, there is an image of the data.

To finish I am leaving here what I have done on Matlab

% TRIAL AND ERROR PROCEDURE
% clear all
% close all


% EL ACT data
% points
x = ([0.06 0.1 0.12 0.2 0.24 0.4 1 2 3 10].*1e-3);
y = ([500 150 100 40 30 20 13 12 12 12].*1e-15).^2;     % now is a PSD?
% plot(x,y)           %trial and error to make the interpol suitable
%
loglog(x,y)
hold on
loglog(x,y,'o')


%%
xlog = log10(x);
ylog = log10(y);

% fitting
N = 1000;
pp = polyfit(xlog,ylog,3);
freq_log = linspace(xlog(1),xlog(end),N);
ELEM_sensing_ASD_log = polyval(pp,freq_log);

ELEM_sensing_ASD = 10.^ELEM_sensing_ASD_log;
freq = 10.^freq_log;
loglog(freq,ELEM_sensing_ASD)  %fitted?

%%%%%%%%%%%%
% Data from Ref
s = tf('s');
H_EL = ((s+2e-5)*(s+3e-5)*(s+3e-4))/((s+2e-6)*(s+5e-6)*(s+1e-5));
Num = cell2mat(H_EL.Numerator);
Den = cell2mat(H_EL.Denominator);

[A,B,C,D] = tf2ss(Num,Den);


%% Simulink
T_filter = 100;
% open('noise_shape_trial_SIM')
sim('noise_shape_trial_SIM')

%%

segment_lenth_sample = 1e4;
overlap = 0;
DFT_points = 4e4;
f_sampling = 1000;
[PSD_color,f] = pwelch(color_time); %,segment_lenth_sample,overlap,DFT_points,f_sampling);
loglog(f,(PSD_color))
grid on