Given the input and output of a system as a discrete function of time, I need to find the spring constant and damping coefficient of an equivalent system. I can not use tfest but as I don't know the number of poles. Is there any other way around this

I have an input u(t) which is a discrete function with a sampling time of .002 seconds. Similarly I have an output of the system of y(t) with the same sampling time of .002 seconds. So what I did was to use spectral analysis of these time domain data using the "spa" function and found the frequency response . I have plotted the frequency response using the bode function. I calculated the 1st peak of magnitude graph and approximated it to natural frequency. Since I knew mass, I found out spring constant and with the help of transfer function obtained through spa function, I found damping coefficient as a function of frequency.

But calculating natural frequency from the bode plot is very inaccurate and upon that the approximation adds to the error. Is there any way other way to solve this?