基于傅里叶变换的电力测频算法

参考文献：《基于傅里叶变换的精确频率测量算法》牟龙华，邢锦磊 

 

说明：该函数利用相邻两周波的傅里叶变换，推导出频率差，由此 

计算得到真实频率并根据计算得到的频率调节采样频率进行迭代。

FFT_DOT = 64;

cos_tab[] 和sin_tab[]是1024点的正余弦值表，在查表时注意数组下标的移位处理

/*******************************************************************************
  * *   函数名称: Freq_Calculate
  * *     函数说明: 基于相邻两周波采样数据的傅里叶变换提取频率算法
  * *    输入参数: databuf 采样数据缓冲区指针；pDataU 数据结构体指针
  * *      返回值: 无
  *****************************************************************************/
void Freq_Calculate(int16_t* databuf,DataItem*pDataU)
{
    float deltaF = 0;
    float a,b,c,d;
    float cost = 0;
    uint8_t i = 0,phase1,phase2;
    short sign = 0;
    static float frq = 50.0f;

    for(i=0,a=0,b=0,c=0,d=0;i<FFT_DOT;i++)
    {
        a += databuf[i]*sin_tab[i<<4];
        b += databuf[i]*cos_tab[i<<4];
        c += databuf[i+FFT_DOT]*sin_tab[i<<4];
        d += databuf[i+FFT_DOT]*cos_tab[i<<4];
    }
    
    cost = (a*b + c*d)/(a*d +c*b);   
    if(cost<=1.0f && cost >=-1.0f)deltaF = acosf(cost)*7.957747f;
    else deltaF = 0;
    
    if(a>0 && b>0)phase1 = 1;
    else if(a>0 && b<0)phase1 = 4;
    else if(a<0 && b>0)phase1 = 2;
    else if(a<0 && b<0)phase1 = 3;
    
        if(c>0 && d>0)phase2 = 1;
    else if(c>0 && d<0)phase2 = 4;
    else if(c<0 && d>0)phase2 = 2;
    else if(c<0 && d<0)phase2 = 3;
    
    if(phase1 == phase2)
    {
        if(phase1 == 1 || phase1 == 2)
        {
            if(a>c)sign = 1;
            else sign =-1;
        }
        else if(phase1 == 2 || phase1 == 3)
        {
           if(a>c)sign = -1;
            else sign =1; 
        }
    }
    else
    {
        if(phase2 == (phase1+1)|| phase2 == (phase1-3))sign = 1;
        else if(phase2 == (phase1-1)|| phase2 == (phase1+3))sign = -1;
    }
    
    #if DEBUG
    spDev[_485_I].printf(&spDev[_485_I],"deltaF: %c%.6f\r\n",(sign<0)?'-':'+',deltaF);
    #endif
    
    if(deltaF<12.5f && deltaF > -12.5f && c_RealtimeData.wDataUa > 1000)
    {
        frq = frq + sign*deltaF;
        if(frq>62.5f)frq = 62.5f;
        else if(frq<37.5f)frq = 37.5f;
        
        AnalogData[F].set_data(AnalogData+F,frq*100);
        AnalogData[F].refresh_data(AnalogData+F,0);
        bsp_timer1_periodSet(frq*100);
    }
}