FFT c代码的收集

来源：互联网发布：什么叫菜鸟网络编辑：程序博客网时间：2024/05/18 11:17

收集的一些FFT 的c实现代码，进过整理改动能够跑通，从测试的波函数采样（采样点的幅值）经过FFT计算得出频率分布结果（该采样段的频率成分（单频波成分）以及对应频率的幅值）。下一步工作（未开始）是大量的准确度测试，效率测试和优化。

1

#include<stdio.h>#include<math.h>#include<stdlib.h>#include<unistd.h>#include<string.h>#define PI 3.1415926double x1[100000],x2[100000],w1[100000],w2[100000];//角标1代表实数部分,2虚数部分int visited[100000];int s[1000];void wnp(int M,int L,int N) //求(WN)^p{    int i,k;    int t1=(int)(pow(2.0,M-L))-1;    double t2=-2*PI/N;    double a,b;    memset(w1,0,sizeof(w1));    memset(w2,0,sizeof(w2));    k=(int)(pow(2.0,L-1));    a=cos(t2*k);    b=sin(t2*k);    w1[0]=1; //当0的情况    w2[0]=0;    for(i=1;i<=t1;i++)    {        w1[i]=a*w1[i-1]-b*w2[i-1];        w2[i]=a*w2[i-1]+b*w1[i-1];    }}void FFT(int N,int M){    int i,j,n=N;    double a,b;    for(i=1;i<=M;i++)  //从第1级到M级    {        n/=2;        memset(visited,0,sizeof(visited));//标记数组清零        wnp(M,i,N);  //调用wnp函数        for(j=0;j<N;j++)   //分别求每一级的当前实数部分和虚数部分        {            if(!visited[j])            {                visited[j]=1;                visited[j+n]=1;                a=x1[j]-x1[j+n]; //此处曾出错 应先计算出其值 因为后面 x1[j]和x2[j]的值会改变                b=x2[j]-x2[j+n];                x1[j]=x1[j]+x1[j+n];                x2[j]=x2[j]+x2[j+n];                int t=j%(N/(int)(pow(2.0,i*1.0)));                x1[j+n]=a*w1[t]-b*w2[t];                x2[j+n]=a*w2[t]+b*w1[t];            }        }    }}void solve(double *x,int N,int M) //数位倒读{    int a,i,j,k;    double t;    for(k=0;k<N;k++)    {        i=k;        a=0;        memset(s,0,sizeof(s));        for(j=0;j<M;j++)        {            s[j]=i%2;            i/=2;        }        for(j=0;j<M;j++)        {            a=a+s[j]*(int)(pow(2.0,M-1-j));        }        if(a<k)        {            t=x[a];            x[a]=x[k];            x[k]=t;        }    }}int main(){    int N,i,M;    printf("input N(N): ");    scanf("%d",&N);    if(N < 2)    {        printf("too small\n");    }    unsigned short amp[N];    int fs = 44100;    M=floor(log10(N*1.0)/log10(2.0)+0.5);    for(i=0;i<N;i++)    {        x1[i] = 7000*cos(1000*2*(double)PI*((double)i/44100) + (double)3/5*PI)                + 200*cos(600*2*(double)PI*((double)i/44100) + (double)1/2*PI);        x2[i] = 0;    }    FFT(N,M);    solve(x1,N,M);    solve(x2,N,M);    printf("得到的频谱值为:\n");    for(i=0;i<N/2;i++)    {        double originamp = sqrt(x1[i]*x1[i] + x2[i]*x2[i]);        if(i == 0)        {            amp[i] = originamp/N;        }        else        {            amp[i] = originamp/N*2;        }        if(amp[i] > 20)        {            printf("f:%.2lf amp:%d\n",(double)i*fs/N,amp[i]);        }    }    return 0;}

2

#include <math.h>#include <stdio.h>#include <stdlib.h>#include <string.h>#include <unistd.h>//精度0.0001弧度typedef struct//复数类型{    float real;       //实部    float imag;       //虚部}complex;#define PI 3.1415926///////////////////////////////////////////void conjugate_complex(int n,complex in[],complex out[]);void c_plus(complex a,complex b,complex *c);//复数加void c_mul(complex a,complex b,complex *c) ;//复数乘void c_sub(complex a,complex b,complex *c); //复数减法void c_div(complex a,complex b,complex *c); //复数除法void fft(int N,complex f[]);//傅立叶变换 输出也存在数组f中void c_abs(complex f[],float out[],int n);//复数数组取模void conjugate_complex(int n,complex in[],complex out[]){    int i = 0;    for(i=0;i<n;i++)    {        out[i].imag = -in[i].imag;        out[i].real = in[i].real;    }}void c_abs(complex f[],float out[],int n){    int i = 0;    float t;    for(i=0;i<n;i++)    {        t = f[i].real * f[i].real + f[i].imag * f[i].imag;        out[i] = sqrt(t);    }}void c_plus(complex a,complex b,complex *c){    c->real = a.real + b.real;    c->imag = a.imag + b.imag;}void c_sub(complex a,complex b,complex *c){    c->real = a.real - b.real;    c->imag = a.imag - b.imag;}void c_mul(complex a,complex b,complex *c){    c->real = a.real * b.real - a.imag * b.imag;    c->imag = a.real * b.imag + a.imag * b.real;}void c_div(complex a,complex b,complex *c){    c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag);    c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag);}#define SWAP(a,b)  tempr=(a);(a)=(b);(b)=temprvoid Wn_i(int n,int i,complex *Wn,char flag){    Wn->real = cos(2*PI*i/n);    if(flag == 1)        Wn->imag = -sin(2*PI*i/n);    else if(flag == 0)        Wn->imag = -sin(2*PI*i/n);}//傅里叶变化void fft(int N,complex f[]){    complex t,wn;//中间变量    int i,j,k,m,n,l,r,M;    int la,lb,lc;    /*----计算分解的级数M=log2(N)----*/    for(i=N,M=1;(i=i/2)!=1;M++);    /*----按照倒位序重新排列原信号----*/    for(i=1,j=N/2;i<=N-2;i++)    {        if(i<j)        {            t=f[j];            f[j]=f[i];            f[i]=t;        }        k=N/2;        while(k<=j)        {            j=j-k;            k=k/2;        }        j=j+k;    }    /*----FFT算法----*/    for(m=1;m<=M;m++)    {        la=pow(2,m); //la=2^m代表第m级每个分组所含节点数        lb=la/2;    //lb代表第m级每个分组所含碟形单元数        //同时它也表示每个碟形单元上下节点之间的距离        /*----碟形运算----*/        for(l=1;l<=lb;l++)        {            r=(l-1)*pow(2,M-m);            for(n=l-1;n<N-1;n=n+la) //遍历每个分组，分组总数为N/la            {                lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号                Wn_i(N,r,&wn,1);//wn=Wnr                c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算                c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr                c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr            }        }    }}int main(){    int N,i,M;    printf("input N(N): ");    scanf("%d",&N);    if(N < 2)    {        printf("too small\n");        return 0;    }    unsigned short amp[N];    int fs = 44100;    complex input[N];    for(i=0;i<N;i++)    {        input[i].real = 7000*cos(1000*2*(double)PI*((double)i/44100) + (double)3/5*PI)                + 200*cos(600*2*(double)PI*((double)i/44100) + (double)1/2*PI);        input[i].imag = 0;    }    fft(N,input);    printf("得到的频谱值为:\n");    for(i=0;i<N/2;i++)    {        double originamp = sqrt(input[i].real*input[i].real + input[i].imag*input[i].imag);        if(i == 0)        {            amp[i] = originamp/N;        }        else        {            amp[i] = originamp/N*2;        }        if(amp[i] > 20)        {            printf("f:%.2lf amp:%d\n",(double)i*fs/N,amp[i]);        }    }    return 0;}

0 0