FFT 要求輸入的 時序 點的個數 是 2 的整數次 方，例如 1024，2048 ...
int get_pow_2(int NP)
{
int j = 1, k,m = 1;
while ( m < NP ) { m = m * 2; j = j + 1; }
return m;
}
妳需要確定，是重新采樣，把原時序點的個數變 2 的整數次 方 或 後面添加，添加方法為 全0 ，或是假定 時序重復，即添加 原時序的點序直到 滿足 2 的整數次 方。
求時序平均值，把所有點去掉平均值（即平移中線）。
時序濾波。濾波窗考慮：矩形，cos形，或其他。
確定時序 dt, 變換後的df, 自由度或通頻帶
做 FFT
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void jfour1(float ya[], unsigned long nn, int isign)
{
unsigned long n,mmax,m,j,istep,i;
double wtemp,wr,wpr,wpi,wi,theta;
float tempr,tempi;
n=nn << 1;
j=1;
for (i=1;i<n;i+=2) {
if (j > i) {
SWAP(ya[j],ya[i]);
SWAP(ya[j+1],ya[i+1]);
}
m=n >> 1;
while (m >= 2 && j > m) {
j -= m;
m >>= 1;
};
j += m;
};
mmax=2;
while (n > mmax) {
istep=mmax << 1;
theta=isign*(6.28318530717959/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 * ya[j]- wi * ya[j+1];
tempi = wr * ya[j+1] + wi * ya[j];
ya[j] = ya[i] - tempr;
ya[j+1] = ya[i+1] - tempi;
ya[i] += tempr;
ya[i+1] += tempi;
};
wr = (wtemp=wr) * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
};
mmax=istep;
};
}
#undef SWAP
void jrealft(float ya[], unsigned long n, int isign)
{
void jfour1(float ya[], unsigned long nn, int isign);
unsigned long i,i1,i2,i3,i4,np3,n05;
float c1=0.5,c2,h1r,h1i,h2r,h2i;
double wr,wi,wpr,wpi,wtemp,theta;
n05 = n >> 1;
theta=3.141592653589793/(double) (n05);
if (isign == 1) {
c2 = -0.5;
jfour1(ya,n05,1);
} else {
c2=0.5;
theta = -theta;
};
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0+wpr;
wi=wpi;
np3=n+3;
for (i=2;i<=(n>>2);i++) {
i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
h1r = c1 * (ya[i1] + ya[i3]);
h1i = c1 * (ya[i2] - ya[i4]);
h2r = -c2* (ya[i2] + ya[i4]);
h2i = c2 * (ya[i1] - ya[i3]);
ya[i1] = h1r + wr * h2r - wi * h2i;
ya[i2] = h1i + wr * h2i + wi * h2r;
ya[i3] = h1r - wr * h2r + wi * h2i;
ya[i4] = -h1i + wr * h2i + wi * h2r;
wr= (wtemp=wr) * wpr - wi * wpi + wr;
wi=wi * wpr + wtemp * wpi + wi;
};
if (isign == 1) {
ya[1] = (h1r=ya[1]) + ya[2];
ya[2] = h1r-ya[2];
} else {
ya[1] = c1 * ((h1r=ya[1]) + ya[2]);
ya[2]=c1 * (h1r - ya[2]);
jfour1(ya,n05,-1);
}
}
頻譜濾波（去掉截止頻率外的頻譜泄漏）。