You could attempt this in frequency domain.
G = frd(OriginalModel,w) % w is a suitably chosen frequency grid
ReducedModel = ssest(G, 4) % reduction (in L2 error since) to 4th order
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!