FFT - 基于STM32 DSP扩充傅里叶变换数量上限

STM32F4\H7等自带DSP内核的系列，其DSP提供了离散快速傅里叶变换的算法，但是其最高上限数据一般是4096，即可以最大分析4096个点的频谱，用在实时性频谱分析完全是足够的，但是如果一次性对很多的数据如64K、128K、512K的数据分析，一般采用kissfft库计算也比较方便，但没有发挥STM32 DSP库的优势。 离散FFT的蝶形算法可以网上查找。

1. 计算数据需要为4096的倍数，下面计算倍数的取2次幂指数。比如8096 -- 1， 65536 --4

int8_t vCheckNPT2Mul(uint32_t NPT, uint32_t base) {
	uint8_t timesOne = 0;	// the bit 1 times
	uint8_t timesOneLocation = 0;
	uint8_t	bit = 0;
	uint32_t	tmp;
	for(uint32_t j = 0; j < 32; j++) {
		tmp = (base >> j) & 0x01;
		if( tmp == 1 ) {
			timesOne++;
			bit = j;
		}
		if( timesOne > 1) {
			return -1;
		}
	}
	// if 0, bit is 0
	timesOne = 0;
	for(uint32_t i = 0; i < 32; i++) {
		tmp = (NPT >> i) & 0x01;
		if( tmp == 1) {
			timesOne++;
			timesOneLocation = i;
		}
		if(timesOne > 0 && i < bit ) {
			return -1;	// NPT is not multiple of 4096
		}
		if(timesOne > 1) {
			return -1;	// NPT is not multiple of 2.
		}
	}
	return timesOneLocation - bit;
}

 2. 蝶形算法的Rader算法，但是这里个数只有上面的倍数计算值, 目前最高支持点数1M，可以根据需要扩充。

// Max suppoert 1M
void radersort(uint32_t space, uint8_t *pLocation)
{
	uint8_t	i=0, j=0, k=0, tmp=0;
	if(space >= 256) {
		printf("MAX support Point is: 1048576\r\n");
		return;
	}
	for(uint8_t i = 0; i < space; i++) {
		pLocation[i] = i;	// fill the 0, 1, 2, 3, 4....
	}
	for (i=0;i<space-1;i++) {
		if(i<j) {
			tmp=pLocation[i];
			pLocation[i]=pLocation[j];
			pLocation[j]=tmp;
		}
		k=space/2;
		while(j>=k) {
			j=j-k;
			k=k/2;
		}
		j=j+k;
	}
}

 3. 离散傅里叶变换

void dsp_fft_bigN(uint32_t NPT, complex_t *pSrc, complex_t *pOut)
{
	// add to big.
	uint32_t	base = 4096;
	complex_t	tmp0, tmp1;
	float	tmpsin, tmpcos;
	uint8_t	count;
	uint8_t multiple2 = vCheckNPT2Mul(NPT, base);
	uint8_t	Space	= (0x01 << multiple2);
	float 	Wn;
	if( multiple2 < 0) {
		printf("[ERROR] number must be multiple of 4096.\r\n");
		return;
	}
	// rearrange the point
	uint8_t	sHead[256];
	radersort(Space, sHead);
	for(uint8_t j = 0; j < Space; j++) {
		for(uint32_t i = 0; i < NPT/Space; i++) {
			pOut[i + j*NPT/Space] = pSrc[sHead[j] + Space*i];
		}
	}
	// get the Space number of 4096 series.
	for(uint8_t i = 0; i < Space; i++) {
		// need modify here.
		arm_cfft_f32(&arm_cfft_sR_f32_len4096, (float *)&pOut[i*base], 0, 1);
	}
	// butterfly algorithm to get the result.
	// first exp calculate is difficult, need dsp.
	//Wnk: cos(2*PI*k/N) - j * sin(2*PI*k/N)
	count = Space;
	while (count != 1) {
		for (uint8_t m = 0; m < count; m += 2) {
			m = m;
			for (uint32_t k = m * base; k < m*base + base; k++) {
				tmp0 = pOut[k];
				tmp1 = pOut[base + k];
				Wn = 2 * PI * k / (2 * base);
				tmpsin = arm_sin_f32(Wn);
				tmpcos = arm_cos_f32(Wn);
				pOut[k].r = tmp0.r + (tmpcos * tmp1.r)
						+ (tmpsin * tmp1.i);
				pOut[k].i = tmp0.i + (tmpcos * tmp1.i)
						- (tmpsin * tmp1.r);
				pOut[base + k].r = tmp0.r - (tmpcos * tmp1.r)
						- (tmpsin * tmp1.i);
				pOut[base + k].i = tmp0.i - (tmpcos * tmp1.i)
						+ (tmpsin * tmp1.r);
			}
		}
		base = 2 * base; 		// iterate until finally.
		count = count / 2; 	// Space
	}
}

一切根据蝶形算法改编。