powerspectol - FFTW@wiki

/**
* Powerspectol calculation
* @author maeda
* @date 2008/09/01
*/
#include "powerspectol.h"
#define _USE_MATH_DEFINES
#include "math.h"
 
/**
* get the log scaled powerspectol of the signal and normalize it to 0-1 range
* @param in input data ( frequency domain )
* @param out output data ( frequency domain )
* @param cols ( width )
* @param rows ( height )
* @return true if succeed
* -------------------------------------------------------------------------
* (1) get the log scaled powerspectol : log(1 + sqrt(re*re + im*im))
* (2) normalize
*/
int maGetPowerSpectol2D( fftw_complex *in, double *out, int cols, int rows )
{
	int i,j;
	int idx; // index of data
	double max, min, scale; // max/min of powerspectol
 
	if( in==NULL || out==NULL ) return false;
	if( rows<0 || cols<0 )      return false;
 
	// (1) get the log scaled powerspectol : log(1 + sqrt(re*re + im*im))
	//     ( powerspectol is defined as re*re + im*im )
	// this is equal or greater than 0
	for( j=0; j<rows; j++ ){
		for( i=0; i<cols; i++ ){
			idx = j*cols + i;
			out[idx] = log10( 1 + sqrt(pow(in[idx][0],2) + pow(in[idx][1],2)) );
		}
	}
 
	// (2) normalization
 
	// (2-1) get the maximum and minimum
	max = min = out[0];
 
	for( j=0; j<rows; j++ ){
		for( i=0; i<cols; i++ ){
			idx = j*cols + i;
			if( out[idx]>max ) max = out[idx];
			if( out[idx]<min ) min = out[idx];
		}
	}
 
	// (2-2) normalize
	scale = 1./(max-min);
	for( j=0; j<rows; j++ ){
		for( i=0; i<cols; i++ ){
			idx = j*cols + i;
			out[idx] = ( out[idx] - min ) * scale;
		}
	}
 
	return true;
}