Xref: utzoo comp.dsp:401 comp.graphics:9082 sci.math:9091 Path: utzoo!utgpu!jarvis.csri.toronto.edu!clyde.concordia.ca!uunet!ingr!kimy From: kimy@ingr.com (Yongsup Kim) Newsgroups: comp.dsp,comp.graphics,sci.math Subject: Correlation and FFT Keywords: Correlation, FFT Message-ID: <8108@ingr.com> Date: 27 Dec 89 15:11:30 GMT Organization: Intergraph Corp. Huntsville, Al Lines: 309 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 #include #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; ir = 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