This is a good question.

The steps to compute Zwicker's loudness are

1. scale the sound segment, so that the signal

represents the corresponding SPL

2. call the lautheitII function implemented in the SI++

system

3. transform the sone envelope to phon

4. calculate the pth percentile of the distribution of

loudness levels

These steps are part of an algorithmic approach measuring Zwicker's loudness in

https://www.tcelectronic.com/media/1014016/skovenborg_2004_loudness_m.pdf

In the same paper, the steps are given in an ISO method to compute Zwicker's loadness: 

1. apply calibration gain, so that the digital

representation of a 1 kHz test tone corresponds to

the appropriate loudness level in phon

2. third-octave filterbank

3. calculate the Leq (in dB) per frequency band

4. perform the three-step procedure specified in ISO

532B, using the C implementation

5. transform the resulting total loudness in sone to

phon

See, also, from RG:

https://www.researchgate.net/publication/233523623_A_Revision_of_Zwicker's_Loudness_Model

https://www.researchgate.net/publication/220692458_Assessment_and_Prediction_of_Speech_Quality_in_Telecommunications

More to the point, see

https://www.researchgate.net/publication/263649325_A_Study_on_BSR_Noise_and_Sound_Quality_Property_for_Vehicle_Interior_Module

Article A Revision of Zwicker's Loudness Model

Book Assessment and Prediction of Speech Quality in Telecommunications

Article A Study on BSR Noise and Sound Quality Property for Vehicle ...

Thanks Mr.James F Peters for your suggestions. I will try to implement it.

Good question and good answer :)

Perhaps a little late, but there is a toolbox available which contains matlab implementations of a number of loudness models, which I've found useful in the past. http://genesis-acoustics.com/en/loudness_online-32.html