I would use the Signal Processing Toolbox invfreqz function. You will have to calculate the transfer functions between two signals at a time. I am not aware of any function that can identify a multiple-input-multiple-output (MIMO) system.
The System Identification Toolbox has a number of other functions for system identification.
