幂幂源码
❶ 我用vb6.0编了个计算器 乘方问题不会编了请高手给我源码
a^b 就是a的b次方
log要用换弟公式
❷ 急！那位大侠有用C语言实现的互相关算法的源代码，谢谢啦
#include <math.h>
#define M_PI 3.14159265358979323846
#define FALSE 0
#define TRUE 1
#define BIG 1e10
#define SMALL 1e-10
typedef struct {
float r, i;
} complex;
/* FAST CORRELATION OF X(0:L) AND Y(0:L). FINDS RXY(0) THRU RXY(NMAX). */
/* L=LAST INDEX IN BOTH X AND Y. MUST BE (POWER OF 2)+1 AND AT LEAST 5. */
/* ITYPE=TYPE OF CORRELATION=0 IF X AND Y ARE THE SAME VECTOR (AUTO- */
/* CORRELATION), OR NOT 0 IF X AND Y ARE DIFFERENT VECTORS. */
/* NMAX=MAXIMUM LAG OF INTEREST IN THE CORRELATION FUNCTION. */
/* FFT LENGTH ,N, USED INTERNALLY, IS L-1. */
/* LET K=INDEX OF FIRST NONZERO SAMPLE IN Y(0)---Y(N-1). THEN X(0) */
/* 到 X(N-1) MUST INCLUDE PADDING OF AT LEAST NMAX-K ZEROS. */
/* CORRELATION FUNCTION, RXY, REPLACES X(0) THRU X(NMAX). */
/* Y(0) THRU Y(L) IS REPLACED BY ITS FFT, COMPUTED USING SPFFTR. */
/* IERROR=0 NO ERROR DETECTED */
/* 1 L-1 NOT A POWER OF 2 */
/* 2 NMAX OUT OF RANGE */
/* 3 INADEQUATE ZERO */
void spcorr(float *x, float *y, long *l, long *type, long *nmax, long *error)
/*
x:序列X；
y：序列Y；
l:序列X与序列Y的长度，不小5，且要为2的幂次方；
type:相关的类型，0：表示X与Y序列相同，其它值：X与Y序列不相同
nmax:相关的最大时延；
error:运行出错提示；0：无错；1：数据长度不是2的幂次方；2：时延超界；3：无足够零填充出错
*/
{
long j, k, m, n;//n:FFT长度；k:序列Y中的首个非零样本的位置序号；在序列Y中必须最少包含有(nmax-k)零填充。
complex cx;
float test;
n = *l - 1;
if (*nmax < 0 || *nmax >= n)
{
*error = 2;
return;
}
test = (float) n;
test /= 2.0;
while ((test - 2.0) > 0.0)
{
test /= 2.0;
}
if ((test - 2.0) == 0)
{
for (k = 0 ; k < n && y[k] == 0.0 ; ++k) ;
for (j = n - 1 ; j >= 0 && x[j] == 0.0 ; --j) ;
if ((n - 1 - j) < (*nmax - k))
{
*error = 3;
return;
}
spfftr(x, &n);//对X序列FFT变换
if (*type != 0)
{
spfftr(y, &n);//如果X、Y是相同序列，则对Y序列也进行FFT
}
for (m = 0 ; m <= (n / 2) ; ++m)
{
cx.r = x[m * 2] * y[m * 2] - -x[(m * 2) + 1] * y[(m * 2) + 1];
cx.i = x[m * 2] * y[(m * 2) + 1] + -x[(m * 2) + 1] * y[m * 2];
x[m * 2] = cx.r / n;
x[(m * 2) + 1] = cx.i / n;
}
spiftr(x, &n);
*error = 0;
}
else if ((test - 2.0) < 0.0)
{
*error = 1;
}
return;
} /* spcorr */
/* SPFFTR 11/12/85 */
/* FFT ROUTINE FOR REAL TIME SERIES (X) WITH N=2**K SAMPLES. */
/* COMPUTATION IS IN PLACE, OUTPUT REPLACES INPUT. */
/* INPUT: REAL VECTOR X(0:N+1) WITH REAL DATA SEQUENCE IN FIRST N */
/* ELEMENTS; ANYTHING IN LAST 2. NOTE: X MAY BE DECLARED */
/* REAL IN MAIN PROGRAM PROVIDED THIS ROUTINE IS COMPILED */
/* SEPARATELY ... COMPLEX OUTPUT REPLACES REAL INPUT HERE. */
/* OUTPUT: COMPLEX VECTOR XX(O:N/2), SUCH THAT X(0)=REAL(XX(0)),X(1)= */
/* IMAG(XX(0)), X(2)=REAL(XX(1)), ..., X(N+1)=IMAG(XX(N/2). */
/* IMPORTANT: N MUST BE AT LEAST 4 AND MUST BE A POWER OF 2. */
//FFT计算函数
void spfftr(complex *x, long *n)
{
/* Builtin functions */
void r_cnjg();
/* Local variables */
void spfftc();
long m, tmp_int;
complex u, tmp, tmp_complex;
float tpn, tmp_float;
tpn = (float) (2.0 * M_PI / (double) *n);
tmp_int = *n / 2;
spfftc(x, &tmp_int, &neg_i1);
x[*n / 2].r = x[0].r;
x[*n / 2].i = x[0].i;
for (m = 0 ; m <= (*n / 4) ; ++m)
{
u.r = (float) sin((double) m * tpn);
u.i = (float) cos((double) m * tpn);
r_cnjg(&tmp_complex, &x[*n / 2 - m]);
tmp.r = (((1.0 + u.r) * x[m].r - u.i * x[m].i)
+ (1.0 - u.r) * tmp_complex.r - -u.i * tmp_complex.i) / 2.0;
tmp.i = (((1.0 + u.r) * x[m].i + u.i * x[m].r)
+ (1.0 - u.r) * tmp_complex.i + -u.i * tmp_complex.r) / 2.0;
tmp_float = ((1.0 - u.r) * x[m].r - -u.i * x[m].i
+ (1.0 + u.r) * tmp_complex.r - u.i * tmp_complex.i) / 2.0;
x[m].i = ((1.0 - u.r) * x[m].i + -u.i * x[m].r
+ (1.0 + u.r) * tmp_complex.i + u.i * tmp_complex.r) / 2.0;
x[m].r = tmp_float;
r_cnjg(&x[*n / 2 - m], &tmp);
}
return;
} /* spfftr */
/* SPIFTR 02/20/87 */
/* INVERSE FFT OF THE COMPLEX SPECTRUM OF A REAL TIME SERIES. */
/* X AND N ARE THE SAME AS IN SPFFTR. IMPORTANT: N MUST BE A POWER */
/* OF 2 AND X MUST BE DIMENSIONED X(0:N+1) (REAL ARRAY, NOT COMPLEX). */
/* THIS ROUTINE TRANSFORMS THE OUTPUT OF SPFFTR BACK INTO THE INPUT, */
/* SCALED BY N. COMPUTATION IS IN PLACE, AS IN SPFFTR. */
//逆FFT变换函数
void spiftr(complex *x, long *n)
{
long m, tmp_int;
complex u, tmp_complex, tmp;
float tpn, tmp_float;
tpn = (float) (2.0 * M_PI / (double) *n);
for (m = 0 ; m <= (*n / 4) ; ++m)
{
u.r = (float) sin((double) m * tpn);
u.i = (float) -cos((double) m * tpn);
r_cnjg(&tmp_complex, &x[*n / 2 - m]);
tmp.r = ((1.0 + u.r) * x[m].r - u.i * x[m].i)
+ ((1.0 - u.r) * tmp_complex.r - -u.i * tmp_complex.i);
tmp.i = ((1.0 + u.r) * x[m].i + u.i * x[m].r)
+ ((1.0 - u.r) * tmp_complex.i + -u.i * tmp_complex.r);
r_cnjg(&tmp_complex, &x[*n / 2 - m]);
tmp_float = ((1.0 - u.r) * x[m].r - -u.i * x[m].i)
+ ((1.0 + u.r) * tmp_complex.r - u.i * tmp_complex.i);
x[m].i = ((1.0 - u.r) * x[m].i + -u.i * x[m].r)
+ ((1.0 + u.r) * tmp_complex.i + u.i * tmp_complex.r);
x[m].r = tmp_float;
r_cnjg(&x[*n / 2 - m], &tmp);
}
tmp_int = *n / 2;
spfftc(x, &tmp_int, &pos_i1);
return;
} /* spiftr *
void r_cnjg(complex *r, complex *z)
{
r->r = z->r;
r->i = -z->i;
}
❸ FFT的源码含义
在C环境下的源码
源码（1）：
注：亲测，这个版本无法运行，作者删改了重要内容 请参考源码（2） //快速傅立叶变换//入口参数：//l:l=0,傅立叶变换;l=1,逆傅立叶变换//il:il=0,不计算傅立叶变换或逆变换模和幅角；il=1,计算模和幅角//n:输入的点数，为偶数，一般为32，64，128，...,1024等//k:满足n=2^k(k>0),实质上k是n个采样数据可以分解为偶次幂和奇次幂的次数//pr[]:l=0时，存放N点采样数据的实部//l=1时,存放傅立叶变换的N个实部//pi[]:l=0时，存放N点采样数据的虚部//l=1时,存放傅立叶变换的N个虚部////出口参数：//fr[]:l=0,返回傅立叶变换的实部//l=1,返回逆傅立叶变换的实部//fi[]:l=0,返回傅立叶变换的虚部//l=1,返回逆傅立叶变换的虚部//pr[]:il=1,i=0时，返回傅立叶变换的模//il=1,i=1时，返回逆傅立叶变换的模//pi[]:il=1,i=0时，返回傅立叶变换的辐角//il=1,i=1时，返回逆傅立叶变换的辐角voidfft(doublepr[],doublepi[],intn,intk,doublefr[],doublefi[],intl,intil){intit,m,is,i,j,nv,l0;doublep,q,s,vr,vi,poddr,poddi;for(it=0;it<=n-1;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/n;pr[1]=cos(p);pi[1]=-sin(p);if(l)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=p-q;pi=s-p-q;}for(it=0;it<=n-2;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/=2;nv<<=1;for(it=0;it<=(m-1)*nv;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*=(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]+=poddr;fi[it+j]+=poddi;}}if(l)for(i=0;i<=n-1;fr/=n,fi[i++]/=n);if(il)for(i=0;i<=n-1;i++){pr=sqrt(fr*fr+fi*fi);if(fabs(fr)<0.000001*fabs(fi))pi=fi*fr>0?90.0-90.0;elsepi=atan(fi/fr)*360.0/6.283185306;}return;}源码（2）
ps:可以运行的 //.htocorrectlydefineM_PI#define_USE_MATH_DEFINES#include<math.h>#include<stdio.h>#include<stdlib.h>#definePIM_PI/*pitomachineprecision,definedinmath.h*/#defineTWOPI(2.0*PI)/*FFT/IFFTroutine.(seepages507-508ofNumericalRecipesinC)Inputs:data[]:arrayofcomplex*datapointsofsize2*NFFT+1.data[0]isunused,*then'thcomplexnumberx(n),for0<=n<=length(x)-1,isstoredas:data[2*n+1]=real(x(n))data[2*n+2]=imag(x(n))iflength(Nx)<NFFT,:FFTorderNFFT.ThisMUSTbeapowerof2and>=length(x).isign:ifsetto1,computestheforwardFFTifsetto-1,computesInverseFFT-/NFFT.Outputs:data[]:,overwritingtheinput.*/voidfour1(doubledata[],intnn,intisign){intn,mmax,m,j,istep,i;doublewtemp,wr,wpr,wpi,wi,theta;doubletempr,tempi;n=nn<<1;j=1;for(i=1;i<n;i+=2){if(j>i){tempr=data[j];data[j]=data[i];data[i]=tempr;tempr=data[j+1];data[j+1]=data[i+1];data[i+1]=tempr;}m=n>>1;while(m>=2&&j>m){j-=m;m>>=1;}j+=m;}mmax=2;while(n>mmax){istep=2*mmax;theta=TWOPI/(isign*mmax);wtemp=sin(0.5*theta);wpr=-2.0*wtemp*wtemp;wpi=sin(theta);wr=1.0;wi=0.0;for(m=1;m<mmax;m+=2){for(i=m;i<=n;i+=istep){j=i+mmax;tempr=wr*data[j]-wi*data[j+1];tempi=wr*data[j+1]+wi*data[j];data[j]=data[i]-tempr;data[j+1]=data[i+1]-tempi;data[i]+=tempr;data[i+1]+=tempi;}wr=(wtemp=wr)*wpr-wi*wpi+wr;wi=wi*wpr+wtemp*wpi+wi;}mmax=istep;}}在C++环境下的源码
boolFFT(complex<double>*TD,complex<double>*FD,intr){//一维快速Fourier变换。//complex<double>*TD——指向时域数组的指针;complex<double>*FD——指向频域数组的指针;r——2的幂数，即迭代次数LONGcount;//Fourier变换点数inti,j,k;//循环变量intbfsize,p;//中间变量doubleangle;//角度complex<double>*W,*X1,*X2,*X;count=1<<r;//计算Fourier变换点数为1左移r位W=newcomplex<double>[count/2];X1=newcomplex<double>[count];X2=newcomplex<double>[count];//分配运算所需存储器//计算加权系数（旋转因子w的i次幂表）for(i=0;i<count/2;i++){angle=-i*PI*2/count;W[i]=complex<double>(cos(angle),sin(angle));}//将时域点写入X1memcpy(X1,TD,sizeof(complex<double>)*count);//采用蝶形算法进行快速Fourier变换for(k=0;k<r;k++){for(j=0;j<1<<k;j++){bfsize=1<<(r-k);for(i=0;i<bfsize/2;i++){p=j*bfsize;X2[i+p]=X1[i+p]+X1[i+p+bfsize/2]*W[i*(1<<k)];X2[i+p+bfsize/2]=X1[i+p]-X1[i+p+bfsize/2]*W[i*(1<<k)];}}X=X1;X1=X2;X2=X;}//重新排序for(j=0;j<count;j++){p=0;for(i=0;i<r;i++){if(j&(1<<i)){p+=1<<(r-i-1);}}FD[j]=X1[p];}//释放内存deleteW;deleteX1;deleteX2;returntrue;}在Matlab环境下的源码
function X=myfft(x)
N=length(x);
if N==1
X=x;
return;
end
%myfft函数 用递归实现
t=log2(N);
t1=floor(t);
t2=ceil(t);
if t1~=t2; %若x的长度N不为2的整数次幂，则补0至最接近的2的整数次幂
x=[x zeros(1,2^t2-N)];
N=2^t2;
end
w0=exp(-1i*2*pi/N);
X=zeros(1,N);
n=1:N/2;
xo(n)=x(2*n-1);
xe(n)=x(2*n);
XO=myfft(xo); %递归调用
XE=myfft(xe);
for n=1:N/2
X(n)=XO(n)+XE(n)*(w0^(n-1));
X(n+N/2)=XO(n)-XE(n)*(w0^(n-1));
end
