求助傅立叶变换的c语言源代码，谢谢了。
本帖最后由 wdhd 于 2016-3-14 14:32 编辑谢谢大家，请大家多帮忙!!
回复：（游云）求助傅立叶变换的c语言源代码，谢谢了。...
本帖最后由 wdhd 于 2016-3-14 14:33 编辑1, DIT:
==============
1.1 Matlab code:
===============
function fft_2dit(n)
t = ;
m= log2(n);
%%%%%%%% raw data %%%%%%%%%%%
x = sin(t)/1 + sin(2*t)/2 - sin(3*t)/3 - sin(4*t)/4 + sin(5*t)/5;
y = t*0;
subplot(3,1,1)
plot(t, x, '-');
title('radix-2 DIT FFT');
xlabel(' t');
ylabel(' x');
grid;
%%%%%%%%%%% FFT %%%%%%%%%%%%%
%inverse bits
j = 1;
n1 = n - 1;
for i = 1:n1
if i < j
xt = x(j);
x(j) = x(i);
x(i) = xt;
xt = y(j);
y(j) = y(i);
y(i) = xt;
end
k = n / 2;
while k<j
j = j - k;
k = k / 2;
end
j = j + k;
end
%calculate FFT
for k = 1:m
n2 = 2.^(k-1);
n1 = 2.^k;
ee = -2*pi/n1;
for j = 1 : n2
a = (j-1)*ee;
c = cos(a);
s = sin(a);
for p = j:n1:n
q = p + n2;
xtt = x(p);
ytt = y(p);
xt = x(q)*c - y(q)*s;
yt = x(q)*s + y(q)*c;
x(p) = xtt + xt;
y(p) = ytt + yt;
x(q) = xtt - xt;
y(q) = ytt - yt;
end
end
end
subplot(3,1,2)
plot(t, x,'-r');
xlabel(' t');
ylabel('x');
grid;
%%%%%%%%%% IFFT %%%%%%%%%%
%inverse bits
j = 1;
n1 = n - 1;
for i = 1:n1
if i < j
xt = x(j);
x(j) = x(i);
x(i) = xt;
xt = y(j);
y(j) = y(i);
y(i) = xt;
end
k = n / 2;
while k<j
j = j - k;
k = k / 2;
end
j = j + k;
end
%calculate IFFT
n2 = 2;
for k = 1:m
n2 = 2.^(k-1);
n1 = 2.^k;
ee = 2*pi/n1;
for j = 1 : n2
a = (j-1)*ee;
c = cos(a);
s = sin(a); for p = j:n1:n
q = p + n2;
xtt = x(p);
ytt = y(p);
xt = x(q)*c - y(q)*s;
yt = x(q)*s + y(q)*c;
x(p) = xtt + xt;
y(p) = ytt + yt;
x(q) = xtt - xt;
y(q) = ytt - yt;
end
end
end
subplot(3,1,3)
plot(t, x/n, '-');
xlabel('t');
ylabel('x');
grid;
======================
1.2 C code:
======================
#define PI 3.1415926535
/**************************************************************
*
*
*
*
*
***************************************************************/
void fft_2dit(int n, double* x, double* y)
{
double a, e, xt, yt, c, s, xtt,ytt;
int n1, n2, i, j, k, m, q;
/*
* Calculate m such that n=2^m
*
* NOTE: If frexp() == 1, then frexp does not conform
* to the ANSI C spec of [0.5, 1)
*/
m = (int)(log10(n) / log10(2));
if( pow(2, m) != n)
{
printf("m must be 2's pow!\n");
}
/* --------------MAIN FFT LOOPS----------------------------- */
/* Parameter adjustments */
--y;
--x;
/* ------------DIGIT REVERSE COUNTER----------------- */
j = 1;
n1 = n - 1;
for(i= 1; i<= n1; i++)
{
if(i < j)
{
xt = x;
x = x;
x = xt;
xt = y;
y = y;
y = xt;
}
k = n / 2;
while(k < j)
{
j = j - k;
k = k / 2;
}
j = j + k;
}
//calculate FFT
for(k = 1; k <= m; k++)
{
n2 = 1<<(k-1);
n1 = 1<<k;
e = -2*PI/n1;
for(j = 1; j <= n2; j++)
{
a = (j-1)*e;
c = cos(a);
s = sin(a);
for(i = j; i <= n; i += n1)
{
q = i + n2;
xtt = x;
ytt = y;
xt = x*c - y*s;
yt = x*s + y*c;
x = xtt + xt;
y = ytt + yt;
x = xtt - xt;
y = ytt - yt;
}
}
}
}
/**************************************************************
*
*
*
*
*
***************************************************************/
void ifft_2dit(int n, double* x, double* y)
{
double a, e, xt, yt, c, s, xtt,ytt;
int n1, n2, i, j, k, m, q;
/*
* Calculate m such that n=2^m
*
* NOTE: If frexp() == 1, then frexp does not conform
* to the ANSI C spec of [0.5, 1)
*/
m = (int)(log10(n) / log10(2));
if( pow(2, m) != n)
{
printf("m must be 2's pow!\n");
}
/* --------------MAIN FFT LOOPS----------------------------- */
/* Parameter adjustments */
--y;
--x;
/* ------------DIGIT REVERSE COUNTER----------------- */
j = 1;
n1 = n - 1;
for(i= 1; i<= n1; i++)
{
if(i < j)
{
xt = x;
x = x;
x = xt;
xt = y;
y = y;
y = xt;
}
k = n / 2;
while(k < j)
{
j = j - k;
k = k / 2;
}
j = j + k;
}
//calculate IFFT
for(k = 1; k <= m; k++)
{
n2 = 1<<(k-1);
n1 = 1<<k;
e = 2*PI/n1;
for(j = 1; j <= n2; j++)
{
a = (j-1)*e;
c = cos(a);
s = sin(a);
for(i = j; i <= n; i += n1)
{
q = i + n2;
xtt = x;
ytt = y;
xt = x*c - y*s;
yt = x*s + y*c;
x = xtt + xt;
y = ytt + yt;
x = xtt - xt;
y = ytt - yt;
}
}
}
//each IFFT-ed data must be divided by 8
for(i = 1; i <= n; i++)
{
x = x/n;
y = y/n;
}
}
===================
2, DIF
===================
2.1 Matlab code:
===================
function fft_2dif(n)
m = log2(n);
t = ;
%%%%%%%%%% raw data %%%%%%%%
x = sin(t)/1 + sin(2*t)/2 - sin(3*t)/3 - sin(4*t)/4 + sin(5*t)/5;
y = t*0;
subplot(3,1,1)
plot(t, x, '-');
title('radix-2 DIF FFT');
xlabel(' t');
ylabel(' x');
grid;
%%%%%%%%%% FFT %%%%%%%%%%%%
n2 = n;
for k = 1:m
n1 = n2;
n2 = n2 / 2;
ee = -2*pi/n1;
a = 0.0;
for j = 1:n2
c = cos(a);
s = sin(a);
a = j*ee;
for p = j:n1:n
q = p + n2;
xt = x(p) - x(q);
yt = y(p) - y(q);
x(p) = x(p) + x(q);
y(p) = y(p) + y(q);
x(q) = c*xt - s*yt;
y(q) = c*yt + s*xt;
end
end
end
j = 1;
n1 = n - 1;
for i = 1:n1
if i < j
xt = x(j);
x(j) = x(i);
x(i) = xt;
xt = y(j);
y(j) = y(i);
y(i) = xt;
end
k = n / 2;
while k<j
j = j - k;
k = k / 2;
end
j = j + k;
end
subplot(3,1,2)
plot(t, x,'-r');
xlabel('t');
ylabel('x');
grid;
%%%%%%%%% IFFT %%%%%%%%%%%%%
n2 = n;
for k = 1:m
n1 = n2;
n2 = n2 / 2;
ee = 2*pi/n1;
a = 0.0;
for j = 1:n2
c = cos(a);
s = sin(a);
a = j*ee;
for p = j:n1:n
q = p + n2;
xt = x(p) - x(q);
yt = y(p) - y(q);
x(p) = x(p) + x(q);
y(p) = y(p) + y(q);
x(q) = c*xt - s*yt;
y(q) = c*yt + s*xt;
end
end
end
j = 1;
n1 = n - 1;
for i = 1:n1
if i < j
xt = x(j);
x(j) = x(i);
x(i) = xt;
xt = y(j);
y(j) = y(i);
y(i) = xt;
end
k = n / 2;
while k<j
j = j - k;
k = k / 2;
end
j = j + k;
end
subplot(3,1,3)
plot(t, x/n, '-');
xlabel(' t');
ylabel(' x');
grid;
======================
2.2 C code:
======================
/*****************************************************************************
*
* fft: FFT implementation with Cooley-Tukey radix-2 DIF
*
* Arguments: n -- data number to be FFT-ed
* x -- pointer to real section of original data
* y -- pointer to imaginary section of original data
*
* Notes: (original comments for the implementation)
* fft_rt.c - Cooley-Tukey radix-2 DIF FFT
* Complex input data in arrays x and y
* C.S. Burrus, Rice University, Sept. 1983
*
* Copyright 1995-2002 The MathWorks, Inc.
* $Revision: 1.15 $ $Date: 2002/04/12 22:18:20 $
*****************************************************************************
*/
void fft_2dif(int n, double *x, double *y)
{
double a, e, xt, yt, c, s;
int n1, n2, i, j, k, m, q;
/*
* Calculate m such that n=2^m
*
* NOTE: If frexp() == 1, then frexp does not conform
* to the ANSI C spec of [0.5, 1)
*/
m = (int)(log10(n) / log10(2));
if( pow(2, m) != n)
{
printf("m must be 2's pow!\n");
}
/* --------------MAIN FFT LOOPS----------------------------- */
/* Parameter adjustments */
--y;
--x;
/* Function Body */
n2 = n;
for (k = 1; k <= m; ++k)
{
n1 = n2;
n2 /= 2;
e = (double)-6.283185307179586476925286766559005768394 / n1;
a = 0.0;
for (j = 1; j <= n2; ++j)
{
c = cos(a);
s = sin(a);
a = j * e;
for (i = j;
i <= n;
i += n1
)
{
q = i + n2;
xt = x - x;
x += x;
yt = y - y;
y += y;
x = c * xt - s * yt;
y = c * yt + s * xt;
}
}
}
/* ------------DIGIT REVERSE COUNTER----------------- */
j = 1;
n1 = n - 1;
for (i = 1; i <= n1; ++i)
{
if (i < j)
{
xt = x;
x = x;
x = xt;
xt = y;
y = y;
y = xt;
}
k = n / 2;
while (k<j)
{
j -= k;
k /= 2;
}
j += k;
}
}
/*****************************************************************************
*
* ifft: IFFT implementation with Cooley-Tukey radix-2 DIF
*
* Arguments: n -- data number to be IFFT-ed
* x -- poiter to the real section of FFT-ed data
* y -- poiter to the imaginary section of FFT-ed data
*
* Notes: the only two difference between fft() and ifft() are:
* (1) e = (double) 6.283185307179586476925286766559005768394 / n1, in
fft()
* changed to
* e = (double)-6.283185307179586476925286766559005768394 / n1, in
iff();
* (2)each IFFT-ed data must be divied by n;
*****************************************************************************
*/
void ifft_2dif(int n, double *x, double *y)
{
double a, e, xt, yt, c, s;
int n1, n2, i, j, k, m, q;
/*
* Calculate m such that n=2^m
*
* NOTE: If frexp() == 1, then frexp does not conform
* to the ANSI C spec of [0.5, 1)
*/
m = (int)(log10(n) / log10(2));
if( pow(2, m) != n)
{
printf("m must be 2's pow!\n");
}
/* --------------MAIN FFT LOOPS----------------------------- */
/* Parameter adjustments */
--y;
--x;
/* Function Body */
n2 = n;
for (k = 1; k <= m; ++k) //step loops
{
n1 = n2;
n2 /= 2;
e = (double)6.283185307179586476925286766559005768394 / n1;
a = 0.0;
for (j = 1; j <= n2; ++j) //butterfly loops within each gr
oup
{
c = cos(a);
s = sin(a);
a = j * e; //theta for calculating scale fa
ctor
for (i = j; //group loops
i <= n;
i += n1
) //top input of a butterfly
{
q = i + n2; //bottom input of a butterfly
xt = x - x;
x += x;
yt = y - y;
y += y;
x = c * xt - s * yt;
y = c * yt + s * xt;
}
}
}
/* ------------DIGIT REVERSE COUNTER----------------- */
j = 1;
n1 = n - 1;
for (i = 1; i <= n1; ++i)
{
if (i < j)
{
xt = x;
x = x;
x = xt;
xt = y;
y = y;
y = xt;
}
k = n / 2;
while (k<j)
{
j -= k;
k /= 2;
}
j += k;
}
//each IFFT-ed data must be divided by n
for(i = 1; i <= n; i++)
{
x = x/n;
y = y/n;
}
}
这是matlab的，你用c实现一下就可以了
回复：（游云）求助傅立叶变换的c语言源代码，谢谢了。...
本帖最后由 wdhd 于 2016-3-14 14:34 编辑128点DIT FFT函数：
/* 采样来的数据放在dataR[ ]数组中，运算前dataI[ ]数组初始化为0 */
void FFT(float dataR[],float dataI[])
{int x0,x1,x2,x3,x4,x5,x6;
int L,j,k,b,p;
float TR,TI,temp;
/********** following code invert sequence ************/
for(i=0;i<128;i++)
{ x0=x1=x2=x3=x4=x5=x6=0;
x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01; x6=(i/64)&0x01;
xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;
dataI=dataR;
}
for(i=0;i<128;i++)
{ dataR=dataI; dataI=0; }
/************** following code FFT *******************/
for(L=1;L<=7;L++) { /* for(1) */
b=1; i=L-1;
while(i>0)
{b=b*2; i--；} /* b= 2^(L-1) */
for(j=0;j<=b-1;j++) /* for (2) */
{ p=1; i=7-L;
while(i>0) /* p=pow(2,7-L)*j; */
{p=p*2; i--；}
p=p*j;
for(k=j;k<128;k=k+2*b) /* for (3) */
{ TR=dataR; TI=dataI; temp=dataR;
dataR=dataR+dataR*cos_tab+dataI*sin_tab;
dataI=dataI-dataR*sin_tab+dataI*cos_tab;
dataR=TR-dataR*cos_tab-dataI*sin_tab;
dataI=TI+temp*sin_tab-dataI*cos_tab;
} /* END for (3) */
} /* END for (2) */
} /* END for (1) */
for(i=0;i<32;i++){ /* 只需要32次以下的谐波进行分析 */
w=sqrt(dataR*dataR+dataI*dataI);
w=w/64；}
w=w/2;
} /* END FFT */
split-radix DIF fft algorithm
本帖最后由 wdhd 于 2016-3-14 14:34 编辑/*----------------------------------------------------------------------
Routine msplfft:to perform the split-radix DIF fft algorithm.
input parameters:
xr,xi : float array. real and image of input signal.
n : the dimension of xr.
isign:if isign=-1 For Forward Transform
if isign=+1 For Inverse Transform.
output parameters:
xr,xi : float array. DFT result.
Notes:
n must be power of 2.
----------------------------------------------------------------------*/
void CAnalyze::MsplFft(float *xr,float *xi,int n,int isign)
{
_complex xt;
float es,e,a,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3;
int m,n2,k,n4,j,is,id,i1,i2,i3,i0,n1,i,nn;
for(m=1;m<=16;m++)
{
nn=(int)pow(2,m);
if(n==nn)break;
}
if(m>16)
{
AfxMessageBox(" N is not a power of 2 ! \n");
return;
}
n2=n*2;
es=(float)(-isign*atan(1.0f)*8.0f);
for(k=1;k<m;k++)
{
n2=n2/2;
n4=n2/4;
e=es/n2;
a=0.0;
for(j=0;j<n4;j++)
{
a3=3*a;
cc1=(float)cos(a);
ss1=(float)sin(a);
cc3=(float)cos(a3);
ss3=(float)sin(a3);
a=(j+1)*e;
is=j;
id=2*n2;
do
{
for(i0=is;i0<n;i0+=id)
{
i1=i0+n4;
i2=i1+n4;
i3=i2+n4;
r1=xr-xr;
s1=xi-xi;
r2=xr-xr;
s2=xi-xi;
xr+=xr;xi+=xi;
xr+=xr;xi+=xi;
if(isign!=1)
{
s3=r1-s2;
r1+=s2;
s2=r2-s1;
r2+=s1;
}
else
{
s3=r1+s2;
r1=r1-s2;
s2=-r2-s1;
r2=-r2+s1;
}
xr=r1*cc1-s2*ss1;
xi=-s2*cc1-r1*ss1;
xr=s3*cc3+r2*ss3;
xi=r2*cc3-s3*ss3;
}
is=2*id-n2+j;
id=4*id;
}while(is<n-1);
}
}
/* ------------ special last stage -------------------------*/
is=0;
id=4;
do
{
for(i0=is;i0<n;i0+=id)
{
i1=i0+1;
xt.x=xr;
xt.y=xi;
xr=(float)(xt.x+xr);
xi=(float)(xt.y+xi);
xr=(float)(xt.x-xr);
xi=(float)(xt.y-xi);
}
is=2*id-2;
id=4*id;
}while(is<n-1);
j=1;
n1=n-1;
for(i=1;i<=n1;i++)
{
if(i<j)
{
xt.x=xr;
xt.y=xi;
xr=xr;
xi=xi;
xr=(float)xt.x;
xi=(float)xt.y;
}
k=n/2;
do
{
if(k>=j)
break;
j-=k;
k/=2;
}while(1);
j+=k;
}
if(isign==-1)
return;
for(i=0;i<n;i++)
{
xr/=(float)n;
xi/=(float)(n);
}
}
页:
[1]