FFT離散傅裡葉音頻分解轉換_C語言實作

離散傅裡葉原理講解:

https://www.bilibili.com/video/av44600709

https://www.bilibili.com/video/BV1g4411a7rb/?spm_id_from=333.788.videocard.5

傅裡葉技術畫圖應用新玩法:

https://www.bilibili.com/video/av49238862 

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