离散Gabor变换算法的C++实现

2019独角兽企业重金招聘Python工程师标准>>>

头文件：

/*
 * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 2 or any later version.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * This program is distributed in the hope that it will 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 is available at:
 * http://www.fsf.org/licensing/licenses
 */


/*****************************************************************************
 *                                    dgt.h
 *
 * Discrete Gabor Transform.
 *
 * These routines are designed for calculating discrete Gabor transform and
 * its inversion of 1D signals. In order to eliminate the border effect, the
 * input signal("signal") is extended by three forms: zeros padded("zpd"),
 * periodized extension("ppd") and symetric extension("sym").
 *
 * The analysis/synthesis function is given by users, and it's daul
 * (synthesis/analysis) function can be computed by "daul" routine. The over
 * sampling rate is equal to N/dM, where N denotes frequency sampling numbers
 * and dM denotes the time sampling interval.
 *
 * N and dM should can be devided evenly by the window length "Lw". The
 * recovered signal just has the elements from 1 to dM*floor(Ls/dM) of the
 * original signal. So you'd better let dM can be deviede evenly by  the
 * original signal length "Ls".
 *
 * Zhang Ming, 2010-03, Xi'an Jiaotong University.
 *****************************************************************************/


#ifndef DGT_H
#define DGT_H


#include <string>
#include <fft.h>
#include <matrix.h>
#include <linequs1.h>
//#include <linequs3.h>
#include <utilities.h>


namespace splab
{


    template<typename Type>
    Vector<Type> daul( const Vector<Type>&, int, int );

    template<typename Type>
    Matrix< complex<Type> > dgt( const Vector<Type>&,
                                 const Vector<Type>&,
                                 int, int,
                                 const string &mode = "zpd" );
    template<typename Type>
    Vector<Type> idgt( const Matrix< complex<Type> >&,
                       const Vector<Type>&,
                       int, int );


    #include <dgt-impl.h>

}
// namespace splab


#endif
// DGT_H

实现文件：

/*
 * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 2 or any later version.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * This program is distributed in the hope that it will 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 is available at:
 * http://www.fsf.org/licensing/licenses
 */


/*****************************************************************************
 *                                dgt-impl.h
 *
 * Implementation for Discrete Gabor Transform.
 *
 * Zhang Ming, 2010-03, Xi'an Jiaotong University.
 *****************************************************************************/


/**
 * Return the daul function of input window "g".
 */
template <typename Type>
Vector<Type> daul( const Vector<Type> &gn, int N, int dM )
{
    int L = gn.size(),
        NL = 2*L-N;

    Vector<Type> hn(L);
    Vector<Type> gg = wextend( gn, NL, "right", "zpd" );

    // coefficient matrix and constant vector
    Matrix<Type> H(NL/N,L/dM);
    Vector<Type> u(NL/N);
    u[0] = Type(1.0/N);

    // evaluate matrix H
    for( int k=0; k<dM; ++k )
    {
        for( int q=0; q<NL/N; ++q )
            for( int p=0; p<L/dM; ++p )
            {
                int index = mod( k+p*dM+q*N, NL );
                H[q][p] = gg[index];
            }

        // calculate the kth part value of h
        Vector<Type> tmp = trMult( H, luSolver( multTr(H,H), u ) );
//        Vector<Type> tmp = tsvd( H, u );
        for( int i=0; i<tmp.dim(); ++i )
            hn[k+i*dM] = tmp[i];
    }
    return hn;
}


/**
 * Compute the discrete Gabor transform of "signal". The coeffitions are
 * stored in "coefs", a Ls+Lw by Lw (Ls: lengthof "signsl", Lw:
 * length of "window") matrix. The row represents frequency ordinate
 * and column represents the time ordinate.
 */
template <typename Type>
Matrix< complex<Type> > dgt( const Vector<Type> &signal,
                             const Vector<Type> &anaWin,
                             int N, int dM, const string &mode )
{
    int Ls = signal.dim();
    int Lw = anaWin.dim();
    int M = (Ls+Lw)/dM;

    Vector<Type> sn = wextend( signal, Lw, "both", mode );

    Matrix< complex<Type> > coefs(N,M);
    Vector<Type> segment(Lw);
    Vector< complex<Type> > segDFT(Lw);
    Vector< complex<Type> > tmp(N);
    complex<Type> W = polar( Type(1), Type(-2*PI/N) );

    for( int m=0; m<M; ++m )
    {
        // intercept signal by window function
        for( int i=0; i<Lw; ++i )
            segment[i] = sn[i+m*dM]*anaWin[i];

        // Fourier transform
        segDFT = fft( segment );

        // calculate the mth culumn coefficients
        for( int n=0; n<N; ++n )
            tmp[n] = pow(W,n*m*dM) * segDFT[n*Lw/N];

        coefs.setColumn( tmp, m );
    }

    return coefs;
}


/**
 * Compute the inverse discrete Gabor transform from "coefs".
 */
template <typename Type>
Vector<Type> idgt( const Matrix< complex<Type> > &coefs,
                   const Vector<Type> &synWin,
                   int N, int dM )
{
    int M = coefs.cols();
    int Lw = synWin.size();
    int Ls = dM*M-Lw;

    // reallocate for signal and initialize it by "0"
    Vector<Type> signal(Ls);

    Matrix<Type> idftCoefs(N,M);
    Vector<Type> sn(N);
    Vector< complex<Type> > Sk(N);

    for( int i=0; i<M; ++i )
    {
        Sk = coefs.getColumn(i);
        sn = ifftc2r(Sk);
        sn *= Type(N);

        idftCoefs.setColumn( sn, i );
    }

    // compulate the ith element of signal
    for( int i=0; i<Ls; ++i )
    {
        int p = ceil(i+1, dM);
        for( int m=p; m<p+Lw/dM; ++m )
        {
            int n = mod(i,N);
            signal[i] += idftCoefs[n][m]*synWin[Lw-m*dM+i];
        }
    }

    return signal;
}

测试代码：

/*****************************************************************************
 *                               dgt_test.cpp
 *
 * Discrete Gabor transform testing.
 *
 * Zhang Ming, 2010-03, Xi'an Jiaotong University.
 *****************************************************************************/


#define BOUNDS_CHECK

#include <iostream>
#include <vectormath.h>
#include <timing.h>
#include <dgt.h>


using namespace std;
using namespace splab;


typedef double  Type;
const   int     Fs = 1000;
const   int     Ls = 10000;
const   int     Lg = 80;
const   int     N  = 40;
const   int     dM = 10;      // over sampling ratio is N/dM


int main()
{

	/******************************* [ signal ] ******************************/
	Type a = 0;
	Type b = Ls-1;
	Vector<Type> t = linspace( a, b, Ls ) / Type(Fs);
	Vector<Type> st = cos( Type(400*PI) * pow(t,Type(2.0)) );

	/******************************* [ widow ] ******************************/
	a = 0;
	b = Type( Lg-1 );
	Type r = sqrt( dM*N / Type(2*PI) );
	Type u = (Lg-1) / Type(2.0);
	t = linspace( a, b, Lg );
	Vector<Type> h = gauss( t, u, r );
	h = h / norm(h);

	/***************************** [ daul function ] *************************/
	Type runtime = 0.0;
	Timing cnt;
	cout << "Compute daul function." << endl;
	cnt.start();
	Vector<Type> g = daul( h, N, dM );
	cnt.stop();
	runtime = cnt.read();
	cout << "The running time = " << runtime << " (ms)" << endl << endl;

	/******************************** [ DGT ] ********************************/
	cout << "Taking discrete Gabor transform." << endl;
	cnt.start();
	Matrix< complex<Type> > C = dgt( st, g, N, dM, "sym" );
	cnt.stop();
	runtime = cnt.read();
	cout << "The running time = " << runtime << " (ms)" << endl << endl;

	/******************************** [ IDGT ] *******************************/
	cout << "Taking inverse discrete Gabor transform." << endl;
	cnt.start();
	Vector<Type> xt = idgt( C, h, N, dM );
	cnt.stop();
	runtime = cnt.read();
	cout << "The running time = " << runtime << " (ms)" << endl << endl;

	cout <<"The relative error : norm(s-x) / norm(s) = "
	     << norm(st-xt)/norm(st) << endl << endl;

	return 0;
}

运行结果：

Compute daul function.
The running time = 8.67362e-019 (ms)

Taking discrete Gabor transform.
The running time = 0.063 (ms)

Taking inverse discrete Gabor transform.
The running time = 0.016 (ms)

The relative error : norm(s-x) / norm(s) = 5.79937e-012


Process returned 0 (0x0)   execution time : 0.172 s
Press any key to continue.

转载于:https://my.oschina.net/zmjerry/blog/7708