有感FFT.------频谱的分析

[lanren_me]

问FFT很多情况都是知其然而不只其然。
看了点资料，不敢独享，特此奉献。

void __fastcall four1(double data[65],
int nn, //复数个数;
int isign)
/*
零频率分量f0的实部和虚部在data[1],data[2]中；最小的非零正频率的实部和虚部在data[3],data[4];
数值上最小的非零频率的实部和虚部在data[2*nn-1]和data[2*nn]
数量上增长的正频率存储在实-虚对data[5]-data[6],直到data[nn-1];data[nn];
............负.................data[2*nn-3],data[2*nn-2];直到data[nn+3];data[nn+4];
最后，data[nn-1]和data[nn+2]存储一个混叠点的实部和虚部;其包含最大的正频率和最小的频率


-----

输入的数组包含N(2的幂次)个事件上的取样，形成结长为2N的实形数组,其实部和虚部交其相隔;

输出数组是包含N个频率值上的复数FFT频谱,实部和虚部也相隔交替;
数组从零频率开始,逐渐排列最大正频率,然后紧接负频率，从次最小的负频率一直升到到恰好位于零以下的负频率;

-----
例如:(1,2)--对应---f=0;
(3,4).........f=1/(N*取样间隔);
(n-1,n).......f=(n/2-1)/(N*取样间隔);
(n+1,n+2).....f=(+/-)1/(2*取样间隔);
(n+3,n+4).....f=-(n/2-1)/(N*取样间隔);
(2n-1,2n).....f=-1/(N*取样间隔);

-----
当输入的数组是下标0开始 ,即范围为0..2*n-1;
调用时必须将data的指针-1; 例如:four1(data-1,1024,1) ,
这时实部返回在data[0],虚部在data[1];以次类推;

----正常情况
isign=1; 则置后的data[1]..data[2*nn] 是输入的离散fft;
isign=-1; ........逆变换..
data是2*nn的实数数组;
nn必须为2的幂;




*/

{
int n,j,i,m,mmax,istep;
double tempr,tempi,theta,wpr,wpi,wr,wi,wtemp;
n = 2 * nn;
j = 1;
for (i = 1; i<=n ; i=i+2) //整理成位序颠倒的次序；
{
if( j > i)
{
tempr = data[j];
tempi = data[j + 1];
data[j] = data[i];
data[j + 1] = data[i + 1];
data[i] = tempr;
data[i + 1] = tempi;
}
m = n / 2;
while (m >= 2 && j > m)
{
j = j - m;
m = m / 2;
}
j = j + m;
}
mmax = 2;
while( n > mmax ) //log2N次的外循环
{
istep = 2 * mmax;
theta = 6.28318530717959 / (isign * mmax); //关于三角递归的初始化;
wpr = -2.0 * sin(0.5 * theta)*sin(0.5 * theta);
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
for( m = 1; m<=mmax; m=m+2)
{
for (i = m; i<=n; i=i+istep)
{
j = i + mmax; //这是丹尼尔深--兰左司公式;
tempr=double(wr)*data[j]-double(wi)*data[j+1];
tempi=double(wr)*data[j+1]+double(wi)*data[j];
data[j] = data[i] - tempr;
data[j + 1] = data[i + 1] - tempi;
data[i] = data[i] + tempr;
data[i + 1] = data[i + 1] + tempi;
}
wtemp = wr; // 三角递归
wr = wr * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
mmax = istep;
}
}

[lanren_me]

1111