冪冪源碼
❶ 我用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
