Matrix user forums

Code: Select all

void DSP_FFT(MX_STRING Real, MX_UINT8 MSZ_Real, MX_STRING Imaginary, MX_UINT8 MSZ_Imaginary)
{
	MX_SINT32 mr = 0, i, j, l, k, c, istep, shift;
	MX_SINT16 qr, qi, tr, ti, wr, wi, b, m;

	m = FFT_SAMPLE_POWER;

	// decimation in time - re-order data
	for (m = 1; m <= FFT_SAMPLE_CNTM1; ++m)
	{
		l = FFT_SAMPLE_COUNT;
		do
		{
			l >>= 1;
		} while (mr+l > FFT_SAMPLE_CNTM1);

		mr = (mr & (l-1)) + l;
		if (mr <= m) continue;

		tr = Real[m];
		Real[m] = Real[mr];
		Real[mr] = tr;
		ti = Imaginary[m];
		Imaginary[m] = Imaginary[mr];
		Imaginary[mr] = ti;
	}

	l = 1;
	k = LOG2_N_WAVE-1;

	while (l < n)
	{
		istep = l << 1;
		for (m=0; m<l; ++m)
		{
			j = m << k;
			// 0 <= j < N_WAVE/2
			wr =  Sinewave[j + (FFT_N_WAVE / 4)];
			wi = -Sinewave[j];

			wr >>= 1;
			wi >>= 1;

			for (i=m; i<n; i+=istep)
			{
				j = i + l;

				// Here I unrolled the multiplications to prevent overhead
				// for procedural calls (we don't need to be clever about
				// the actual multiplications since the pic has an onboard
				// 8x8 multiplier in the ALU):

				// tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
				c = ((long int)wr * (long int)Real[j]);
				c = c >> 14;
				b = c & 0x01;
				tr = (c >> 1) + b;

				c = ((long int)wi * (long int)Imaginary[j]);
				c = c >> 14;
				b = c & 0x01;
				tr = tr - ((c >> 1) + b);

				// ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
				c = ((long int)wr * (long int)Imaginary[j]);
				c = c >> 14;
				b = c & 0x01;
				ti = (c >> 1) + b;

				c = ((long int)wi * (long int)Real[j]);
				c = c >> 14;
				b = c & 0x01;
				ti = ti + ((c >> 1) + b);

				qr = Real[i];
				qi = Imaginary[i];
				qr >>= 1;
				qi >>= 1;

				Real[j] = qr - tr;
				Imaginary[j] = qi - ti;
				Real[i] = qr + tr;
				Imaginary[i] = qi + ti;
			}
		}

		--k;
		l = istep;
	}
}