write a python code for this:

The steps involved in the energy sorted matrix pencil method are as follows [34]–[36]:
(1) Obtain free vibration response measurements with N samples [37]:
y ¨(kΔt)=x ¨(kΔt)≈∑_(j=1)^M R_j z_j^k, for k = 0, 1, 2, ..., N-1	(18)

where Δt is the sample period, R_j is the residues or complex amplitudes, and z_j=e^(λ_j Δt).
(2) Arrange a single response measurement into the Hankel matrices Y1 and Y2 of dimension (N-L) × L as follows:
[Y_1 ]=[■(x ̈(0)&x ̈(1)&x ̈(2)&⋯&x ̈(L-1)@x ̈(1)&x ̈(2)&x ̈(3)&⋯&x ̈(L)@⋮&⋮&⋮&⋱&⋮@x ̈(N-L-2)&x ̈(N-L-1)&x ̈(N-L)&⋯&x ̈(N-3)@x ̈(N-L-1)&x ̈(N-L)&x ̈(N-L+1)&⋯&x ̈(N-2))]
and,
[Y_2 ]=[■(x ̈(1)&x ̈(2)&x ̈(3)&⋯&x ̈(L)@x ̈(2)&x ̈(3)&x ̈(4)&⋯&x ̈(L+1)@⋮&⋮&⋮&⋱&⋮@x ̈(N-L-1)&x ̈(N-L)&x ̈(N-L+1)&⋯&x ̈(N-2)@x ̈(N-L)&x ̈(N-L+1)&x ̈(N-L+2)&⋯&x ̈(N-1))]	(19)

The parameter L is the pencil parameter which is chosen as N/2 for efficient noise filtering.
(3) Form the combined data matrix Y with dimension of (N-L) × (L+1). This augments Y1 and Y2 into a single matrix. Apply SVD to matrix Y as follows:
[Y]=[■(y ̈(0)&y ̈(1)&y ̈(2)&⋯&y ̈(L)@y ̈(1)&y ̈(2)&y ̈(3)&⋯&y ̈(L+1)@⋮&⋮&⋮&⋱&⋮@y ̈(N-L-2)&y ̈(N-L-1)&y ̈(N-L)&⋯&y ̈(N-2)@y ̈(N-L-1)&y ̈(N-L)&y ̈(N+1)&⋯&y ̈(N-1))]
and,
[Y]=[U][Σ][V]^H	(20)

where the superscript H denotes the conjugate transpose. Notice that [Y1] is obtained from [Y] by deleting the last column, and [Y2] is obtained from [Y] be deleting the first column.
(4) Construct the reduced matrix V' containing only the M dominant right singular vectors of V.  For reducing noise from the noisy data, the M dominant singular values of Σ are selected as 〖σ_d≥e〗^(-p) σ_max, where σ_d is the dominant singular value, σ_max is the maximum singular value, and p is the filtering factor. The singular values, which are below e^(-p) σ_max, are recognized as noise singular value and are removed from the data recovery algorithm [38].
(5) Use V' to form the filtered Hankel matrices Y1 and Y2. This filters the noise in the data:
[Y_1 ]=[U][Σ^' ] [V_1^' ]^H=[V_1^' ]^H [V_1^' ]
[Y_2 ]=[U][Σ^' ] [V_2^' ]^H=[V_2^' ]^H [V_1^' ]	(21)

The matrix V_1^' is derived from the matrix V' by eliminating its last row, the matrix V_2^' is obtained from V' by removing its first row, and [Σ^' ] is obtained from the M columns of [Σ] corresponding to the M dominant singular values.
(6) Solve the generalized eigenvalue problem ([Y_1 ]^† [Y_2 ]x=λx) to extract eigenvalues or poles λ=z_i, where using Moore-Penrose pseudoinverse [Y_1 ]^†=([Y_1 ]^H [Y_1 ])^(-1) [Y_1 ]^H. 
(7) Compute the natural circular frequencies ω_i and damping ratios ξ_i as follows:
z_i=e^(λ_i Δt)⟹λ_i=(ln⁡(z_i ))/Δt
ω_i=√((real(λ_i ))^2+(imag(λ_i ))^2 )
ξ_i=abs(real(λ_i ))/√((real(λ_i ))^2+(imag(λ_i ))^2 )	(22)

The poles are also used to calculate the mode shapes (R_j) as follows:
[■(z_1^0&z_2^0&⋯&z_M^0@z_1^1&z_2^1&⋯&z_M^1@⋮&⋮&⋱&⋮@z_1^(N-1)&z_2^(N-1)&⋯&z_M^(N-1) )][■(R_1@R_2@⋮@R_M )]=[■(y ̈(0)@y ̈(1)@⋮@y ̈(N-1))]	(23)

The modal shapes can be determined by least-square method.
(8) Calculate the energy Ei for each mode (each zi):
E_i=∑_(j=1)^N▒(real(R_i z_i^j ))^2 , for i = 1, 2, 3, …, M – 1, M	(24)

(9) Sort the modes based on their energy content Ei, with highest energy mode listed first.
(10) The modes with highest energy content are the dominant modes in the system.