C + + implementation of the basic---convolution and its fast algorithm for image processing

Header file:

/* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), [email protected] * * This program is a free software; You can redistribute it and/or modify it * under the terms of the GNU general public License as published by the * Free so Ftware Foundation, either version 2 or any later version. * * Redistribution and use in source and binary forms, with or without * modification, is permitted provided that the FOL Lowing conditions is met: * * 1. Redistributions of source code must retain the above copyright notice, * This list of conditions and the following disc Laimer. * * 2. Redistributions in binary form must reproduce the above copyright * Notice, this list of conditions and the following D Isclaimer in the * documentation and/or other materials provided with the distribution. * * This program was distributed in the hope that it'll be useful, but without * any WARRANTY; Without even the implied warranty of merchantability or * FITNESS for A particular PURPOSE. See the GNU GeneRAL Public License for * more details. A copy of the GNU general public License are available at: * http://www.fsf.org/licensing/licenses *//********************* * Convolution.h * * Linear Convol Ution and polynomial multiplication. * * The convolution routine "conv" is implemented by it's definition in time * domain. If the sequence to being convoluted is a long, you should use the * fast convolution algorithm "Fastconv", which is Implemente d in frequency * Domain by usin FFT. * * Zhang Ming, 2010-01, Xi ' an Jiaotong University. /#ifndef Convolution_h#define Convolution_h#include <vector.h> #include <fft.h> #include <utilities.h>namespace splab{template                                               <typename type> vector<type> Conv (const vector<type>&   Const vector<type>&); Template<typename type> vector<type> convolution (const vector<type>&    Const vector<type>&);                                                   Template<typename type> vector<type> fastconv (const vector<type>&    Const vector<type>&); #include <convolution-impl.h>}//namespace splab#endif//convolution_h

Implementation file:

/* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), [email protected] * * This program is a free software; You can redistribute it and/or modify it * under the terms of the GNU general public License as published by the * Free so Ftware Foundation, either version 2 or any later version. * * Redistribution and use in source and binary forms, with or without * modification, is permitted provided that the FOL Lowing conditions is met: * * 1. Redistributions of source code must retain the above copyright notice, * This list of conditions and the following disc Laimer. * * 2. Redistributions in binary form must reproduce the above copyright * Notice, this list of conditions and the following D Isclaimer in the * documentation and/or other materials provided with the distribution. * * This program was distributed in the hope that it'll be useful, but without * any WARRANTY; Without even the implied warranty of merchantability or * FITNESS for A particular PURPOSE. See the GNU GeneRAL Public License for * more details. A copy of the GNU general public License are available at: * http://www.fsf.org/licensing/licenses *//********************* * convolution-impl.h * * Implementati On for linear convolution. * * Zhang Ming, 2010-01, Xi ' an Jiaotong University. * * Convolution and ploynonal Multiplication. */template <typename type>vector<type> Conv (const vector<type> &signal, const VECTOR&LT;TYPE    > &filter) {if (Signal.dim () < Filter.dim ()) Return convolution (filter, signal); else return convolution (signal, filter);} Template <typename type>vector<type> convolution (const vector<type> &signal, const Vector<    Type> &filter) {int siglength = Signal.dim ();    int fillength = Filter.dim (); ASSERT (Siglength >= FillenGTH);    int length = siglength + filLength-1;    vector<type> x (length);        for (int i=1; i<=length; ++i) {x (i) = 0;        if (I < fillength) for (int j=1; j<=i; ++j) x (i) + = Filter (j) * Signal (I-J+1); else if (i <= siglength) for (int j=1; j<=fillength; ++j) x (i) + = Filter (j) * Signal (I-J        +1);    else for (int j=i-siglength+1; j<=fillength; ++j) x (i) + = Filter (j) * Signal (I-J+1); } return x;} /** * Fast convolution by FFT. */template<typename type>vector<type> fastconv (const vector<type> &AMP;XN, const VECTOR&LT;TYPE    > &yn) {int M = Xn.dim (), N = Yn.dim (); vector<type> xnpadded = wextend (xn, N-1, "right", "ZPD"), ynpadded = Wextend (yn, M-1, "right", "Z    PD "); Return IFFTC2R (FFT (xnpadded) * FFT (ynpadded));//vector< complex<type> > Zk = FFT (xnpadded) * FFT (Ynpadde D);/return IFFTC2R (ZK);//Return IFFTC2R (FFT (Wextend (xn,n-1, "right", "ZPD")) * FFT (Wextend (yn,m-1, "right", "ZPD")));} 

Test code:

/***************************************************************************** * Convolution.cpp * * convolution testing. * * Zhang Ming, 2010-01, Xi ' an Jiaotong University. /#define Bounds_check#include <  Iostream> #include <convolution.h>using namespace std;using namespace splab;typedef double Type;const int M    = 3;const int N = 5;int Main () {vector<type> xn (M), yn (n);    Vector<type> Zn;    for (int i=0; i<m; ++i) xn[i] = i;    for (int i=0; i<n; ++i) yn[i] = I-N/2;    Convolution zn = conv (xn, yn);    cout << "xn:" << xn << endl << "yn:" << yn << endl;    cout << "Convolution of xn and yn:" << Zn << Endl;    Zn = Fastconv (xn, yn);    cout << "Fast convolution of xn and yn:" << Zn << Endl; return 0;} 

Operation Result:

Xn:  size:3 by 1012yn:  size:5 by 1-2-1012convolution of Xn and yn:   size:7 by 10-2-5-2144fast convolution of Xn and yn:   size:7 by 1-2.53765e-016-2-5-2144process returned 0 (0x0)   execution time:0.078 spress any key to con Tinue.

http://my.oschina.net/zmjerry/blog/3671

C + + implementation of the basic---convolution and its fast algorithm for image processing