[STM32F7] 基于STM32F4/F7的图形应用演示方案

#include "math.h"

#include "fft.h"

//精度0.0001弧度

//复数的交换 

void conjugate_complex(int n,complex in[],complex out[])

{

  int i = 0;

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

  {

    out[i].imag = -in[i].imag;

    out[i].real = in[i].real;

  }        

}

//求所有复数的模

void c_abs(complex f[],float out[],int n)

{

  int i = 0;

  float t;

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

  {

    t = f[i].real * f[i].real + f[i].imag * f[i].imag;

    out[i] = sqrt(t);

  }        

}



//求复数的和 

void c_plus(complex a,complex b,complex *c)

{

  c->real = a.real + b.real;

  c->imag = a.imag + b.imag;

}

//求复数的差  

void c_sub(complex a,complex b,complex *c)

{

  c->real = a.real - b.real;

  c->imag = a.imag - b.imag;        

}

//求复数的积

void c_mul(complex a,complex b,complex *c)

{

  c->real = a.real * b.real - a.imag * b.imag;

  c->imag = a.real * b.imag + a.imag * b.real;        

}

//求复数的商 

void c_div(complex a,complex b,complex *c)

{

  c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag);

  c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag);

}

#define SWAP(a,b)  tempr=(a);(a)=(b);(b)=tempr

void Wn_i(int n,int i,complex *Wn,char flag)

{

  Wn->real = cos(2*PI*i/n);

  if(flag == 1)

  Wn->imag = -sin(2*PI*i/n);

  else if(flag == 0)

  Wn->imag = -sin(2*PI*i/n);

}

//傅里叶变化

void fft(int N,complex f[])

{

  complex t,wn;//中间变量

  int i,j,k,m,n,l,r,M;

  int la,lb,lc;

  /*----计算分解的级数M=log2(N)----*/

  for(i=N,M=1;(i=i/2)!=1;M++); 

  /*----按照倒位序重新排列原信号----*/

  for(i=1,j=N/2;i<=N-2;i++)

  {

    if(i<j)

    {

      t=f[j];

      f[j]=f[i];

      f[i]=t;

    }

    k=N/2;

    while(k<=j)

    {

      j=j-k;

      k=k/2;

    }

    j=j+k;

  }

 

  /*----FFT算法----*/

  for(m=1;m<=M;m++)

  {

    la=pow(2,m); //la=2^m代表第m级每个分组所含节点数                

    lb=la/2;    //lb代表第m级每个分组所含碟形单元数

                 //同时它也表示每个碟形单元上下节点之间的距离

    /*----碟形运算----*/

    for(l=1;l<=lb;l++)

    {

      r=(l-1)*pow(2,M-m);        

      for(n=l-1;n<N-1;n=n+la) //遍历每个分组，分组总数为N/la

      {

        lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号     

        Wn_i(N,r,&wn,1);//wn=Wnr

        c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算

        c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr

        c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr

      }

    }

  }

}

//傅里叶逆变换

void ifft(int N,complex f[])

{

  int i=0;

  conjugate_complex(N,f,f);

  fft(N,f);

  conjugate_complex(N,f,f);

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

  {

    f[i].imag = (f[i].imag)/N;

    f[i].real = (f[i].real)/N;

  }

}

#ifndef __FFT_H__

#define __FFT_H__

 

typedef struct complex //复数类型

{

  float real;                //实部

  float imag;                //虚部

}complex;

 

#define PI 3.1415926535897932384626433832795028841971

///////////////////////////////////////////

void conjugate_complex(int n,complex in[],complex out[]);

void c_plus(complex a,complex b,complex *c);//复数加

void c_mul(complex a,complex b,complex *c) ;//复数乘

void c_sub(complex a,complex b,complex *c);        //复数减法

void c_div(complex a,complex b,complex *c);        //复数除法

void fft(int N,complex f[]);//傅立叶变换 输出也存在数组f中

void ifft(int N,complex f[]); // 傅里叶逆变换

void c_abs(complex f[],float out[],int n);//复数数组取模

////////////////////////////////////////////

#endif

#include "sys.h"

#include "delay.h"

#include "usart.h"

#include "led.h"

#include "math.h"

#include "fft.h"

#define  N    256          //采样点数

#define  Fs   44800        //采样频率

#define  F    175          //分辨率

#define  P1   30           //测试相位

#define  P2   60

#define  P3   90



//FFT测试数据集 输入数组

complex  FFT_256PointIn[N];

//FFT测试数据集 输出数组

float   FFT_256PointOut[N/2];                                                                                

//填入数组                                                                                                                                

void InitBufInArray()

{

 unsigned short i;

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

        {

       FFT_256PointIn[i].real  = 1500 * sin(2*PI * i * 350.0 / Fs+(PI*P1/180)) 

                                             +2700 * sin(2*PI * i * 8400.0 / Fs+(PI*P2/180)) 

                                             +4000 * sin(2*PI * i * 18725.0 / Fs+(PI*P3/180));

                   FFT_256PointIn[i].imag = 0;

    }        

}



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

函数名称:GetPowerMag()

函数功能:计算各次谐波幅值

参数说明:

备　　注:先将FFT_256PointIn分解成实部(X)和虚部(Y)，

         然后计算幅值:(sqrt(X*X+Y*Y)*2/N

         然后计算相位:atan2(Y/X)

作　　者:土耳其冰激凌

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

void GetPowerMag()

{

    unsigned short i;

          float  X,Y,P,Mag;

                 c_abs(FFT_256PointIn,FFT_256PointOut,N/2);

    for(i=0; i<N/2; i++)

    {

                          X = FFT_256PointIn[i].real/N;    //计算实部

                          Y = FFT_256PointIn[i].imag/N;    //计算虚部

                          Mag = FFT_256PointOut[i]*2/N;    //计算幅值

                          P = atan2(Y,X)*180/PI;           //计算相位

                                printf("%d      ",i);

                                printf("%d      ",F*i); 

                                printf("%f      ",Mag);

                          printf("%f      ",P);

                                printf("%f      ",X);

                                printf("%f      \r\n",Y);                        

    }

}

int main(void)

 {        

        delay_init();            //延时函数初始化          

        LED_Init();                          //初始化与LED连接的硬件接口

        uart_init(9600);  //初始化串口 9600波特率 

        printf("这是一个FFT 测试实验\r\n");  

  InitBufInArray(); 

        fft(N,FFT_256PointIn);

        printf("点数   频率  幅值   实部  虚部\n"); 

        GetPowerMag();

        while(1)

                 {

                                LED0=0;

                                LED1=1;

                                delay_ms(300);         //延时300ms

                                LED0=1;

                                LED1=0;

                                delay_ms(300);        //延时300ms

                        }

 }