|||
一、Python,使用numpy.fft.fft或scipy的fftpack.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语言不容易调试,一但编写成功,速度和可用性还是有优势的。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-11-25 21:49
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社