kimy@ingr.com (Yongsup Kim) (12/27/89)
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!JDurand@cup.portal.com (Gerald Joseph Durand) (12/29/89)
>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 Hi Kim ------ You may be helped by the following reference book "Multidimensional Digital Signal Processing" by DAN E. DUDGEON & RUSSEL M. MERSEREAU PRENTICE-HALL 1984 thank you Jerry Durand. 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 */