FFT离散傅里叶音频分解转换_C语言实现

离散傅里叶原理讲解:

离散傅里叶变换零基础入门-中文1（针对工科生，无需连续傅立叶变换知识）_哔哩哔哩_bilibili

离散傅里叶变换零基础入门-中文2（针对工科生，无需连续傅立叶变换知识）_哔哩哔哩_bilibili

傅里叶技术画图应用新玩法:

如何用傅里叶级数绘制出任意图像来?双语_哔哩哔哩_bilibili

傅里叶一笔画图开源工程：

Fourier Series - Zoom / Mike Bostock | Observable

以下demo是通过FFT算法分析 PC录制的音频文件，输出该音频的各个频率的分布；

#include "stdio.h"
#include "math.h"

#define SR_NUM 8000 //设置采样率8000

#define pi    3.1415926
struct complex
{
	float re; //real实部
	float im; //image虚部
};

complex complexAdd(complex a, complex b) { //复数加
	complex rt;
	rt.re = a.re + b.re;
	rt.im = a.im + b.im;
	return rt;
}

complex complexMult(complex a, complex b) { //复数乘
	complex rt;
	rt.re = a.re*b.re - a.im*b.im;
	rt.im = a.im*b.re + a.re*b.im;
	return rt;
}

void complexSet(complex *a, short *b, int Num)//复数设置
{
	for (int i = 0; i < Num; i++)
	{
		a[i].re = b[i];
		a[i].im = 0;
	}
	return;
}

//离散傅里叶变换
//X[]标识变换后频域(in_out)，x[]为时域采样信号(in)
void dft(complex Xf[], complex xt[], int len)
{ 
	complex temp;
	for (int i = 0; i < len; i++)
	{
		Xf[i].re = 0;
		Xf[i].im = 0;
		for (int n = 0; n < len; n++)
		{
			temp.re = (float)cos(2 * pi*i*n/len);
			temp.im = -(float)sin(2 * pi*i*n/ len);
			Xf[i] = complexAdd(Xf[i], complexMult(xt[n], temp));
		}

		if (i >= 6000)//去除高频信号>6K
		{
			Xf[i].re = 0.0;
			Xf[i].im = 0.0;
		}


		int frqValue = int(2.0/SR_NUM*sqrt(Xf[i].re*Xf[i].re + Xf[i].im*Xf[i].im));//转换后的分频幅度
		if( frqValue>0){  //打印输出非0的幅度值
			printf("freq[%d] = %d  \n",i,frqValue );
		}


	}
}
//离散傅里叶逆变换
//Xin[]标识变换后频域，xOut[]为时域采样信号
void idft(complex Xin[], complex xOut[], int len) {
	complex temp;
	//int k, n;
	for (int k = 0; k < len; k++)
	{
		xOut[k].re = 0;
		xOut[k].im = 0;
		for (int n = 0; n < len; n++)
		{
			temp.re = (float)cos(2 * pi*k*n / len );
			temp.im = (float)sin(2 * pi*k*n / len );
			xOut[k] = complexAdd(xOut[k], complexMult(Xin[n], temp));
		}
		xOut[k].re /= len;
		xOut[k].im /= len;
	}
}


int main() {

	printf("----testing fft data ,samples rate=%d -----\n",SR_NUM );
	complex samples[SR_NUM], XFreq[SR_NUM], xTime[SR_NUM]; //原始数据，变换后的频域数据，逆变换后的时域数据
	FILE *fp;
	FILE *fp2;
	short buf[SR_NUM];
	fp=fopen("d://fft_test.pcm", "rb"); //在d盘放一个pcm文件，用cooledit等软件制作8000采样率的PCM文件
	fp2 = fopen("d://fft_out.pcm", "wb"); //在d盘新建一个pcm文件
	if (!fp ||!fp2) {
		printf("PCM file is not ready! \n");
		getchar();
		return 1;
	}
	else
	{
		while (fread(buf, sizeof(short)*SR_NUM,1,fp) > 0)//末尾数据忽略
		{
			//因为每秒采样SR_NUM个数据，所以这里每次取1秒的数据
			complexSet(samples, buf, SR_NUM);

			dft(XFreq, samples, SR_NUM);
			idft(XFreq, xTime, SR_NUM);
			for (int i = 0; i < SR_NUM; i++)
			{
				buf[i] = xTime[i].re;
			}
			fwrite(buf, sizeof(short)*SR_NUM,1 , fp2);
		}
	}
	fclose(fp);
	fclose(fp2);

	printf("-----end-----\n");
	getchar();

	return 0;
}