cambaluc的个人博客分享 http://blog.sciencenet.cn/u/cambaluc

博文

用Python、Matlab、C实现傅立叶变换FFT()

已有 4553 次阅读 2020-2-2 11:56 |个人分类:算法|系统分类:科研笔记| FFT, Python

一、Python,使用numpy.fft.fftscipyfftpack.fft:

import numpy as np

t=np.array([0.300, 0.550, 1.119, 0.538, 0.857, 1.069, 0.369, 0.206, 1.000,

            0.550, 1.119, 1.538, 1.557, 1.069, 1.369, 0.606])

ff=np.fft.fft(t)

print(ff)


from scipy import fftpack

f=fftpack.fft(t)

for i in f:

print("{0.real:.4f}{0.imag:+.4f}j".format(i)) 

输出结果:

[ 1.38160000e+01+0.j         -6.02483817e-03+2.48405969j

 -2.74175981e+00-0.65980613j -3.30153876e-01-0.00602484j

 -2.62000000e-01-0.35j       -2.48405969e+00-0.0202384j

  5.13759810e-01+0.34019387j  2.02384005e-02-0.33015388j

  1.56400000e+00+0.j          2.02384005e-02+0.33015388j

  5.13759810e-01-0.34019387j -2.48405969e+00+0.0202384j

 -2.62000000e-01+0.35j       -3.30153876e-01+0.00602484j

 -2.74175981e+00+0.65980613j -6.02483817e-03-2.48405969j]

13.8160+0.0000j

-0.0060+2.4841j

-2.7418-0.6598j

-0.3302-0.0060j

-0.2620-0.3500j

-2.4841-0.0202j

0.5138+0.3402j

0.0202-0.3302j

1.5640+0.0000j

0.0202+0.3302j

0.5138-0.3402j

-2.4841+0.0202j

-0.2620+0.3500j

-0.3302+0.0060j

-2.7418+0.6598j

-0.0060-2.4841j

 

二、用Matlab

A=[0.300; 0.550; 1.119; 0.538; 0.857; 1.069; 0.369; 0.206; 1.000;0.550; 1.119; 1.538; 1.557; 1.069; 1.369; 0.606]

B=fft(A)

输出:

B =

  13.8160         

  -0.0060 + 2.4841i

  -2.7418 - 0.6598i

  -0.3302 - 0.0060i

  -0.2620 - 0.3500i

  -2.4841 - 0.0202i

   0.5138 + 0.3402i

   0.0202 - 0.3302i

   1.5640         

   0.0202 + 0.3302i

   0.5138 - 0.3402i

  -2.4841 + 0.0202i

  -0.2620 + 0.3500i

  -0.3302 + 0.0060i

  -2.7418 + 0.6598i

  -0.0060 - 2.4841i

 

三、用C语言

// fft20200202.cpp : kangjian, console application.

#include "stdafx.h"

#include <math.h>

#include <malloc.h>

#define PI 3.1415926535

typedef struct{   double R;   double I;

}COMPLEX;

COMPLEX add(COMPLEX a,COMPLEX b)

{    COMPLEX c;

     c.R=a.R+b.R;     c.I=a.I+b.I;       return c;

}

COMPLEX sub(COMPLEX a,COMPLEX b)

{    COMPLEX c;

     c.R=a.R-b.R;     c.I=a.I-b.I;         return c;

}

COMPLEX mul(COMPLEX a,COMPLEX b)

{    COMPLEX c;

     c.R=a.R*b.R-a.I*b.I;     c.I=a.R*b.I+b.R*a.I;

          return c;

}

int nx(int k,int p)

{   int q,t=0;

         for(q=0;q<p;q++)  

                            if((1<<q)&k)   

                                     t+=1<<(p-q-1);

         return t;

}

void kj_FFT(COMPLEX *T,COMPLEX *F,int p)

{  int N,q,k,r,j,i,t;

  COMPLEX *A1,*A2, *W;

  double a;

    N=1<<p;

    A1=(COMPLEX *)malloc(sizeof(COMPLEX)*N);

    A2=(COMPLEX *)malloc(sizeof(COMPLEX)*N);

    W =(COMPLEX *)malloc(sizeof(COMPLEX)*N/2);

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

      A1[i]=T[i];

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

                   {  a=-i*2*PI/N; 

                      W[i].R=cos(a);  

                      W[i].I=sin(a);

                   }

   for(q=0;q<p;q++)

            {for(k=0;k<1<<q;k++)

                      { r=1<<(p-q-1);

                              for(j=0;j<r;j++)

                                     { t=k*( 1<<(p-q)  );

                                       A2[j+t]  =add(A1[j+t], mul(A1[j+t+r],W[ nx((j+t)>>(p-q-1),p)]) );//

                                       A2[j+t+r]=sub(A1[j+t], mul(A1[j+t+r],W[ nx((j+t)>>(p-q-1),p)]) );

                                      }

                      }

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

                            A1[i]=A2[i];

            }

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

   { 

            F[i].R=A2[nx(i,p)].R/N;  

       F[i].I=A2[nx(i,p)].I/N;

   }

}

int testfft()

{           int n=16,i,p=4;

        double y[16]={0.300, 0.550, 1.119, 0.538, 0.857, 1.069, 0.369, 0.206, 1.000,

            0.550, 1.119, 1.538, 1.557, 1.069, 1.369, 0.606};

            COMPLEX  *pT,*pF;

       pT=(COMPLEX *)malloc(sizeof(COMPLEX)*n);//no free

       pF=(COMPLEX *)malloc(sizeof(COMPLEX)*n);

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

            {   pT[i].R=y[i];

                      pT[i].I=0.;

            }

             kj_FFT(pT,pF,p);

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

                            {   pF[i].R=pF[i].R*n; //  python,no /N,

                                     pF[i].I=pF[i].I*n;

                printf("%.4lf,%+.4lfj\n" ,pF[i].R, pF[i].I );

                            }

       return 0;

}

int main(int argc, char* argv[])

{           testfft();

       return 0;

}

 

输出结果:

13.8160,+0.0000j

-0.0060,+2.4841j

-2.7418,-0.6598j

-0.3302,-0.0060j

-0.2620,-0.3500j

-2.4841,-0.0202j

0.5138,+0.3402j

0.0202,-0.3302j

1.5640,+0.0000j

0.0202,+0.3302j

0.5138,-0.3402j

-2.4841,+0.0202j

-0.2620,+0.3500j

-0.3302,+0.0060j

-2.7418,+0.6598j

-0.0060,-2.4841j

 

C语言不容易调试,一但编写成功,速度和可用性还是有优势的。




https://blog.sciencenet.cn/blog-797552-1216569.html

上一篇:让EXCEL帮你念成绩单,录入校对
下一篇:准备上网课时用到的小程序段
收藏 IP: 114.242.250.*| 热度|

1 郑泽

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-11-25 21:49

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部