笔者学识有限。 本博客旨在对mitre_sfr核心代码进行简单注解， 详细内容读者仁者见仁智者见智。

SFR算法对ESF曲线进行差分得到LSF；

给LSF加上汉明窗， 降低干扰；

对加窗后的LSF曲线进行DFT， 得到SFR。

/*****************************************************************************/

/* This has been modified from Annex A, to more closely match Annex D and

reduce finite difference errors. Now allows either [-1 1] derivative

(when separation = 0) or [-1/2 0 1/2] derivative (when separation=1)

Inputs: len length of ESF array

AveEdge array of ESF values

separation type of derivative

0 = [-1 1]

1 = [-1/2 0 1/2]

Outputs: AveTmp array of original ESF values

AveEdge array of derivative (LSF) values

centroid centroid of the derivative

*/

// 计算ESF的差分图像LSF， 并计算LSF的质心

void calculate_derivative(unsigned int len, double *AveTmp, double *AveEdge,

double *centroid,int separation)

{

unsigned long i;

double dt, dt1;

dt = 0.0;

dt1 = 0.0;

for (i=0; i< len; i++)

AveTmp[i] = AveEdge[i];

for (i=1; i< len-separation; i++) {

/* Not wasting time with division by 2 since constant factors don't

change SFR computation */

AveEdge[i] = (AveTmp[i+separation]-AveTmp[i-1]);

if (separation == 1)

AveEdge[i] /= 2.0;

dt += AveEdge[i] * (double)i;

dt1 += AveEdge[i];

}

*centroid = dt/dt1;

AveEdge[0] = AveEdge[1];

if (separation == 1) AveEdge[len-1] = AveEdge[len-2];

}

/*****************************************************************************/

void apply_hamming_window( unsigned short alpha,

unsigned int oldlen,

unsigned short newxwidth,

double *AveEdge, long *pcnt2)

{

long i,j,k, begin, end, edge_offset;

double sfrc;

/* Shift the AvgEdge[i] vector to centre the lsf in the transform window */

// 将LSF的最大值移到中心位置， 两边由于移动造成的空位由0补齐

edge_offset = (*pcnt2) - (oldlen/2);

if (edge_offset != 0) {

if (edge_offset < 0 ) {

for (i=oldlen-1; i > -edge_offset-1; i--)

AveEdge[i] = (AveEdge[i+edge_offset]);

for (i=0; i < -edge_offset; i++)

AveEdge[i] = 0.00; /* last operation */

} else {

for (i=0; i < oldlen-edge_offset; i++)

AveEdge[i] = (AveEdge[i+edge_offset]);

for (i=oldlen-edge_offset; i < oldlen; i++)

AveEdge[i] = 0.00;

}

}

/* Multiply the LSF data by a Hamming window of width NEWXWIDTH*alpha */

// 在当前上下文中， 下面这段代码应该不会起任何作用， 不知道为什么会有这么复杂而冗余的逻辑

begin = (oldlen/2)-(newxwidth*alpha/2);

if (begin < 0) begin = 0;

end = (oldlen/2)+(newxwidth*alpha/2);

if (end > oldlen ) end = oldlen;

for (i=0; i< begin; i++)

AveEdge[i] = 0.0;

for (i=end; i< oldlen; i++)

AveEdge[i] = 0.0;

// 给begin和end之间的数据加上汉明窗

// 汉明窗 W(n,α ) = (1 -α ) - α cos(2*PI*n/(N-1)) ，(0≤n≤N-1)

// 一般情况下，α取0.46

// 下面计算方法等于窗发生了平移(故符号发生了变化)， 结果是一样的

for (i=begin,j = -newxwidth*alpha/2; i < end; i++,j++) {

sfrc = 0.54 + 0.46*cos( (M_PI*(double)j)/(newxwidth*alpha/2) );

AveEdge[i] = (AveEdge[i])*sfrc;

}

// 将有效数据平移到起始位置

// 在当前上下文中， 下面这段代码应该不会起任何作用， 不知道为什么会有这么复杂而冗余的逻辑

if (begin != 0) /* Shift LSF to begin at index 0 (rather than begin) */

for (k=0, i=begin; k

/*****************************************************************************/

/* This is the DFT magnitude code */

unsigned short ftwos(long number, double dx, double *lsf,

long ns, double ds, double *sfr)

{

double a, b, twopi, g;

long i,j;

// n-1 k

// DFT ==> X[k] = Σ x[n]e^(-j2π - n)

// n=0 N

twopi = 2.0 * M_PI;

for (j = 0; j < ns; j++){

g = twopi * dx * ds * (double)j;

for (i = 0, a = 0, b = 0; i < number; i++) {

a += lsf[i] * cos(g * (double)(i));

b += lsf[i] * sin(g * (double)(i));

}

sfr[j] = sqrt(a * a + b * b);

}

return 0;

}