FFT c代码的收集

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

1

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<unistd.h>
#include<string.h>

#define PI 3.1415926

double 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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

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)=tempr

void 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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174