1. Create PMREM(Prefiltered Mipmaped Radiance Environment map).
a) Get the hdr-map, as an example pisa.hdr(0) (bunch of here (1)).

b) Download HDR Shop (2) free edition.
c) Open HDR Shop.

• open your .hdr – go to Image -> Panorama -> Panoramic Transform

• settings Source : Longitude,  Dest : Cube Env(Vertical Cross) – Convert

• Save as .. *.HDR -> pisa_cross.hdr

d) Download ModifiedCubeMapGen-1_66 (3), many thanks to Sebastien Lagarde (4).

• Load Cube Cross (choose pisa_cross.hdr)
• Click on Filter Cubemap(to recieve better IBL resutl try to change settings)

•  Click CHECK on Save MipChain and SaveCubeMap(.dds) – pisa_cross.dds

2. Generate Sh-Coeff.
a) Upload to engine our PMREM(pisa_cross.dds)
b) Calculate coefficients of spherical harmonics from cube texture

Example with functions:restore lighting value from sh-coeff, for direction,  add and scale.

void sphericalHarmonicsEvaluateDirection(float * result, int order,
const Math::Vector3 & dir)
{
result[0] = 0.282095;
result[1] = 0.488603 * dir.y;
result[2] = 0.488603 * dir.z;
result[3] = 0.488603 * dir.x;
result[4] = 1.092548 * dir.x*dir.y;
result[5] = 1.092548 * dir.y*dir.z;
result[6] = 0.315392 * (3.f*dir.z*dir.z - 1.f);
result[7] = 1.092548 * dir.x * dir.z;
result[8] = 0.546274 * (dir.x*dir.x - dir.y*dir.y);
}

void sphericalHarmonicsAdd(float * result, int order,
const float * inputA, const float * inputB)
{
const int numCoeff = order * order;
for (int i = 0; i < numCoeff; i++)
{
result[i] = inputA[i] + inputB[i];
}
}

void sphericalHarmonicsScale(float * result, int order,
const float * input, float scale)
{
const int numCoeff = order * order;
for (int i = 0; i < numCoeff; i++)
{
result[i] = input[i] * scale;
}
}

void sphericalHarmonicsFromTexture(GLuint cubeTexture,
std::vector<Math::Vector3> & output, const uint order)
{
const uint sqOrder = order*order;

// allocate memory for calculations
output.resize(sqOrder);
std::vector<float> resultR(sqOrder);
std::vector<float> resultG(sqOrder);
std::vector<float> resultB(sqOrder);

// variables that describe current face of cube texture
GLubyte* data;
GLint width, height;
GLint internalFormat;
GLint numComponents;

// initialize values
float fWt = 0.0f;
for (uint i = 0; i < sqOrder; i++)
{
output[i].x = 0;
output[i].y = 0;
output[i].z = 0;
resultR[i] = 0;
resultG[i] = 0;
resultB[i] = 0;
}
std::vector<float> shBuff(sqOrder);
std::vector<float> shBuffB(sqOrder);

const GLenum cubeSides[6] = {
GL_TEXTURE_CUBE_MAP_POSITIVE_Y, // Top
GL_TEXTURE_CUBE_MAP_NEGATIVE_X, // Left
GL_TEXTURE_CUBE_MAP_POSITIVE_Z, // Front
GL_TEXTURE_CUBE_MAP_POSITIVE_X, // Right
GL_TEXTURE_CUBE_MAP_NEGATIVE_Z, // Back
GL_TEXTURE_CUBE_MAP_NEGATIVE_Y // Bottom
};

// bind current texture
glBindTexture(GL_TEXTURE_CUBE_MAP, cubeTexture);

int level = 0;
// for each face of cube texture
for (int face = 0; face < 6; face++)
{
// get width and height
glGetTexLevelParameteriv(cubeSides[face], level, GL_TEXTURE_WIDTH, &width);
glGetTexLevelParameteriv(cubeSides[face], level, GL_TEXTURE_HEIGHT, &height);

if (width != height)
{
return;
}

// get format of data in texture

glGetTexLevelParameteriv(cubeSides[face], level,
GL_TEXTURE_INTERNAL_FORMAT, &internalFormat);

// get data from texture
if (internalFormat == GL_RGBA)
{
numComponents = 4;
data = new GLubyte[numComponents * width * width];
}
else if (internalFormat == GL_RGB)
{
numComponents = 3;
data = new GLubyte[numComponents * width * width];
}
else
{
return;
}
glGetTexImage(cubeSides[face], level, internalFormat, GL_UNSIGNED_BYTE, data);

// step between two texels for range [0, 1]
float invWidth = 1.0f / float(width);
// initial negative bound for range [-1, 1]
float negativeBound = -1.0f + invWidth;
// step between two texels for range [-1, 1]
float invWidthBy2 = 2.0f / float(width);

for (int y = 0; y < width; y++)
{
// texture coordinate V in range [-1 to 1]
const float fV = negativeBound + float(y) * invWidthBy2;

for (int x = 0; x < width; x++)
{
// texture coordinate U in range [-1 to 1]
const float fU = negativeBound + float(x) * invWidthBy2;

// determine direction from center of cube texture to current texel
Math::Vector3 dir;
switch (cubeSides[face])
{
case GL_TEXTURE_CUBE_MAP_POSITIVE_X:
dir.x = 1.0f;
dir.y = 1.0f - (invWidthBy2 * float(y) + invWidth);
dir.z = 1.0f - (invWidthBy2 * float(x) + invWidth);
//dir = -dir;
break;
case GL_TEXTURE_CUBE_MAP_NEGATIVE_X:
dir.x = -1.0f;
dir.y = 1.0f - (invWidthBy2 * float(y) + invWidth);
dir.z = -1.0f + (invWidthBy2 * float(x) + invWidth);
//dir = dir;
break;
case GL_TEXTURE_CUBE_MAP_POSITIVE_Y:
dir.x = -1.0f + (invWidthBy2 * float(x) + invWidth);
dir.y = 1.0f;
dir.z = -1.0f + (invWidthBy2 * float(y) + invWidth);
//dir = dir;
break;
case GL_TEXTURE_CUBE_MAP_NEGATIVE_Y:
dir.x = -1.0f + (invWidthBy2 * float(x) + invWidth);
dir.y = -1.0f;
dir.z = 1.0f - (invWidthBy2 * float(y) + invWidth);
//dir = dir; //!
break;
case GL_TEXTURE_CUBE_MAP_POSITIVE_Z:
dir.x = -1.0f + (invWidthBy2 * float(x) + invWidth);
dir.y = 1.0f - (invWidthBy2 * float(y) + invWidth);
dir.z = 1.0f;
break;
case GL_TEXTURE_CUBE_MAP_NEGATIVE_Z:
dir.x = 1.0f - (invWidthBy2 * float(x) + invWidth);
dir.y = 1.0f - (invWidthBy2 * float(y) + invWidth);
dir.z = -1.0f;
break;
default:
return;
}

// normalize direction
dir = Math::Normalize(dir);
// scale factor depending on distance from center of the face
const float fDiffSolid = 4.0f / ((1.0f + fU*fU + fV*fV) *
sqrtf(1.0f + fU*fU + fV*fV));
fWt += fDiffSolid;

// calculate coefficients of spherical harmonics for current direction
sphericalHarmonicsEvaluateDirection(shBuff.data(), order, dir);

// index of texel in texture
uint pixOffsetIndex = (x + y * width) * numComponents;
// get color from texture and map to range [0, 1]
Math::Vector3 clr(
float(data[pixOffsetIndex]) / 255,
float(data[pixOffsetIndex + 1]) / 255,
float(data[pixOffsetIndex + 2]) / 255
);

//clr.x = pow(clr.x, 1.1f);
//clr.y = pow(clr.y, 1.1f);
//clr.z = pow(clr.z, 1.1f);

clr.x = clr.x*0.5f;
clr.y = clr.y*0.5f;
clr.z = clr.z*0.5f;

// scale color and add to previously accumulated coefficients
sphericalHarmonicsScale(shBuffB.data(), order,
shBuff.data(), clr.x * fDiffSolid);
sphericalHarmonicsAdd(resultR.data(), order,
resultR.data(), shBuffB.data());

sphericalHarmonicsScale(shBuffB.data(), order,
shBuff.data(), clr.y * fDiffSolid);
sphericalHarmonicsAdd(resultG.data(), order,
resultG.data(), shBuffB.data());

sphericalHarmonicsScale(shBuffB.data(), order,
shBuff.data(), clr.z * fDiffSolid);
sphericalHarmonicsAdd(resultB.data(), order,
resultB.data(), shBuffB.data());
}
}

delete[] data;
}

// final scale for coefficients
const float fNormProj = (4.0f * Math::PI<float>()) / fWt;
sphericalHarmonicsScale(resultR.data(), order, resultR.data(), fNormProj);
sphericalHarmonicsScale(resultG.data(), order, resultG.data(), fNormProj);
sphericalHarmonicsScale(resultB.data(), order, resultB.data(), fNormProj);

// save result
for (uint i = 0; i < sqOrder; i++)
{
output[i].x = resultR[i];
output[i].y = resultG[i];
output[i].z = resultB[i];
}

//glBindTexture(GL_TEXTURE_CUBE_MAP, 0);
}

Thanks to Reev(5) from gamedev.ru

Example in-game SH-coeff shader:

const float PhongRandsConsts[32] =
{
0,
188,
137,
225,
99,
207,
165,
241,
71,
198,
151,
233,
120,
216,
177,
248,
50,
193,
145,
229,
110,
212,
171,
244,
86,
203,
158,
237,
129,
220,
182,
252
};

inline float tosRGBFloat(float rgba)
{
float srgb = (rgba*rgba)*(rgba*0.2848f + 0.7152f);
return srgb;
}

Math::Vector4* GetPhongRands()
{
static Math::Vector4 rands[32];
float r1 = 0.032f;

for (int it = 0, end = 32; it != end; ++it)
{
float r2 = tosRGBFloat(PhongRandsConsts[it] / 255.f);

//float r2 = Platform::Random::RandFloat(0.f, 1.0f);

rands[it].x = r1;
rands[it].y = r2;
rands[it].z = Math::Cos(2.f * Math::PI<float>() * r2);
rands[it].w = Math::Sin(2.f * Math::PI<float>() * r2);

r1 += 0.032f;
}

return rands;
}