频谱分析：c和python对比FFT的效率并画出幅度谱

一、c语言编写的FFT程序

c语言程序

正弦波表达式为： s(t) = 0.6 sin( 2π 50t ) 和s(t) = 0.6 sin( 2π 500t )
频率为8000Hz,近似为8192
FFT程序来源于徐士良老师算法书中
1.FFT.c

#include "stdio.h"
#include "math.h"
#include "G:\C code\FFTDuiBi\FFTDuiBi\kfft.c"//注意路径
#include <time.h>

#define PI 3.1415926535

main()
{
	int i, j;
	double pr[8192], pi[8192], fr[8192], fi[8192], t[8192];
	clock_t begin, end;
	double cost1, cost2, persent;
	for (i = 0; i <8192; i++)  //生成输入信号
	{
		t[i] = i*0.001;
		pr[i] = 0.6*sin(2 * PI * 500 * i) + 0.6*sin(2 * PI * 50 * i); pi[i] = 0.0; //0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i)
	}
	begin = clock();	
	//begin=time(NULL);								//开始记录时间
	kfft(pr, pi, 8192, 13, fr, fi);  //调用FFT函数

	end = clock();										//结束记录时间
	cost1 = (double)(end - begin ) / CLOCKS_PER_SEC;
	//cost1 = (double)(end - begin);
	for (i = 0; i<8192; i++)
	{
		printf("%d\t%lf\n", i, pr[i]); //输出结果
	}
	printf("花费时间为：%lfs",cost1);

	
}

2.kfft.c

#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;
}

编译结果并绘图

在TCC下进行编译，并用gnuplot绘图

 plot "<FFT.exe" w l

二、python编写

python程序

import numpy as np#导入一个数据处理模块
import time

import matplotlib.pyplot as plt#导入一个绘图模块

# 依据快速傅里叶算法得到信号的频域
def test_fft():
    sampling_rate = 8192  # 采样率
    fft_size = 8192  # FFT取样长度
    t = np.arange(0, 8.192, 1.0 / sampling_rate)
    #np.arange(起点，终点，间隔)产生8.192s长的取样时间
    x=0.6*np.sin(2*np.pi*500*t)+0.6*np.sin(2*np.pi*50*t)
    # 两个正弦波叠加，500HZ和50HZ
    # N点FFT进行精确频谱分析的要求是N个取样点包含整数个取样对象的波形。
    # 因此N点FFT能够完美计算频谱对取样对象的要求是n*Fs/N（n*采样频率/FFT长度），
    # 因此对8KHZ和512点而言，完美采样对象的周期最小要求是8000/512=15.625HZ,
    # 所以156.25的n为10,234.375的n为15。

    t1=time.perf_counter()
    xs = x[:fft_size]# 从波形数据中取样fft_size个点进行运算

    xf = np.fft.rfft(xs) / fft_size  # 返回fft_size/2+1 个频率

    t2 = time.perf_counter()
    #利用np.fft.rfft()进行FFT计算，rfft()是为了更方便对实数信号进行变换，
    # 由公式可知 / fft_size为了正确显示波形能量
    # rfft函数的返回值是N/2+1个复数，分别表示从0(Hz)到sampling_rate/2(Hz)的分。

    # 于是可以通过下面的np.linspace计算出返回值中每个下标对应的真正的频率：
    freqs = np.linspace(0, sampling_rate*10, fft_size/2+1 )  # 表示频率
    #freqs = np.linspace(0, sampling_rate/2 , fft_size/2  + 1)  # 表示频率

    xfp = 20 * np.log10(np.clip(np.abs(xf), 1e-20, 1e100))

    #xfp = np.abs(xf)   # 代表信号的幅值，即振幅
    # 最后我们计算每个频率分量的幅值，并通过 20*np.log10()将其转换为以db单位的值。
    # 为了防止0幅值的成分造成log10无法计算，我们调用np.clip对xf的幅值进行上下限处理

    plt.figure(figsize=(8, 4))
    plt.subplot(211)
    plt.plot(t[:fft_size], xs)
    plt.xlabel(u"时间(秒)", fontproperties='FangSong')
    plt.title(u"500Hz和50Hz的波形和频谱", fontproperties='FangSong')

    plt.subplot(212)
    plt.plot(freqs, xfp)
    plt.xlabel(u"频率(Hz)", fontproperties='FangSong')
    #字体FangSong
    plt.ylabel(u'幅值', fontproperties='FangSong')
    plt.subplots_adjust(hspace=0.4)
    '''subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)
    有六个可选参数来控制子图布局。值均为0~1之间。其中left、bottom、right、top围成的区域就是子图的区域。
    wspace、hspace分别表示子图之间左右、上下的间距。实际的默认值由matplotlibrc文件控制的。
    '''
    plt.show()
    print(f"python的FFT花费的时间为：{t2-t1}")



test_fft()

结果显示：

花费时间对比

python运行FFT花费时间：

python运行FFT花费时间为：3.1e-6

c运行FFT花费时间：

c运行FFT花费时间为：0.002000s