快速傅里叶变换(fft)及其逆变换(iff)的c代码实现

Posted 2020-12-21 zyy-summary

#define float sample_t

// data的长度为n,必须是2的指数倍，result的长度为2n，其中奇数项保存实数，偶数项保存的是虚数
int fft(sample_t *data, int sample_number, sample_t *result)
{
	// 需要给奇数部分填充虚数0
	for(int i = 0; i < sample_number; ++i)
	{
		result[2*i] = data[i];
		result[2*i+1] = 0;
	}
	int flag = 1;
	flag = fft_ifft_implement(result, sample_number, flag);
	return flag;
}

// data的长度是2n，result的长度为n,n必须是2的指数倍
int ifft(sample_t *data, int sample_number, sample_t *result)
{
	int flag = -1;
	flag = fft_ifft_implement(data, sample_number, flag);
	// 奇数部分是虚数，需要舍弃
	for (int i = 0; i < sample_number; i++)
	{
		result[i] = data[2*i] / sample_number;
	}
	return flag;

}

static int fft_ifft_implement(sample_t *data, int sample_number, int flag)
{
	// 判断样本个数是不是2的指数倍，如果不是能否补零成指数倍？
	sample_t number_log = log(sample_number) / log(2);
	int mmax = 2, j=0;
	int n = sample_number<<1;
	int istep,  m;
	sample_t theta, wtemp, wpr, wpi, wr, wi, tempr, tempi;
	if (((int)number_log - number_log) != 0)
	{
		return 0;
	}
	for(int i = 0; i < n-1; i=i+2)
	{
		if(j > i)
		{
			swap(data, j ,i);
			swap(data, j + 1 ,i + 1);
		}
		m = n / 2;
		while(m >= 2 && j >= m)
		{
			j = j - m;
			m = m / 2;
		}
		j = j + m;
	}
	while(n > mmax)
	{
		istep = mmax<<1;
		theta = -2 * PI / (flag * mmax);
		wtemp = sin(0.5 * theta);
		wpr = -2.0 * wtemp * wtemp;
		wpi = sin(theta);
		wr = 1.0;
		wi = 0.0;
		for(int m = 1; m < mmax; m = m + 2)
		{
			for(int i = m; i < n + 1; i = i + istep)
			{
				int j = i + mmax;
				tempr = wr * data[j-1] - wi * data[j];
				tempi = wr * data[j] + wi * data[j-1];
				data[j-1] = data[i-1] - tempr;
				data[j] = data[i] - tempi;
				data[i-1] += tempr;
				data[i] += tempi;
			}
			wtemp = wr;
			wr += wr * wpr - wi * wpi;
			wi += wi * wpr + wtemp * wpi;
		}
		mmax = istep;
	}
	return 1;
}

static void swap(sample_t *data ,int m, int n)
{
	sample_t temp = data[m];
	data[m] = data[n];
	data[n] = temp;
}