C语言使用FFT得到信号幅度谱

最新推荐文章于 2024-04-07 23:54:57 发布

*pan

最新推荐文章于 2024-04-07 23:54:57 发布

阅读量8.5k

点赞数 7

分类专栏： C语言学习文章标签： C语言 FFT DFT与FFT计算量比较幅度谱

本文链接：https://blog.csdn.net/qq_33523925/article/details/91048758

版权

C语言学习专栏收录该内容

28 篇文章 10 订阅

订阅专栏

前面的一篇文章我们介绍了使用DFT得到信号的幅度谱的方法，现在我们来看一下FFT实现信号幅度谱。这里我们使用的FFT源程序是徐士良老师的C语言算法程序——快速傅里叶变换。另外，本文也会对DFT以及FFT实现信号幅度谱计算时，所需要的计算量以及花费的时间进行比较。

1.FFT原理

FFT（Fast Fourier Transformation）是离散傅氏变换（DFT）的快速算法，也就是快速傅里叶变换。

我们知道，DFT的计算式为：
在这里插入图片描述
可以看出，DFT需要计算大约N ^2次乘法和N ^2次加法。当N较大时，这个计算量是很大的。我们对DFT进行分解，分解的大致思路是：

设一个多项式为A(x)：
在这里插入图片描述
按下标的奇偶分为两部分：

再设两个式子A1(x)和A2(x)：

那么：

当N为2^m，可以看到：

它们两个式子的计算结果只有后半部分的符号不同，在计算机上实现时，我们计算得到前面一半的结果，后面的结果基本就可以得出来了。并且，如果N=2^m，还可以一直分下去，分到最后都是2点离散傅里叶变换计算。

典型的8点FFT运算图如下：
在这里插入图片描述

2.FFT实现信号幅度谱

前面介绍了FFT的实现基本思路，实际上FFT计算时就是将前面的过程倒过来就可以，也就是说，计算分解后的各个算式，然后再进行相加、相乘，就得到了最后的结果。下面我们以徐士良老师的C语言程序来实现FFT计算。

我们需要分析的信号为：幅度为0.6，频率分别为50Hz与500Hz，采样率为8000的正弦波信号的相加。由前面的分析我们可以知道，n=2^m次的时候，可以取点数为8192，也就是2 ^13。

下面的程序中，满足n=2^k，pr表示采样数据的实部，pi表示采样数据的虚部，fr表示离散傅里叶变换结果的实部，fi表示离散傅里叶变换结果的虚部。

程序代码

  #include "math.h"
  #include "stdio.h"
  
  #define PI 3.1415926535	
      
  void kfft(pr,pi,n,k,fr,fi,l,il)
  int n,k,l,il;
  double pr[],pi[],fr[],fi[];
  { 
	int it,m,is,i,j,nv,l0;
    double p,q,s,vr,vi,poddr,poddi;
    for (it=0; it<=n-1; it++)
      { m=it; is=0;
        for (i=0; i<=k-1; i++)
          { j=m/2; is=2*is+(m-2*j); m=j;}
        fr[it]=pr[is]; fi[it]=pi[is];
      }
    pr[0]=1.0; pi[0]=0.0;
    p=6.283185306/(1.0*n);
    pr[1]=cos(p); pi[1]=-sin(p);
    if (l!=0) pi[1]=-pi[1];
    for (i=2; i<=n-1; i++)
      { p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];
        s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
        pr[i]=p-q; pi[i]=s-p-q;
      }
    for (it=0; it<=n-2; it=it+2)
      { vr=fr[it]; vi=fi[it];
        fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];
        fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];
      }
    m=n/2; nv=2;
    for (l0=k-2; l0>=0; l0--)
      { m=m/2; nv=2*nv;
        for (it=0; it<=(m-1)*nv; it=it+nv)
          for (j=0; j<=(nv/2)-1; j++)
            { p=pr[m*j]*fr[it+j+nv/2];
              q=pi[m*j]*fi[it+j+nv/2];
              s=pr[m*j]+pi[m*j];
              s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
              poddr=p-q; poddi=s-p-q;
              fr[it+j+nv/2]=fr[it+j]-poddr;
              fi[it+j+nv/2]=fi[it+j]-poddi;
              fr[it+j]=fr[it+j]+poddr;
              fi[it+j]=fi[it+j]+poddi;
            }
      }
    if (l!=0)
      for (i=0; i<=n-1; i++)
        { fr[i]=fr[i]/(1.0*n);
          fi[i]=fi[i]/(1.0*n);
        }
    if (il!=0)
      for (i=0; i<=n-1; i++)
        { pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
          if (fabs(fr[i])<0.000001*fabs(fi[i]))
            { if ((fi[i]*fr[i])>0) pi[i]=90.0;
              else pi[i]=-90.0;
            }
          else
            pi[i]=atan(fi[i]/fr[i])*360.0/6.283185306;
        }
    return;
  }
 
  int main()
  { 
	  int i;
	  double t;
	  double pr[8192],pi[8192],fr[8192],fi[8192];
      
      //产生抽样序列
	  for (i=0; i<8192; i++)
      { 
		 t=i/8192.0;
		 pr[i]=0.6*(sin(2*PI*500*t)+sin(2*PI*50*t));
		 pi[i]=0.0;
	  }
     
	  kfft(pr,pi,8192,13,fr,fi,0,1);  //调用fft函数
 
      for (i=0; i<8192; i++)
      {
		  t=i/8192.0;
		  printf("%.6f %.16f\n",t,pr[i]); //输出数据
	  }

      return 0;
  }

作图命令

  cd path
  tcc fft.c -o fft.exe
  fft.exe>fft.dat
  gnuplot
  plot "fft.dat" w l

运行结果
在这里插入图片描述
可以得到，DFT和FFT最终的计算结果是一样的。

3.DFT与FFT的计算量比较

前面已经了解过，FFT是DFT改进后的快速计算方法，N点DFT共需要N^ 2次复数乘法和N(N-1)次复数加法，共4N^2次实数乘法和(2N ^2+2N*(N-1))次实数加法。当N很大时，这是一个非常大的计算量。

利用FFT算法之后，任何一个N为2的整数幂（即N= 2 ^M）的DFT，都可以通过M次分解，最后成为2点的DFT来计算。M次分解构成了从x(n)到X(k)的M级迭代计算，每级由N/2个蝶形运算组成。完成一个蝶形计算需一次乘法和两次复数加法。因此，完成N点的时间抽选FFT计算的总运算量为：

复数乘法次数：
在这里插入图片描述
复数加法次数：

接下来我们从程序的运行时间来看一下两种算法的计算情况，方法是分别计算两种算法进行计算时所需的时间（我这里计算的是进行DFT或FFT计算并输出结果的时间），具体实现如下。

DFT算法时间计算

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>

#define PI 3.1415926535	

struct data
{
	double real,imag;
};

int main()
{
	int N=8000;        //序列长度

	double Squence[8000];       //序列
	double Amplitude[8000];     //幅值

	struct data one[8000];       //单个数据
	struct data sum;               //求和的结果

	int n;        
	int k;        //表示第k个数

	double t,f;   //信号的时间 频率

	double time_start,time_end;

	//将采样信号存到Sequence[k]里面 
	for(k=0;k<N;k++)
	{
		t=k/8000.0;
		f=0.6*(sin(2*PI*500*t)+sin(2*PI*50*t));
		Squence[k]=f;
	}

	 time_start=clock();   //开始计时

	//进行DFT计算，并输出结果
	for(k=0;k<N;k++)
	{
		//对sum赋初值
		sum.real=0;               
	    sum.imag=0;

		//计算第k个点处的幅值
		for(n=0;n<N;n++)
		{
			one[n].real = cos(2*PI/N*k*n)*Squence[n];  //实部
			one[n].imag = -sin(2*PI/N*k*n)*Squence[n]; //虚部

			sum.real += one[n].real;                   //实部求和
			sum.imag += one[n].imag;                   //虚部求和
		}
		Amplitude[k] = sqrt(sum.real*sum.real+sum.imag*sum.imag);  //计算幅值 sqrt（a^2+b^2）
 
		printf("%.6f %.16f\n",k/8000.0,Amplitude[k]); //输出计算结果
	}

	time_end=clock();    //计时结束

	printf("DFT得到信号幅度谱的时间为%f ms\n",time_end-time_start);  //输出所需时间

	return 0;
}

运算结果：
在这里插入图片描述
FFT算法时间计算

在前面的FFT计算中加入计时，只需要加入time.h头文件，然后修改main函数如下：

  int main()
  { 
	  int i;
	  double t;
	  double pr[8192],pi[8192],fr[8192],fi[8192];
      
	  double time_start,time_end;


	  //产生抽样序列
	  for (i=0; i<8192; i++)
      { 
		 t=i/8192.0;
		 pr[i]=0.6*(sin(2*PI*500*t)+sin(2*PI*50*t));
		 pi[i]=0.0;
	  }
      
	  time_start=clock();   //开始计时
	  
	  kfft(pr,pi,8192,13,fr,fi,0,1);  //调用fft函数
 
      for (i=0; i<8192; i++)
      {
		  t=i/8192.0;
		  printf("%.6f %.16f\n",t,pr[i]); //输出
	  }
	  
	  time_end=clock();    //计时结束

	  printf("FFT得到信号幅度谱的时间为%.2f ms\n",time_end-time_start);  //输出所需时间

      return 0;
  }

运算结果：
在这里插入图片描述

FFT与DFT计算时间分析

根据上面的运行结果，我们可以看到，FFT计算明显比DFT计算所需的时间短，这与理论分析是一致的。

*pan

关注

7
点赞
踩
39

收藏

觉得还不错? 一键收藏
0
评论
C语言使用FFT得到信号幅度谱

前面的一篇文章我们介绍了使用DFT得到信号的幅度谱的方法，现在我们来看一下FFT实现信号幅度谱。这里我们使用的FFT源程序是徐士良老师的C语言算法程序——快速傅里叶变换。另外，本文也会对DFT以及FFT实现信号幅度谱计算时，所需要的计算量以及花费的时间进行比较。1.FFT原理FFT（Fast Fourier Transformation）是离散傅氏变换（DFT）的快速算法，也就是快速傅里叶变换...
复制链接

扫一扫