A method of estimating the MFB50 combustion index, namely the crank angle at which 50% of the fuel mass has been burnt in a cylinder of an internal combustion engine provided with a drive shaft coupled to at least one pair of position sensors (5,5*; 5A, 5B), each of which being arranged at a respective end of the drive shaft ; the estimation method comprising the steps of: capturing the signals coming from the two position sensors (5,5*; 5A, 5B); determining the angular torsion (Δθ) of the drive shaft based on the signals coming from the two position sensors (5,5*; 5A, 5B); determining, by means of a frequency analysis of the angular torsion (Δθ) of the drive shaft , a first harmonic (Ω 1 ) of the angular torsion (Δθ) signal characterized by its own module (|Ω 1 |) and by its own phase (Arg(Ω 1 )); and estimating, during normal operation,the MFB50 combustion index of each cylinder based on the phase (Arg(Ω 1 ) of the first harmonic (Ω 1 ) of the angular torsion (Δθ) signal and a first correlation coefficient (P i ) of each cylinder between the MFB50 combustion index of each cylinder and the phase (Arg(Ω 1 ) of the first harmonic (Ω 1 ) of the angular torsion (Δθ) signal; wherein the first correlation coefficient (P i ) is determined for each cylinder , during a preliminary set-up phase, as follows: P i = MFB 50 i / Arg Ω 1 i with i = 1 ... n MFB50 i MFB50 combustion index based on the cylinder ; Arg (Ω 1 ) 1 phase Arg (Ω 1 ) of the first harmonic (Ω 1 ) of the angular torsion (Δθ) signal; P i first correlation coefficient between the MFB50 combustion index and the phase Arg (Ω 1 ) of the first harmonic (Ω 1 ) of the angular torsion (Δθ) signal; and n number of cylinders . capturing the signals coming from the two position sensors (5,5*; 5A, 5B); determining the angular torsion (Δθ) of the drive shaft based on the signals coming from the two position sensors (5,5*; 5A, 5B); determining, by means of a frequency analysis of the angular torsion (Δθ) of the drive shaft , a first harmonic (Ω 1 ) of the angular torsion (Δθ) signal characterized by its own module (|Ω 1 |) and by its own phase (Arg(Ω 1 )); and estimating, during normal operation,the MFB50 combustion index of each cylinder based on the phase (Arg(Ω 1 ) of the first harmonic (Ω 1 ) of the angular torsion (Δθ) signal and a first correlation coefficient (P i ) of each cylinder between the MFB50 combustion index of each cylinder and the phase (Arg(Ω 1 ) of the first harmonic (Ω 1 ) of the angular torsion (Δθ) signal; wherein the first correlation coefficient (P i ) is determined for each cylinder , during a preliminary set-up phase, as follows: P i = MFB 50 i / Arg Ω 1 i with i = 1 ... n MFB50 i MFB50 combustion index based on the cylinder ; Arg (Ω 1 ) 1 phase Arg (Ω 1 ) of the first harmonic (Ω 1 ) of the angular torsion (Δθ) signal; P i first correlation coefficient between the MFB50 combustion index and the phase Arg (Ω 1 ) of the first harmonic (Ω 1 ) of the angular torsion (Δθ) signal; and n number of cylinders . MFB50 i MFB50 combustion index based on the cylinder ; Arg (Ω 1 ) 1 phase Arg (Ω 1 ) of the first harmonic (Ω 1 ) of the angular torsion (Δθ) signal; P i first correlation coefficient between the MFB50 combustion index and the phase Arg (Ω 1 ) of the first harmonic (Ω 1 ) of the angular torsion (Δθ) signal; and n number of cylinders . <EOS>