Reworking FFT code

[twobob]

[root@kindle sound]# ./ffttest-arm

Testing a 64-point FFT.
1:1 mapping; no reorder required.

SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 8.90 us per 64-point FFT, or 14.39 Mbps with QPSK,
or 30.75 insns per point on a 220MHz SA-1100.

Yeah baby. happy days.

Ill post a how-to once I test it out

[twobob]

Also noted for the future:

http://www.nr.com/public-domain.html

Public Domain Software

In general, the Numerical Recipes routines are NOT in the public domain, and are NOT available for use without a valid license. However, a few small parts of Numerical Recipes have been put into the public domain. These can be freely copied and redistributed:


Utility files for C:

* nrutil.c
* nrutil.h
* complex.c
* complex.h


Wavelet routines in C:

* wt1.c
* daub4.c
* pwt.c
* pwtset.c
* wtn.c

[twobob]

Hmm.

bash-4.1$ arm-none-linux-gnueabi-gcc -O3 -march=armv6j -mtune=arm1136jf-s -mfpu=vfp -fomit-frame-pointer -o ffttest-arm radix4fft.c testfft.c testmain.c -lm

Any way to squeeze a bit more performance out of that in terms of compilation?
*Goes and googles things.*

Well -03 from http://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html looks to be as good as it gets - I vaguely recall someone posting a nice set of args but that could have been for the K5.

Cheers.

[geekmaster]

Originally Posted by twobob

Hmm.

bash-4.1$ arm-none-linux-gnueabi-gcc -O3 -march=armv6j -mtune=arm1136jf-s -mfpu=vfp -fomit-frame-pointer -o ffttest-arm radix4fft.c testfft.c testmain.c -lm

Any way to squeeze a bit more performance out of that in terms of compilation?
*Goes and googles things.*

Did you make changes to any of those .c files?

[twobob]

Originally Posted by geekmaster

Did you make changes to any of those .c files?

Nope siree. ran as it.

Did some reading for the processor. read the readme.

Started to make a MAKE file and then thought the better of it and just compiled on command line.

Just like that.

[twobob]

Okay so some perfomance feedback.

Using my current Load test program

[root@kindle sound]# ./Balsamic noise.wav

Spoiler:

#include <sched.h>
#include <errno.h>
#include <getopt.h>
#include <asoundlib.h>
#include <math.h>
#include <pthread.h>
// HEAD WELD

//====================================================
// gmplay 1.5a - geekmaster's kindle video player
// Copyright (C) 2012 by geekmaster, with MIT license:
// http://www.opensource.org/licenses/mit-license.php
//----------------------------------------------------
// Tested on DX,DXG,K3,K4main,K4diags,K5main,K5diags.
//----------------------------------------------------
#include <sys/ioctl.h>     // ioctl
#include <sys/mman.h>     // mmap, munmap
#include <stdio.h>       // printf
#include <stdlib.h>     // malloc, free
#include <linux/fb.h>  // screeninfo
#include <sys/time.h> // gettimeofday
#include <unistd.h>  // usleep
#include <string.h> // memset, memcpy
#include <fcntl.h> // open, close, write
#include <time.h> // time



#include <math.h>
#include <memory.h>

// THIS VERSION JUST FEEDS THE VIDEO IN FROM THE AUDIO FOR LOAD TESTING.

// USE A LONG AUDIO FILE. THE SEGFAULT AT THE END OF THE AUDIO FILE IS EXPECTED...
// AS IS THE NASTY SCREETCHY NOISE :)


typedef unsigned long u64;
typedef unsigned int u32;
typedef unsigned short u16;
typedef unsigned char u8;
u32 __invalid_size_argument_for_IOC; // ioctl.h bug fix for tcc
//----- eink definitions from eink_fb.h and mxcfb.h -----
typedef signed char prog_int8_t;
#define EU3 0x46dd
#define EU50 0x4040462e
#define EU51 0x4048462e
#define NUM_THREADS 3
#define VOLUME_BOUND 100
struct update_area_t {int x1,y1,x2,y2,which_fx;u8 *buffer;};
struct mxcfb_rect {u32 top,left,width,height;};
struct mxcfb_alt_buffer_data {u32 phys_addr,width,height;
struct mxcfb_rect alt_update_region;};
struct mxcfb_update_data {struct mxcfb_rect update_region;
u32 waveform_mode,update_mode,update_marker;int temp;uint flags;
struct mxcfb_alt_buffer_data alt_buffer_data;};
struct mxcfb_update_data51 {struct mxcfb_rect update_region;
u32 waveform_mode,update_mode,update_marker;
u32 hist_bw_waveform_mode,hist_gray_waveform_mode;
int temp;uint flags;struct mxcfb_alt_buffer_data alt_buffer_data;};
//----- function prototypes -----
void gmplay4(void);
//void gmplay8(void);
int getmsec(void);
int gmlib(int);
//----- gmlib global vars -----
enum GMLIB_op {GMLIB_INIT,GMLIB_CLOSE,GMLIB_UPDATE,GMLIB_VSYNC};

typedef enum {
	AUDIO_VOLUME_SET,
	AUDIO_VOLUME_GET,
} audio_volume_action;

u8 *fb0=NULL;     // framebuffer pointer
int fdFB=0;      // fb0 file descriptor
int teu=0;      // eink update time
u32 fs=0;      // fb0 stride
u32 MX=0;     // xres (visible)
u32 MY=0;    // yres (visible)
u32 VY=0;   // (VY>MY): mxcfb driver
u8 ppb=0;  // pixels per byte
u32 fc=0; // frame counter
#define FBSIZE (600/8*800)

//====================================
// gmlib - geekmaster function library
// op (init, update, vsync, close)
//------------------------------------
int gmlib(int op) {
	static struct update_area_t ua={0,0,600,800,21,NULL};
	static struct mxcfb_update_data ur={
			{0,0,600,800},257,0,1,0x1001,0,{0,0,0,{0,0,0,0}}};
	static struct mxcfb_update_data51 ur51={
			{0,0,600,800},257,0,1,0,0,0x1001,0,{0,0,0,{0,0,0,0}}};
	static int eupcode; static void *eupdata=NULL;
	struct fb_var_screeninfo screeninfo;
	if (GMLIB_INIT==op) { teu=getmsec(); fdFB=open("/dev/fb0",O_RDWR);
	ioctl(fdFB,FBIOGET_VSCREENINFO,&screeninfo);
	ppb=8/screeninfo.bits_per_pixel; fs=screeninfo.xres_virtual/ppb;
	VY=screeninfo.yres_virtual; MX=screeninfo.xres; MY=screeninfo.yres;
	ua.x2=MX; ua.y2=MY; ur.update_region.width=MX; ur.update_region.height=MY;
	fb0=(u8 *)mmap(0,MY*fs,PROT_READ|PROT_WRITE,MAP_SHARED,fdFB,0); // map fb0
	if (VY>MY) { eupcode=EU50; eupdata=&ur; ur.update_mode=0;
	if (ioctl(fdFB,eupcode,eupdata)<0) { eupcode=EU51; eupdata=&ur51; }
	} else { eupcode=EU3; eupdata=&ua; }
	system("eips -f -c;eips -c"); sleep(1);
	} else if (GMLIB_UPDATE==op) {
		if (ioctl(fdFB,eupcode,eupdata)<0) system("eips ''");  // 5.1.0 fallback
	} else if (GMLIB_VSYNC==op) { while (teu>getmsec()) usleep(1000); // fb0 busy
	} else if (GMLIB_CLOSE==op) { gmlib(GMLIB_UPDATE); sleep(1); // last screen
	system("eips -f -c;eips -c"); munmap(fb0,MY*fs); close(fdFB);
	} else { return -1; }
	return 0;
}
//====================================
// getmsec - get msec since first call
// (tick counter wraps every 12 days)
//------------------------------------
int getmsec(void) {
	int tc; static int ts=0; struct timeval tv;
	gettimeofday(&tv,NULL); tc=tv.tv_usec/1000+1000*(0xFFFFF&tv.tv_sec);
	if (0==ts) ts=tc;
	return tc-ts;
}

// Points to loaded WAVE file's data
unsigned char			*WavePtr;

// Size (in frames) of loaded WAVE file's data
snd_pcm_uframes_t		WaveSize;

// Sample rate
unsigned short			WaveRate;

// Bit resolution
unsigned char			WaveBits;

// Number of channels in the wave file
unsigned char			WaveChannels;

// For WAVE file loading
static const unsigned char Riff[4]	= { 'R', 'I', 'F', 'F' };
static const unsigned char Wave[4] = { 'W', 'A', 'V', 'E' };
static const unsigned char Fmt[4] = { 'f', 'm', 't', ' ' };
static const unsigned char Data[4] = { 'd', 'a', 't', 'a' };

snd_pcm_t *playback_handle;

unsigned int i,j;
int playcount;
int indexer;
int err;





#pragma pack (1)
/////////////////////// WAVE File Stuff /////////////////////
// An IFF file header looks like this
typedef struct _FILE_head
{
	unsigned char	ID[4];
	// could be {'R', 'I', 'F', 'F'} or {'F', 'O', 'R', 'M'}
	unsigned int	Length;
	// Length of subsequent file (including remainder of header).
	//This is in Intel reverse byte order if RIFF, Motorola format if FORM.
	unsigned char	Type[4];
	// {'W', 'A', 'V', 'E'} or {'A', 'I', 'F', 'F'}
} FILE_head;


// An IFF chunk header looks like this
typedef struct _CHUNK_head
{
	unsigned char ID[4];	// 4 ascii chars that is the chunk ID
	unsigned int	Length;	// Length of subsequent data within this chunk. This is in Intel reverse byte
	// order if RIFF, Motorola format if FORM. Note: this doesn't include any
	// extra byte needed to pad the chunk out to an even size.
} CHUNK_head;

// WAVE fmt chunk
typedef struct _FORMAT {
	short				wFormatTag;
	unsigned short	wChannels;
	unsigned int	dwSamplesPerSec;
	unsigned int	dwAvgBytesPerSec;
	unsigned short	wBlockAlign;
	unsigned short	wBitsPerSample;
	// Note: there may be additional fields here, depending upon wFormatTag
} FORMAT;
#pragma pack()

/********************** compareID() *********************
 * Compares the passed ID str (ie, a ptr to 4 Ascii
 * bytes) with the ID at the passed ptr. Returns TRUE if
 * a match, FALSE if not.
 */

static unsigned char compareID(const unsigned char * id, unsigned char * ptr)
{
	register unsigned char i = 4;

	while (i--)
	{
		if ( *(id)++ != *(ptr)++ ) return(0);
	}
	return(1);
}




/********************** waveLoad() *********************
 * Loads a WAVE file.
 *
 * fn =			Filename to load.
 *
 * RETURNS: 0 if success, non-zero if not.
 *
 * NOTE: Sets the global "WavePtr" to an allocated buffer
 * containing the wave data, and "WaveSize" to the size
 * in sample points.
 */

// END WELD

FILE_head				head;

FORMAT	format;

static unsigned char waveLoad(const char *fn)
{
	const char				*message;
	FILE_head				head;
	register int			inHandle;

	if ((inHandle = open(fn, O_RDONLY)) == -1)
		message = "didn't open";

	// Read in IFF File header
	else
	{
		if (read(inHandle, &head, sizeof(FILE_head)) == sizeof(FILE_head))
		{
			// Is it a RIFF and WAVE?
			if (!compareID(&Riff[0], &head.ID[0]) || !compareID(&Wave[0], &head.Type[0]))
			{
				message = "is not a WAVE file";
				goto bad;
			}

			// Read in next chunk header
			while (read(inHandle, &head, sizeof(CHUNK_head)) == sizeof(CHUNK_head))
			{
				// ============================ Is it a fmt chunk? ===============================
				if (compareID(&Fmt[0], &head.ID[0]))
				{


					// Read in the remainder of chunk
					if (read(inHandle, &format.wFormatTag, sizeof(FORMAT)) != sizeof(FORMAT)) break;

					// Can't handle compressed WAVE files
					if (format.wFormatTag != 1)
					{
						message = "compressed WAVE not supported";
						goto bad;
					}

					WaveBits = (unsigned char)format.wBitsPerSample;
					WaveRate = (unsigned short)format.dwSamplesPerSec;
					WaveChannels = format.wChannels;
				}

				// ============================ Is it a data chunk? ===============================
				else if (compareID(&Data[0], &head.ID[0]))
				{
					// Size of wave data is head.Length. Allocate a buffer and read in the wave data
					if (!(WavePtr = (unsigned char *)malloc(head.Length)))
					{
						message = "won't fit in RAM";
						goto bad;
					}

					if (read(inHandle, WavePtr, head.Length) != head.Length)
					{
						free(WavePtr);
						break;
					}

					// Store size (in frames)
					WaveSize = (head.Length * 8) / ((unsigned int)WaveBits * (unsigned int)WaveChannels);

					//


					close(inHandle);
					return(0);
				}

				// ============================ Skip this chunk ===============================
				else
				{
					if (head.Length & 1) ++head.Length;  // If odd, round it up to account for pad byte
					lseek(inHandle, head.Length, SEEK_CUR);
				}
			}
		}

		message = "is a bad WAVE file";
		bad:	close(inHandle);
	}

	printf("%s %s\n", fn, message);
	return(1);
}


/*********************** free_wave_data() *********************
 * Frees any wave data we loaded.
 *
 * NOTE: A pointer to the wave data be in the global
 * "WavePtr".
 */

static void free_wave_data(void)
{
	if (WavePtr) free(WavePtr);
	WavePtr = 0;
}



float highpass(float* sig, float freq, float *del, int vecsize, float sr){
	double costh, coef; costh = 2. - cos(2*M_PI*freq/sr);
	coef = costh - sqrt(costh*costh - 1.);

	for(i =0; i < vecsize; i++){
		sig[i] = (float) (sig[i]*(1 - coef) - *del*coef);
		*del = sig[i];
	}

	return *sig;
}

// Do this last.

int here;
int frameCount;
int totalCount;
//==================================
int holderx, holdery;

// gmplay4 - play video on 4-bit fb0
//----------------------------------
void gmplay4(void) {

	// We should scoop info from the unified file - and possibly via an indirect buffer.
	// Time indexing lookup - perhaps via a mutex provsioned RO value might be an idea
	// They are both TODO

	frameCount = 0;
	u32 i,x,y,b,p,off=(MY/2-400)*fs+MX/4-150,fbsize=FBSIZE; u8 fbt[FBSIZE];


	while (1)
	{


		int pure = (WavePtr[frameCount+64]) +1 / 8; // Don't div by 0

	//	printf ("DVAL: %f",&delayed);
	//	printf (" PURE: %u\n",pure);



		if (frameCount %11==0)
		{
			for (y=0;y<800;y++) for (x=0;x<600;x+=(8+(frameCount%pure))) {

				b=WavePtr[frameCount];  // "Random" 0-254 val;

				i=y*fs+x/2+off;


				// Paint it black.

				fb0[i] = fb0[i+1] = fb0[i+2] = fb0[i+3] = b;


			} fc++; gmlib(GMLIB_UPDATE);
		}

		else


		if (frameCount %4==0)
		{
			for (y=0;y<800;y++) for (x=0;x<600;x+=(8+(frameCount%4))) {

				// Or something

				//  FFT RESULTS(THIS)

		//		b=WavePtr[((frameCount)+600/8*y+x/(8))%(WaveSize-1)]; i=y*fs+x/2+off;

				b=WavePtr[((frameCount)+600/(8))*y+(x/8)%(WaveSize-1)]; i=y*fs+x/2+off;

				//	  b=result; i=y*fs+x/2+off;
				p=(b&1)*240; b>>=1; fb0[i]=p|(b&1)*15; b>>=1;
				p=(b&1)*240; b>>=1; fb0[i+1]=p|(b&1)*15; b>>=1;
				p=(b&1)*240; b>>=1; fb0[i+2]=p|(b&1)*15; b>>=1;
				p=(b&1)*240; b>>=1; fb0[i+3]=p|(b&1)*15;
			} fc++; gmlib(GMLIB_UPDATE);
		}

		else
		{
			usleep(10);
		//	system("eips -f -c;eips -c");
		}
		frameCount++;
	}
}



int
playback_callback (snd_pcm_sframes_t nframes)
{


	playcount++;
	if (WaveBits < 15 )
	{
		indexer = 4096*(playcount);
	}
	else
	{
		indexer = 4096*(playcount * 2 );
	}

	if ((err = snd_pcm_writei (playback_handle, WavePtr+indexer, nframes)) < 0) {
		fprintf (stderr, "write failed (%s)\n", snd_strerror (err));
	}
	return err;

}


// Showtime

main (int argc, char *argv[])
{

	int i;
	gmlib(GMLIB_INIT);

	char *wavename = argv[1];

	if (!waveLoad(wavename))
	{
		register int		err;
	}
	else
	{

		printf("Loaded: %s\n", wavename);

	}

	snd_pcm_hw_params_t *hw_params;
	snd_pcm_sw_params_t *sw_params;
	snd_pcm_sframes_t frames_to_deliver;
	int nfds;
	int err;
	struct pollfd *pfds;

	if ((err = snd_pcm_open (&playback_handle, "plughw", SND_PCM_STREAM_PLAYBACK, 0)) < 0) {
		fprintf (stderr, "cannot open audio device plughw (%s)\n",
				snd_strerror (err));
		exit (1);
	}

// There are far better ways of doing this...
	if (NULL==argv[2])
	{
		if ((err = snd_pcm_set_params(playback_handle,
				SND_PCM_FORMAT_S16_LE,
				SND_PCM_ACCESS_RW_INTERLEAVED,
				1,
				WaveRate,
				0,
				900000)) < 0) {   /* 0.5sec */
			printf("Playback open error: %s\n", snd_strerror(err));
			exit(EXIT_FAILURE);
		}
	}
	else
	{
		if ((err = snd_pcm_set_params(playback_handle,
				SND_PCM_FORMAT_S16_LE,
				SND_PCM_ACCESS_RW_INTERLEAVED,
				1,
				atoi(argv[2]),
				0,
				900000)) < 0) {   /* 0.5sec */
			printf("Playback open error: %s\n", snd_strerror(err));
			exit(EXIT_FAILURE);
		}
	}


	/* tell ALSA to wake us up whenever 4096 or more frames
		   of playback data can be delivered. Also, tell
		   ALSA that we'll start the device ourselves.
	 */

	if ((err = snd_pcm_sw_params_malloc (&sw_params)) < 0) {
		fprintf (stderr, "cannot allocate software parameters structure (%s)\n",
				snd_strerror (err));
		exit (1);
	}
	if ((err = snd_pcm_sw_params_current (playback_handle, sw_params)) < 0) {
		fprintf (stderr, "cannot initialize software parameters structure (%s)\n",
				snd_strerror (err));
		exit (1);
	}
	if ((err = snd_pcm_sw_params_set_avail_min (playback_handle, sw_params, 4096)) < 0) {
		fprintf (stderr, "cannot set minimum available count (%s)\n",
				snd_strerror (err));
		exit (1);
	}
	if ((err = snd_pcm_sw_params_set_start_threshold (playback_handle, sw_params, 0U)) < 0) {
		fprintf (stderr, "cannot set start mode (%s)\n",
				snd_strerror (err));
		exit (1);
	}
	if ((err = snd_pcm_sw_params (playback_handle, sw_params)) < 0) {
		fprintf (stderr, "cannot set software parameters (%s)\n",
				snd_strerror (err));
		exit (1);
	}

	/* the interface will interrupt the kernel every 4096 frames, and ALSA
		   will wake up this program very soon after that.
	 */

	if ((err = snd_pcm_prepare (playback_handle)) < 0) {
		fprintf (stderr, "cannot prepare audio interface for use (%s)\n",
				snd_strerror (err));
		exit (1);
	}


	// SETUP THE SCREEN

	pthread_t threads[NUM_THREADS];
	int rc;
	long t;

	gmlib(GMLIB_INIT);

	// Make a thread


	printf("Balsamic Main: creating Video thread %ld\n", t);
	printf("Wave size %d\n",WaveSize);
	char *message1 = "Video Thread";

	// This is the thread.

	rc = pthread_create(&threads[t], NULL, gmplay4, (void*) message1);

	if (rc){
		printf("ERROR; return code from pthread_create() is %d\n", rc);
		exit(-1);
		i=getmsec()/100; printf("%d frames in %0.1f secs = %2.1f FPS\n",
				fc,(double)i/10.0,(double)fc*10.0/i);
		gmlib(GMLIB_CLOSE);

	}



	while (1) {

		//  printf("bits %d\n",WaveBits);

		// HERE IS OUR UPDATE LOOP SEED

		// We should create a buffer area to scoop audio into from the main feed -
		// as audio is bloody picky about actually playing
		// it would be nice to consider some type of time indexing for this I guess to give us a sense of info unity.


		/* wait till the interface is ready for data, or 1 second
			   has elapsed.
		 */

		if ((err = snd_pcm_wait (playback_handle, 1000)) < 0) {
			fprintf (stderr, "poll failed (%s)\n", strerror (errno));
			break;
		}

		/* find out how much space is available for playback data */

		if ((frames_to_deliver = snd_pcm_avail_update (playback_handle)) < 0) {
			if (frames_to_deliver == -EPIPE) {
				fprintf (stderr, "an xrun occured\n");
				break;
			} else {
				fprintf (stderr, "unknown ALSA avail update return value (%d)\n",
						frames_to_deliver);
				break;
			}
		}

		if (	frames_to_deliver > 4096 ){ frames_to_deliver = 4096;}

		/* deliver the data */

		if (playback_callback (frames_to_deliver) != frames_to_deliver) {
			fprintf (stderr, "playback callback failed\n");
			break;
		}
	}

	//	close(inHandle);

	i=getmsec()/100; printf("%d frames in %0.1f secs = %2.1f FPS\n",
			fc,(double)i/10.0,(double)fc*10.0/i);
	gmlib(GMLIB_CLOSE);
	snd_pcm_close (playback_handle);
	exit (0);
}

I get results:

Before and after load applied during results


Ill format in a bit

Spoiler:

SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed...^[[A 8.91 us per 64-point FFT, or 14.37 Mbps with QPSK,
or 30.79 insns per point on a 220MHz SA-1100.
[root@kindle sound]# ./ffttest-arm 
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 8.89 us per 64-point FFT, or 14.40 Mbps with QPSK,
or 30.71 insns per point on a 220MHz SA-1100.
[root@kindle sound]# ./ffttest-arm 
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 8.89 us per 64-point FFT, or 14.40 Mbps with QPSK,
or 30.73 insns per point on a 220MHz SA-1100.
[root@kindle sound]# ./ffttest-arm 
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 8.89 us per 64-point FFT, or 14.39 Mbps with QPSK,
or 30.74 insns per point on a 220MHz SA-1100.
[root@kindle sound]# ./ffttest-arm 
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 8.91 us per 64-point FFT, or 14.37 Mbps with QPSK,
or 30.79 insns per point on a 220MHz SA-1100.

LOAD APPLIED TO 100% CPU USE FROM OTHER SHELL

[root@kindle sound]# ./ffttest-arm 
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 19.66 us per 64-point FFT, or 6.51 Mbps with QPSK,
or 67.96 insns per point on a 220MHz SA-1100.
[root@kindle sound]# ./ffttest-arm 
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 19.43 us per 64-point FFT, or 6.59 Mbps with QPSK,
or 67.14 insns per point on a 220MHz SA-1100.
[root@kindle sound]# ./ffttest-arm 
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 20.85 us per 64-point FFT, or 6.14 Mbps with QPSK,
or 72.07 insns per point on a 220MHz SA-1100.
[root@kindle sound]# ./ffttest-arm 
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 19.64 us per 64-point FFT, or 6.52 Mbps with QPSK,
or 67.88 insns per point on a 220MHz SA-1100.
[root@kindle sound]# ./ffttest-arm 
Testing a 64-point FFT.
1:1 mapping; no reorder required.
SNR: 14167370.493333 (71.512893 dB)
Mean energy in: 0.989857 out_test: 0.989554 out_ref: 0.989579
Timing FFT speed... 20.81 us per 64-point FFT, or 6.15 Mbps with QPSK,
or 71.91 insns per point on a 220MHz SA-1100.

I thought it would be worse...

[geekmaster]

Now you need to DISPLAY that FFT on the screen. Where is the spectrographic waterfall display?

[knc1]

Originally Posted by twobob

Hmm.

bash-4.1$ arm-none-linux-gnueabi-gcc -O3 -march=armv6j -mtune=arm1136jf-s -mfpu=vfp -fomit-frame-pointer -o ffttest-arm radix4fft.c testfft.c testmain.c -lm

Any way to squeeze a bit more performance out of that in terms of compilation?
*Goes and googles things.*

Well -03 from http://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html looks to be as good as it gets - I vaguely recall someone posting a nice set of args but that could have been for the K5.

Cheers.

Someone did. I did and they are indexed by the site index.
The ones above should be a "one size fits all" set, executable by all.

Question: Are you running that under the emulator I posted or on the real hardware?

Question: GCC does not generate mixed 32/16 bit code for arm by default, have you tried enabling mixed code generation (thumb)?
Note: the K3 only does thumb, the K4&K5 does thumb-2

Question: Do you want me to post the armv7l (K4/K5) emulator along with the armv6l that is already posted?

[knc1]

Originally Posted by twobob

In the interim I am looking at this. http://www.lartmaker.nl/projects/fft-arm/

I'm Happy with the smallest size - we just want some results.

the code is 11 years old. I'm hoping maths hasn't moved on too much in that time.

That code spends a significant amount of its time generating the look-up tables when doing such a small fft sample.

Extract the code to generate the lookup tables, write them (the tables) to disk.
Then just mmap them in when needed since they are constant data.

If that helps, then let me know, I will lookup the objcopy magic-foo to turn the binary tables into an elf-object that you can link into the code section of the program.
Bad idea, seems that bin2c is the currently recommended way.

****** Viral Warning *******
That code carries the GPL virus. The other code I suggested is BSD licensed.

[knc1]

Attached, source and static-arm build of bin2c.

If I type any more in this post, the post will be longer than the source code.

Build notes:

(armv6l:1) /home # gcc --static --static-libgcc -O2 -o bin2c bin2c.c
(armv6l:1) /home # ldd bin2c
	not a dynamic executable
(armv6l:1) /home # strip --strip-unneeded bin2c
(armv6l:1) /home # tar -czf bin2c-static-arm.tar.gz bin2c bin2c.c

[geekmaster]

Regarding gcc optimization, I read a report awhile back about an automated test that tried to find the best combination of gcc options to get the best optimization (by analyzing code build time, code size, and code speed), and the surprising results were that MOST of the optimization options did absolutely NOTHING (no changes to generated code). In the end after all that study, they recommended just using -O3, as I recall.

Of course, that was for x86 code, so gcc code generation for other architectures may vary. That test was also version-specific, and things may have changed since then. The main point was to not expect much performance gain from most of the optimization settings available in gcc (in that particular test environment).

Also, the ARM FPU is not noted for speed, so integer math may outperform it in some cases. Vector processing looks interesting though (and there are nice co-processors in the i.MX508 used in the K4/K5 that could help as well).

[knc1]

Originally Posted by geekmaster

Of course, that was for x86 code, so gcc code generation for other architectures may vary.

Under the specific combinations of conditions that apply to the above post -O2 produces the smallest code.
Both -O3 and -Os where larger.
Might have been different if I fiddled with the options, but did not feel the cost/benefit was worth the effort on something a person uses twice a decade.

[twobob]

I concur about the GPL. I'll see if I can grab/implement something else then. It was as much the actual implementation of something that works (rather than theoretical ports) that attracted me. Looking at something end - to - end may give me inspiration / education - to get the other ports done.

I have no less than 10 different - not - quite - working - implementations to sift through and analyse - this was a nice solution to shove some numbers in for now.

If the GPL becomes an issue then sobeit. Worst ways I will contact the author and see what shakes after a decade.
*notes drool on keyboard* - note to self - remember to go to bed before passing out.
Thanks for the feedback guys. and the tool. Much appreciated. bedtime

[knc1]

Originally Posted by twobob

I concur about the GPL. I'll see if I can grab/implement something else then. It was as much the actual implementation of something that works (rather than theoretical ports) that attracted me. Looking at something end - to - end may give me inspiration / education - to get the other ports done.

I have no less than 10 different - not - quite - working - implementations to sift through and analyse - this was a nice solution to shove some numbers in for now.

If the GPL becomes an issue then sobeit. Worst ways I will contact the author and see what shakes after a decade.
*notes drool on keyboard* - note to self - remember to go to bed before passing out.
Thanks for the feedback guys. and the tool. Much appreciated. bedtime

I think all of the ones we have looked at today are the same butterfly algorithm.
One working should make the others easier to understand.

We don't want to lose track that the original Ardino code was for both fft and ifft.
If we are going to have that 64 band equalizer on our Kindles, we will need both.
Or many other special effects that can be done on the fft before being put back into the time domain with the ifft.
Like voice disguising, so we can make crank calls over VoIP on our e-book.

One of the (late in the day) references mentions that ifft only differs in the sign of e - meaning the fft only needs one more lookup table to get an ifft out of it.

When I wake up tomorrow, will see if I can spot that reference again.

[twobob]

Originally Posted by knc1

One of the (late in the day) references mentions that ifft only differs in the sign of e - meaning the fft only needs one more lookup table to get an ifft out of it.

When I wake up tomorrow, will see if I can spot that reference again.

FFT ARM 0.0.1 - README

"Want an inverse FFT ? If you look closely at the formula for the FFT,
you'll see that the only difference between an FFT and its inverse is
the sign in the exponent of e. This means that *by reordering the output
alone* you can turn an FFT into its inverse. Again, it's been over three
years since I wrote this code, but have a look at TransTable in testmain.c
for more hints. My gut feeling is that you need the "In-place IFFT table".

thanks M. Yeah that was pretty much my thought process - look for the gaps to get the more relevant code happening.

The very dark room is made considerably brighter by even the dimmest of lights.
Have a great night