FFT快速傅里叶变换C语言实现

fft.c 测试文件（改进自某篇文章，但找不到了）

#include "stdio.h"
#include "math.h"
#include "kfft.c"
#include <time.h>

#define PI 3.1415926535

int main()
{ 
	int i,j;
    int Fs = 1024;  //采样速率
    int L = 32768;  //采样点个数
    int K = 15;     //2^K=L
    double pr[L],pi[L],fr[L],fi[L],t[L];
	clock_t begin, end;
	double cost1 ,cost2, persent;
    for (i=0; i<L; i++)  //生成输入信号
    { 
		t[i] = (double)i/Fs;
		pr[i]=0.6*sin(2*PI*50*t[i])+sin(2*PI*10*t[i]); 
        pi[i]=0.0; //0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i)
	}
	begin = clock();									//开始记录时间
    kfft(pr,pi,L,K,fr,fi);  //调用FFT函数
	end = clock();										//结束记录时间
	cost1 = (double)(end - begin) / CLOCKS_PER_SEC;
	// for (i=0; i<L/2+1; i++)
    // { 
    //     printf("%d\t%f\n", Fs/L*i, pr[i]*2/L); //输出结果
    // }
    printf("%f", cost1);
    return 0;
}

kfft.c fft函数源文件

  #include "math.h"
  void kfft(pr,pi,n,k,fr,fi)
  int n,k;
  double pr[],pi[],fr[],fi[];
  { 
	int it,m,is,i,j,nv,l0;
    double p,q,s,vr,vi,poddr,poddi;
    for (it=0; it<=n-1; it++)  //将pr的实部和虚部循环赋值给fr[]和fi[]
    { 
		m=it; 
		is=0;
		for(i=0; i<=k-1; i++)
        { 
			j=m/2; 
			is=2*is+(m-2*j); 
			m=j;
		}
        fr[it]=pr[is]; 
        fi[it]=pi[is];
    }
    pr[0]=1.0; 
    pi[0]=0.0;
    p=6.283185306/(1.0*n);
    pr[1]=cos(p); //将w=e^-j2pi/n用欧拉公式表示
    pi[1]=-sin(p);

    for (i=2; i<=n-1; i++)  //计算pr[]
    { 
		p=pr[i-1]*pr[1]; 
		q=pi[i-1]*pi[1];
		s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
		pr[i]=p-q; pi[i]=s-p-q;
    }
    for (it=0; it<=n-2; it=it+2)  
    { 
		vr=fr[it]; 
		vi=fi[it];
		fr[it]=vr+fr[it+1]; 
		fi[it]=vi+fi[it+1];
		fr[it+1]=vr-fr[it+1]; 
		fi[it+1]=vi-fi[it+1];
    }
	m=n/2; 
	nv=2;
    for (l0=k-2; l0>=0; l0--) //蝴蝶操作
    { 
		m=m/2; 
		nv=2*nv;
        for (it=0; it<=(m-1)*nv; it=it+nv)
          for (j=0; j<=(nv/2)-1; j++)
            { 
				p=pr[m*j]*fr[it+j+nv/2];
				q=pi[m*j]*fi[it+j+nv/2];
				s=pr[m*j]+pi[m*j];
				s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
				poddr=p-q; 
				poddi=s-p-q;
				fr[it+j+nv/2]=fr[it+j]-poddr;
				fi[it+j+nv/2]=fi[it+j]-poddi;
				fr[it+j]=fr[it+j]+poddr;
				fi[it+j]=fi[it+j]+poddi;
            }
    }
    for (i=0; i<n; i++)
       { 
		  pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);  //幅度的计算
       }
    return;
  }

函数说明-FROTRAN常用算法程序集（徐士良）：

 

matlab对比代码：

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector
%50Hz和120Hz频率信号
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
%添加随机噪声
X = S + 2*randn(size(t));
figure;plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')
% %加噪声后信号X长度1500，进行1500点FFT
t1 = clock;
Y = fft(X);
t2 = clock;
% %为何除以信号长度L?
P2 = abs(Y/L);
P1 = P2(1:L/2+1);%单边谱
P1(2:end-1) = 2*P1(2:end-1);%由于P1(1)是直流吧
f = Fs*(0:(L/2))/L;%采样频率Fs,因此只看fs/2内的信号
figure;plot(f, P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')
disp("time:" + etime(t2, t1))
% % %原始信号S长度1500，进行1500点FFT
% t1 = clock;
% Y = fft(S);
% Y1 = abs(Y);
% Y2 = Y1(1:L/2+1);
% t2 = clock;
% P2 = abs(Y/L);
% P1 = P2(1:L/2+1);
% P1(2:end-1) = 2*P1(2:end-1);
% f = Fs*(0:(L/2))/L;%采样频率Fs,因此只看fs/2内的信号
% figure;plot(f,Y2) 
% title('Single-Sided Amplitude Spectrum of S(t)')
% xlabel('f (Hz)')
% ylabel('|P1(f)|')
% disp("time:" + etime(t2, t1))