🚀 go-pugleaf

Hello World!, Signal/Image Proc. Programmers or Simply Hackers!

I ask your favor to solve a well-known problem in the signal/image 
processing fields.  I want to compute the correlation of two signals 
using the Fourier Transform.  The two discrete signals are as follows:
 
	f(n) = {1,1,1,1,2,2,2,2,8,8,8,8,9,9,9,9},    n=0,1,....,15.
	h(n) = {1,2}				     n=0,1.

The result should have the highest correlation coefficients at 1,2 and 
8,9 transitions.  In other words I want to do a pattern matching using FFT.
I was able to get the expected results in the spatial domain, but not in
the frequency domain. I haven't figured out the normalization steps to take
care of the scale changes when I used the FFT.  
Note that the lengths of the two signals are not the same, but we can 
assume that they are a power of two. I've never posted a news before.
You're more than welcome to call me if you have answers, questions and 
problems.

Thanks in advance.

Name:	Youngsup Kim 
Phone:	1-800-633-7207, 1-205-772-1999 USA
mail: 	kimy@ingr.com 			
	b15!ysk!youngsup@ingr.com		

!!HappyNewYear!!HappyNewYear!!HappyNewYear!!HappyNewYear!!HappyNewYear!!Happy
!!HappyNewYear!!HappyNewYear!!HappyNewYear!!HappyNewYear!!HappyNewYear!!Happy

==========
BACKGROUND
==========

The discrete correlation function is defined as

In the spatial domain:
----------------------

		      M-1
	f(x) o h(x) = SUM f(m) h(x+m)		for x=0,1,..,M-1.
		      m=0

The normalization should be carried out to overcome the scaling
problem.  Here is the formular to compute the correlation coefficients. 

		     SUM [f(x) - f'(x)] [h(x-m) - h']
	r(m) = ----------------------------------------------
		sqrt{ SUM [f(x) - f'(x)] SUM [h(x-m) - h'] }

	where m=0,1,...M-1, h' is the average value of the signal, 
f'(x) is the average value of f(x) in the region coincident with h(x),
and the summations are taken over the coordinates common to both f and h.
The coefficient r(m) is scaled to the range from -1 to 1, independent of
scale changes in the amplitude of f(x) and h(x).  This works great.
Once again I want to do the same thing in the frequency domain or using
the Fourier tranform. 

In the freqency domain: 
----------------------
                        -1         *
	f(x) o h(x) = FFT  [ F(n) H(n) ]

       *
where H(n) is the complex conjugate of H(n).  F(n) and H(n) are the 
Fourier Transform of f(n) and h(n), respectively.


For your convenience, I include the forward/inverse FFT routines and
complex.h.  The FFT routines are not optimized but working. 


/* 
 * File: complex.h 
 *
 *	define the complex number structure and operations
 *
 * History:
 * 11/01/89	YSK 	Creation
 * 
 */

#ifndef IM_CMPLX
#define IM_CMPLX

struct IM_sd_cmplx
 {
  float	r;	/* real part 		*/
  float	i;	/* imaginary part	*/
 };

typedef struct IM_sd_cmplx IM_S_CMPLX;
typedef struct IM_sd_cmplx *IM_p_CMPLX;


/* 
	IDEA: 
	** Computing a complex product (A + iB) * (X + iY) can be done 
        ** with only three real multiplications at a price of one additional 
        ** addition and assignment.
        ** (A+B)(X-Y)+AY-BX + i(BX+AY)
	** When you do the massive complex multiplications, you use the idea
	** described above instead of using CMPLX_mul(). 
 */

#define 	CMPLX_cnj(a)		(a.i = -a.i)
#define		CMPLX_add(a,b,c)	(a.r = b.r + c.r, a.i = b.i + c.i)
#define		CMPLX_addto(a,b)	(a.r += b.r, a.i += b.i)
#define 	CMPLX_mul(a,b,c)	((a.r = b.r * c.r - b.i * c.i), \
					 (a.i = b.r * c.i + b.i * c.r))
#define		CMPLX_abs(a,b)		(a = sqrt(b.r * b.r + b.i * b.i))
#endif






/***************************************************************************
*
*  file: fft.c
*
*	FFT()			Forward Fast Fourier Transform
*	FFT_inverse()		Inverse Fast Fourier Transform
*	is_pow_2 ()		Check if an integer is a power of 2.
*
*  Author: Youngsup Kim
*
*  History:
*
*	11/08/89	YSK	Creation 
*
****************************************************************************/

#include <stdio.h>
#include <math.h>
#include "complex.h" 

#define Debug 0

extern double sin(), cos();
	
/*-----------------------------------------------------------------------*\
				FFT
Abstract: 
 	This routine performs 1-dimensional Fast Fourier Transform.
	This implementation of the successive-doubling FFT algorithm
	is adapted from Cooley and Tukey's paper.

			1   M-1  
		F(u) = ---  SUM f(x) exp(-j2 PI(ux/M))
		        N   x=0 

	The size of the input array or the number of point to compute
	should be a power of 2.
	Also it assumes that n = 2**ln (NO checking will be done).
	The output is overwritten on the input.

Bug:	

History:
	11/1/89:	YSK	Initial implementation
\*-----------------------------------------------------------------------*/




int FFT (f, n, ln)
IM_p_CMPLX 	f;		/* input array				*/
int		n;		/* size of input array			*/
int 		ln;		/* n = 2**ln 				*/
 {
  IM_S_CMPLX 	u, w, t;
  register int	i, j;
  int		k, l, nv2, nm1, ip, le, le1;
  float		ay, bx;

  /* for (i=1, n=2; i<ln; n*=2, i++); Since n is already known, just pass it */

  nv2 = n / 2;
  nm1 = n - 1;
  j = 0;

  for (i=0; i<nm1; i++)
   {
    if (i <= j)
     {
      t = f[j];		f[j] = f[i];		f[i] = t;
     }
    k = nv2;
    while (k <= j)
     {
      j -= k;
      k /= 2;
     }
    j += k;
   }

  for (l=0; l<ln; l++)
   {

    le = 2;
    for (i=1; i<=l; le*=2, i++);
    le1 = le / 2;

    u.r = 1.0; 				u.i = 0.0;

    w.r = cos(M_PI/le1); 		w.i = -sin(M_PI/le1); 

    for (j=0; j<le1; j++)
     {
      for (i=j; i<n; i+=le)
       {
	ip = i + le1;
	CMPLX_mul (t, f[ip], u);
        CMPLX_add (f[ip], f[i], -t);
	CMPLX_addto (f[i], t);
       }

       CMPLX_mul (t, u, w); 
       u = t;

      /**
       ** Computing a complex product (A + iB) * (X + iY) can be done 
       ** with only three real multiplications at a price of one additional 
       ** addition and assignment.
       ** (A+B)(X-Y)+AY-BX + i(BX+AY)
       **/
#if 0
      ay = u.r * w.i;
      bx = u.i * w.r;
      u.r = (u.i + u.r) * (w.r - w.i) + ay - bx;
      u.i = ay + bx;
#endif
     } 
   }

  /* Normalize */
 
  for (i=0; i<n; f[i].r /= (float) n, f[i++].i /= (float) n);
 }





/*-----------------------------------------------------------------------*\
				FFT_inverse
Abstract: 
 	This routine performs 1-dimensional Inverse Fast Fourier 
	Transform in place.
	This routine simply calls the one-dimensional forward FFT after
	it conjugates and scales the input array.

	The size of the input array or the number of point to compute
	should be a power of 2.
	The input is scaled down by n.
	Also it assumes that n = 2**ln (NO checking will be done).
	The output is overwritten on the input or this is in-place
	operation.

Bug:	

History:
	11/1/89:	YSK	Initial implementation
\*-----------------------------------------------------------------------*/


int FFT_inverse (me, n, ln)
IM_p_CMPLX 	me;	/* input one-dimensional array			*/
int 		n;	/* size of "me"					*/
int 		ln;	/* n = 2 ** ln					*/ 
{
  int i;
  IM_p_CMPLX ptr;
  ptr=me;
  for (i=0; i<n; i++, ptr++) 
   {
    ptr->r =  ptr->r * n;
    ptr->i = -ptr->i * n;
   }
  return (FFT (me, n, ln));
}


/*
** is_pow_2: checks if the number is a power of 2 and returns the power,
**	or a -1 if the number isn't
*/

int is_pow_2 (k)
unsigned int	k;	/* is this number a power of 2? */
{
int		i, j, l;

for (i=1; i<sizeof(int) * 8; i++)
	{
	if (k == 0) break;

	j = 1 << i;

	if (j == k) break;
	}
if (i == (sizeof (int) * 8)) i = -1;
return (i);
}

!HappyNewYear!!HappyNewYear!!HappyNewYear!!HappyNewYear!!HappyNewYear!