谁有FFT源代码????

[grgz]

因为时间紧迫,没有时间慢慢研究,谁有FFT源代码?先谢啦.

[grgz]

小波算法也可以

[grgz]

采集到一段波形,要取出14~30Hz的波形,请问如何做?

[会思考的草]

/****************************************************
File : fft.c

Description:
process the FFT and Inverse FFT.

Function:
global:
fft(): process the FFT and Inverse FFT.
****************************************************/
#include <math.h>
#include "typedef.h"
#include "const.h"

/****************************************************
Function: fft()

Description:process the FFT and Inverse FFT.

Inputs:
FLOAT * datam1: Data to be used in FFT or IFFT
int nn: the FFT SIZE
int isign: the sign index of FFT or IFFT
-1:----- FFT
1 :----- IFFT

Outputs:
FLOAT * datam1: Data after FFT or IFFT transform

Return Value:
none

Algorithm Description:
Subroutine FFT: Fast Fourier Transform
* Replaces data by its DFT, if isign is 1, or replaces data *
* by inverse DFT times nn if isign is -1. data is a complex *
* array of length nn, input as a real array of length 2*nn. *
* nn MUST be an integer power of two. This is not checked *
* The real part of the number should be in the zeroeth *
* of data , and the imaginary part should be in the next *
* element. Hence all the real parts should have even indeces *
* and the imaginary parts, odd indeces. *

* Data is passed in an array starting in position 0, but the *
* code is copied from Fortran so uses an internal pointer *
* which accesses position 0 as position 1, etc. *

Improtance:!!!
* This code uses e+jwt sign convention, so isign should be *
* reversed for e-jwt. *
****************************************************/

#define SWAP(a,b) tempr = (a);(a) = (b); (b) = tempr

void fft(FLOAT *datam1,int nn,int isign)

{
int n,mmax,m,j,istep,i;
double wtemp,wr,wpr,wpi,wi,theta;
FLOAT register tempr,tempi;
FLOAT *data;

/* Use pointer indexed from 1 instead of 0 */
data = &datam1[-1];

n = nn << 1;
j = 1;
for( i = 1; i < n; i+=2 ) {
if ( j > i) {
SWAP(data[j],data[i]);
SWAP(data[j+1],data[i+1]);
}
m = n >> 1;
while ( m >= 2 && j > m ) {
j -= m;
m >>= 1;
}
j += m;
}
mmax = 2;
while ( n > mmax) {
istep = 2 * mmax;
theta = 6.28318530717959/(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;
}
}

[grgz]

我总共采集了1024个点,放在OnTime[1024]里面,调用的时候是不是:
fft(OnTime,1024,1)
呢?但是运行到
SWAP(data[j],data[i]);
的时候有错,是数据越界.后来把nn改为256才可以.
还用,最后的结果存放在OnTime里面,里面的数据是表示什么的呢?

[grgz]

up

[CSwain]

帮你up, 我也在找这方面的资料。

http://expert.csdn.net/Expert/TopicView3.asp?id=2490057

[wqs6]

#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>
#include<stdio.h>

void four1(double data[65], int nn, int isign)
{
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 )
{
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;
}
}


void twofft(double data1[], double data2[], double fft1[],
double fft2[], int& n)
{
int j,n2,j2;
double c1r,c1i,c2r,c2i,conjr,conji,h1r,h1i,h2r,h2i;
c1r = 0.5;
c1i = 0.0;
c2r = 0.0;
c2i = -0.5;
for (j = 1; j<=n; j++)
{
fft1[2 * j - 1] = data1[j];
fft1[2 * j] = data2[j];
}
four1(fft1, n, 1);
fft2[1] = fft1[2];
fft2[2] = 0.0;
fft1[2] = 0.0;
n2 = 2 * (n + 2);
for (j = 2; j<=n / 2 + 1; j++)
{
j2 = 2 * j;
conjr = fft1[n2 - j2 - 1];
conji = -fft1[n2 - j2];
h1r=c1r*(fft1[j2-1]+conjr)-c1i*(fft1[j2]+conji);
h1i=c1i*(fft1[j2-1]+conjr)+c1r*(fft1[j2]+conji);
h2r=c2r*(fft1[j2-1]-conjr)-c2i*(fft1[j2]-conji);
h2i=c2i*(fft1[j2-1]-conjr)+c2r*(fft1[j2]-conji);
fft1[j2 - 1] = h1r;
fft1[j2] = h1i;
fft1[n2 - j2 - 1] = h1r;
fft1[n2 - j2] = -h1i;
fft2[j2 - 1] = h2r;
fft2[j2] = h2i;
fft2[n2 - j2 - 1] = h2r;
fft2[n2 - j2] = -h2i;
}
}


int cint(double x)
{
int temp;
double iprt;
if (x>0)
{
x=modf(x,&iprt);
if(fabs(x)<0.5)
temp=int(iprt);
else
temp=int(iprt+1);
}
else if(x==0)
temp=0;
else
{
x=modf(x,&iprt);
if(fabs(x)<0.5)
temp=int(iprt);
else
temp=int(iprt)-1;
}
return temp;
}

void prntft(double data[], double nn2)
{
int n,m,mm;
cout<<"n real(n) imag.(n) real(N-n) imag.(N-n)"<<endl;
cout<<setw(1)<<"0";
cout<<setw(14)<<data[1];
cout<<setw(14)<<data[2];
cout<<setw(14)<<data[1];
cout<<setw(14)<<data[2]<<endl;
for (n = 3; n<=(nn2 / 2) + 1; n=n+2)
{
m = (n - 1) / 2;
mm = nn2 + 2 - n;
cout<<setiosflags(ios::fixed);
cout<<setprecision(0)<<setw(1)<<m;
cout<<setprecision(6)<<setw(14)<<data[n];
cout<<setprecision(6)<<setw(14)<<data[n+1];
cout<<setprecision(6)<<setw(14)<<data[mm];
cout<<setprecision(6)<<setw(14)<<data[mm+1]<<endl;
}
}

void main()
{
//program d12r2
//driver for routine twofft
int n,i,n2,isign;
double data1[33], data2[33], fft1[65], fft2[65],per,x;
n = 32;
n2 = 2 * n;
per = 8.0;
const double pi = 3.1415926;
for( i = 1; i<=n; i++)
{
x = 2.0 * pi * i / per;
data1[i] = cint(cos(x));
data2[i] = cint(sin(x));
}
twofft(data1, data2, fft1, fft2, n);
cout<<setw(1)<<"Fourier transform of first function:"<<endl;
prntft(fft1, n2);
cout<<setw(1)<<"Fourier transform of second function:"<<endl;
prntft(fft2, n2);
//invert transform
isign = -1;
four1(fft1, n, isign);
cout<<setw(1)<<"Inverted transform = first function:"<<endl;
prntft(fft1, n2);
four1(fft2, n, isign);
cout<<setw(1)<<"Inverted transform = second function:"<<endl;
prntft(fft2, n2);
}