原创 STM32学习日志(24)----使用dsp库的FFT函数测相位

attachment download

/**
  ******************************************************************************
  * @file FFT_Demo/src/main.c
  * @author  MCD Application Team
  * @version  V2.0.0
  * @date  04/27/2009
  * @brief  Main program body
  ******************************************************************************
  * @copy
  *
  * THE PRESENT FIRMWARE WHICH IS FOR GUIDANCE ONLY AIMS AT PROVIDING CUSTOMERS
  * WITH CODING INFORMATION REGARDING THEIR PRODUCTS IN ORDER FOR THEM TO SAVE
  * TIME. AS A RESULT, STMICROELECTRONICS SHALL NOT BE HELD LIABLE FOR ANY
  * DIRECT, INDIRECT OR CONSEQUENTIAL DAMAGES WITH RESPECT TO ANY CLAIMS ARISING
  * FROM THE CONTENT OF SUCH FIRMWARE AND/OR THE USE MADE BY CUSTOMERS OF THE
  * CODING INFORMATION CONTAINED HEREIN IN CONNECTION WITH THEIR PRODUCTS.
  *
  * <h2><center>&copy; COPYRIGHT 2009 STMicroelectronics</center></h2>
  */
/* Includes ------------------------------------------------------------------*/
#include "stm32f10x.h"
#include <math.h>
#include "stm3210b_lcd.h"
#include "stm32_dsp.h"
#include "table_fft.h"   
//------------------------------------------------------------------------------

#define FreqSample      (50*NPT) // 50Hz*64=3200
#define ScalingFactor   20000  // 幅值系数,最大为32767
#define PERCENT_DC  0.06  // 直流分量对应的百分比,对应直流分量lBUFMAG[0]=1200
#define PI              3.1415926   
#define PI2             6.28318530717959   // 2*3.14
#define NPT             256                /* 采样点数，必须为64,256,1024之一*/

#if (NPT==64)
     #define cr4_fft_stm32 cr4_fft_64_stm32
#elif (NPT==256)
     #define cr4_fft_stm32 cr4_fft_256_stm32
#else
     #define cr4_fft_stm32 cr4_fft_1024_stm32
        error!!!!
        // 使用的ram超过了22K，超过了EK-STM32F仿真板上的STM32F103VBT6的ram
#endif
//------------------------------------------------------------------------------
extern volatile uint32_t TimingDelay ;
long lBUFIN[NPT];           /* Complex input vector */
long lBUFOUT[NPT];          /* Complex output vector */
long lBUFMAG[NPT + NPT/2];  /* Magnitude vector */
float lBUFPHASE[NPT + NPT/2];/* phase vector,单位为度 */

/* Private function prototypes -----------------------------------------------*/
void MyDualSweep(uint32_t freqinc1,uint32_t freqinc2);
void MygSin(long nfill, long Fs, long Freq1, long Freq2, long Ampli);
void powerMag(long nfill, char* strPara);
void power_Phase_Radians(long nfill);
void DSPDemoInit(void);
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int main(void)
{
  DSPDemoInit();

  while (1)
  {
      MyDualSweep(FreqSample/NPT,6*FreqSample/NPT);     // 测试fft
  }
}
//-----------------------------------------------------------------------------

/**
  * @brief  Produces a combination of two sinewaves as input signal
  * @param freq1: frequency increment for 1st sweep，频率步进值1
  *   Freq2: frequency increment for 2nd sweep，频率步进值2
  * @retval : None
  */
void MyDualSweep(uint32_t freqinc1,uint32_t freqinc2)
{
    uint32_t freq;
    
    freq = freqinc1;
    MygSin(NPT, FreqSample, freq, 0, ScalingFactor);       
    cr4_fft_stm32(lBUFOUT, lBUFIN, NPT);
    powerMag(NPT,"");       // lBUFIN[0]=1192 ,lBUFIN[1]=20008
    power_Phase_Radians(NPT);   
         
    freq = freqinc1;
    MygSin(NPT, FreqSample, freq, freqinc2, ScalingFactor);       
    cr4_fft_stm32(lBUFOUT, lBUFIN, NPT);
    powerMag(NPT,"");        // lBUFIN[0]=1192 ,lBUFIN[1]=20008, lBUFIN[6]=4000   
    power_Phase_Radians(NPT);
}
//------------------------------------------------------------------------------
/**
  * @brief  Produces a combination of two sinewaves as input signal
  * @param ill: length of the array holding input signal
  *   Fs: sampling frequency
  *   Freq1: frequency of the 1st sinewave
  *   Freq2: frequency of the 2nd sinewave
  *   Ampli: scaling factor
  * @retval : None
  */
void MygSin(long nfill, long Fs, long Freq1, long Freq2, long Ampli)
{
  uint32_t i;
  float fFs, fFreq1, fFreq2, fAmpli;
  float fZ,fY;

  fFs = (float) Fs;
  fFreq1 = (float) Freq1;
  fFreq2 = (float) Freq2;
  fAmpli = (float) Ampli;

  for (i=0; i < nfill; i++)
  {
    fY =PERCENT_DC/2            // 相位范围0°----360°
            + sin((196.0*PI/180.0)+PI2 * i * (fFreq1/fFs)) ;
    if (fFreq2)
            fY += 0.2*sin((293.0*PI/180.0)+PI2 * i * (fFreq2/fFs));
   
    fZ = fAmpli * fY;
 
// lBUFIN数组中，每个单元数据高字(高16位)中存储采样数据的实部，低字(低16位)存储采样数据的虚部(总是为0)
    lBUFIN= ((short)fZ) << 16 ;  /* sine_cosine  (cos=0x0) */
         // 高16bit为实部，低16bit为虚部，没有的话，就用0代替
  }
}
//------------------------------------------------------------------------------
/**
  * @brief  Removes the aliased part of the spectrum (not tested)
  * @param ill: length of the array holding power mag
  * @retval : None
  */
void onesided(long nfill)
{
  uint32_t i;
 
  lBUFMAG[0] = lBUFMAG[0];
  lBUFMAG[nfill/2] = lBUFMAG[nfill/2];
  for (i=1; i < nfill/2; i++)
  {
    lBUFMAG = lBUFMAG + lBUFMAG[nfill-i];
    lBUFMAG[nfill-i] = 0x0;
  }
}
//------------------------------------------------------------------------------
/**
  * @brief  Compute power magnitude of the FFT transform  进行FFT变换，并计算各次谐波幅值
  * @param ill: length of the array holding power mag
  *   : strPara: if set to "1SIDED", removes aliases part of spectrum (not tested)
  * @retval : None
  */
void powerMag(long nfill, char* strPara)
{
  int32_t lX,lY;
  uint32_t i;