Python toolboxes for state-space estimation via subspace estimation
You can use pyvib to do frequency based subspace identification. Beware that there is no estimation of the initial state. It is possible to do optimization of the identified model, if the data is not perfectly linear. See the implemenentation, maybe you can use it, in case you want to do your own implementation.

Somewhat incomplete example. Take a look at the example to get a working code:

git clone --depth 1 https://github.com/pawsen/pyvib.git
export PYTHONPATH=pyvib

from pyvib.subspace import Subspace
from pyvib.signal import Signal
# partion the data so the format is (npp,inputs/outputs,R,P)
sig = Signal(uest,yest,fs=fs)
sig.lines = lines  # which freq lines should be used
sig.bla()  # computer best linear approx.
nvec = [2,3]  # model size to scan over
maxr = 5  # max number of rows
model = Subspace(sig)
models, infodict = model.scan(nvec, maxr)
errvec = model.extract_model(yval, uval)  # extract best model on fresh data
model.estimate(n, r)  # or do direct estimation if you know sys size


I have written the code. You are more than welcome to criticize and suggest improvements.