头文件：
/*
* 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
*/
/*****************************************************************************
* utilities.h
*
* Some usable routines converted from "Matlab", which are used in wavelet
* transform and time-frequency analysis, such as "filp"(same to reverse),
* "shift", "circshift", "fftshift", "dyadup", "wkeep", "wextend" and so on.
*
* Zhang Ming, 2010-01, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef UTILITIES_H
#define UTILITIES_H
#include <string>
#include <vector.h>
#include <fft.h>
namespace splab
{
int mod( int, int );
int ceil( int, int );
template<typename Type> Vector<Type> reverse( const Vector<Type>& );
template<typename Type> Vector<Type> flip( const Vector<Type>& );
template<typename Type> Vector<Type> shift( const Vector<Type>& );
template<typename Type> Vector<Type> cirshift( const Vector<Type>& );
template<typename Type> Vector<Type> fftshift( const Vector<Type>& );
template<typename Type> Vector<Type> dyadUp( const Vector<Type>&, int );
template<typename Type> Vector<Type> dyadDown( const Vector<Type>&, int );
template<typename Type> Vector<Type> fftInterp( const Vector<Type>&, int );
template<typename Type>
Vector< complex<Type> > fftInterp( const Vector< complex<Type> >&, int );
template<typename Type>
Vector<Type> wkeep( const Vector<Type>&, int, int );
template<typename Type>
Vector<Type> wkeep( const Vector<Type>&, int,
const string &direction="center" );
template<typename Type>
Vector<Type> wextend( const Vector<Type>&, int,
const string &direction="both",
const string &mode="zpd" );
#include <utilities-impl.h>
}
// namespace splab
#endif
// UTILITIES_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
*/
/*****************************************************************************
* utilities-impl.h
*
* Implementation for utilities.
*
* Zhang Ming, 2010-01 (revised 2010-08), Xi'an Jiaotong University.
*****************************************************************************/
/**
* Modulus after division. return a integer in the range of 0 to n-1.
* e.g. -1%5=-1, mod(-1,5)=4
*/
int mod( int m, int n )
{
if( n != 0 )
{
int r = m % n;
if( r < 0 )
r += n;
return r;
}
else
{
cerr << "The dividend shouldn't be zero." << endl;
return 0;
}
}
/**
* Rounds the elements of a/b to the nearest integers
* greater than or equal to a/b.
* e.g. ceil(10,2) = 5, ceil(10,3)=4.
*/
int ceil( int m, int n )
{
if( n != 0 )
{
int q = m / n;
if( m%n != 0 )
q += 1;
return q;
}
else
{
cerr << "The dividend shouldn't be zero." << endl;
return 0;
}
}
/**
* Flip vector left to right.
*/
template <typename Type>
Vector<Type> reverse( const Vector<Type> &v )
{
Vector<Type> tmp( v.size() );
typename Vector<Type>::iterator ib = tmp.begin();
typename Vector<Type>::const_iterator ie = v.end();
while( ib != tmp.end() )
*ib++ = *--ie;
return tmp;
}
template <typename Type>
inline Vector<Type> flip( const Vector<Type> &v )
{
return reverse( v );
}
/**
* Vector's shift.
*/
template <typename Type>
Vector<Type> shift( const Vector<Type> &v, int shiftsize )
{
Vector<Type> tmp( v.dim() );
if( shiftsize >= 0 )
{
typename Vector<Type>::iterator itrL = tmp.begin()+shiftsize;
typename Vector<Type>::const_iterator itrR = v.begin();
while( itrL != tmp.end() )
*itrL++ = *itrR++;
}
else
{
typename Vector<Type>::iterator itrL = tmp.begin();
typename Vector<Type>::const_iterator itrR = v.begin()-shiftsize;
while( itrR != v.end() )
*itrL++ = *itrR++;
}
return tmp;
}
/**
* Vector's circulary shift.
*/
template <typename Type>
Vector<Type> circshift( const Vector<Type> &v, int shiftsize )
{
Vector<Type> tmp( v.dim() );
if( shiftsize >= 0 )
{
typename Vector<Type>::iterator itrL = tmp.begin()+shiftsize;
typename Vector<Type>::const_iterator itrR = v.begin();
while( itrL != tmp.end() )
*itrL++ = *itrR++;
itrL = tmp.begin();
while( itrR != v.end() )
*itrL++ = *itrR++;
}
else
{
typename Vector<Type>::iterator itrL = tmp.begin();
typename Vector<Type>::const_iterator itrR = v.begin()-shiftsize;
while( itrR != v.end() )
*itrL++ = *itrR++;
itrR = v.begin();
while( itrL != tmp.end() )
*itrL++ = *itrR++;
}
return tmp;
}
/**
* Vector's fft shift.
*/
template <typename Type>
inline Vector<Type> fftshift( const Vector<Type> &v )
{
int shiftsize = v.dim() - v.dim()/2 - 1;
return circshift( v, shiftsize );
}
/**
* dyadic upsampling
* w = dyadup(v, evenodd), where v is a vector, returns an extended
* copy of vector v obtained by inserting zeros. Whether the zeros
* are inserted as even- or odd-indexed elements of w depends on the
* value of positive integer evenodd:
* If evenodd is even, then w[2k]=0, w[2k+1]=v[k], w.size()=2*v.size()+1.
* If evenodd is odd, then w[2k]=v[k], w[2k+1]=0. w.size()=2*v.size()-1.
*/
template <typename Type>
Vector<Type> dyadUp( const Vector<Type> &v, int evenodd )
{
int length = v.dim();
Vector<Type> tmp;
if( evenodd%2 == 0 )
{
tmp.resize( 2*length+1 );
for( int i=0; i<length; ++i )
{
tmp[2*i] = 0;
tmp[2*i+1] = v[i];
}
tmp[2*length] = 0;
}
else
{
tmp.resize( 2*length-1 );
for( int i=0; i<length-1; ++i )
{
tmp[2*i] = v[i];
tmp[2*i+1] = 0;
}
tmp[2*length-2] = v[length-1];
}
return tmp;
}
/**
* dyadic downsampling
* w = dyadup(v, evenodd), where v is a vector, returns a version of v
* that has been downsampled by 2. Whether w contains the even or odd
* indexed samples of v depends on the value of positive integer evenodd:
* If evenodd is even, then w[k]=v[2*k], w.size()=(v.size()+1)/2.
* If evenodd is odd, then w[k]=v[2*k+1], w.size()=v.size()/2.
*/
template <typename Type>
Vector<Type> dyadDown( const Vector<Type> &v, int evenodd )
{
int length = v.dim();
Vector<Type> tmp;
if( evenodd%2 == 0 )
{
tmp.resize( (length+1)/2 );
for( int i=0; i<tmp.dim(); ++i )
tmp[i] = v[2*i];
}
else
{
tmp.resize( length/2 );
for( int i=0; i<tmp.dim(); ++i )
tmp[i] = v[2*i+1];
}
return tmp;
}
/**
* Real signal interpolation by the method of padding zeros in frequency domain.
* The interpolation factor should be >= 1.
*/
template <typename Type>
Vector<Type> fftInterp( const Vector<Type> &sn, int factor )
{
int N = sn.size(),
halfN = N/2,
offset = (factor-1)*N;
Vector< complex<Type> > Sk = fft(sn);
Vector< complex<Type> > Xk(factor*N);
for( int i=0; i<=halfN; ++i )
Xk[i] = Type(factor)*Sk[i];
for( int i=halfN+1; i<N; ++i )
Xk[offset+i] = Type(factor)*Sk[i];
return ifftc2r(Xk);
}
/**
* Complex signal interpolation by the method of padding zeros in frequency domain.
* The interpolation factor should be >= 1.
*/
template <typename Type>
Vector< complex<Type> > fftInterp( const Vector< complex<Type> > &sn,
int factor )
{
int N = sn.size(),
halfN = N/2,
offset = (factor-1)*N;
Vector< complex<Type> > Sk = fft(sn);
Vector< complex<Type> > Xk(factor*N);
for( int i=0; i<=halfN; ++i )
Xk[i] = Type(factor)*Sk[i];
for( int i=halfN+1; i<N; ++i )
Xk[offset+i] = Type(factor)*Sk[i];
return ifft(Xk);
}
/**
* Keep part of vector.
* For a vector, w = wkeep(v,L,opt) extracts the vector w from the vector v.
* The length of w is L. If direction = "center" ("left", "rigth",
* respectively), w is the central (left, right, respectively) part of v.
* w = wkeep(x,L) is equivalent to w = wkeep(v,L,"center").
* w = wkeep(v,L,first) returns the vector v[first] to v[first+L-1].
*/
template <typename Type>
Vector<Type> wkeep( const Vector<Type> &v, int length, int first )
{
Vector<Type> tmp(length);
if( ( 0 < length ) && ( length <= v.dim()-first ) )
{
for( int i=0; i<length; ++i )
tmp[i] = v[first+i];
return tmp;
}
else
{
cerr << "Invalid length input." << endl;
return tmp;
}
}
template <typename Type>
Vector<Type> wkeep( const Vector<Type> &v, int length,
const string &direction )
{
int lv = v.dim();
Vector<Type> tmp(length);
if( ( 0 <= length ) && ( length <= lv ) )
{
if( direction == "right" )
for( int i=0; i<length; ++i )
tmp[i] = v[lv-length+i];
else if( direction == "left" )
for( int i=0; i<length; ++i )
tmp[i] = v[i];
else
{
int first = (lv-length)/2;
for( int i=0; i<length; ++i )
tmp[i] = v[first+i];
}
return tmp;
}
else
{
cerr << "Invalid length input." << endl;
return tmp;
}
}
/**
* extend vector
* The extension types are specified by the string "direction", include
* "left", "right" and "both". The default type is "both". The valid
* extension modes, which specified by strint "mode" are: zero padding
* ("zpd"), periodized extension("ppd") and symetirc extension("sym").
* The default mode is "zpd".
*/
template <typename Type>
Vector<Type> wextend( const Vector<Type> &v, int extLength,
const string &direction, const string &mode )
{
if( extLength >= 0 )
{
Vector<Type> tmp;
int lv = v.dim();
if( direction == "right" )
{
tmp.resize( lv+extLength );
for( int i=0; i<lv; ++i )
tmp[i] = v[i];
if( mode == "sym" )
for( int i=0; i<extLength; ++i )
tmp[lv+i] = v[lv-1-i];
else if( mode == "ppd" )
for( int i=0; i<extLength; ++i )
tmp[lv+i] = v[i];
else
for( int i=0; i<extLength; ++i )
tmp[lv+i] = 0;
}
else if( direction == "left" )
{
tmp.resize( lv+extLength );
if( mode == "sym" )
for( int i=0; i<extLength; ++i )
tmp[i] = v[extLength-1-i];
else if( mode == "ppd" )
for( int i=0; i<extLength; ++i )
tmp[i] = v[lv-extLength+i];
else
for( int i=0; i<extLength; ++i )
tmp[i] = 0;
for( int i=0; i<lv; ++i )
tmp[i+extLength] = v[i];
}
else
{
tmp.resize( lv+2*extLength );
for( int i=0; i<lv; ++i )
tmp[i+extLength] = v[i];
if( mode == "sym" )
for( int i=0; i<extLength; ++i )
{
tmp[i] = v[extLength-1-i];
tmp[lv+extLength+i] = v[lv-1-i];
}
else if( mode == "ppd" )
for( int i=0; i<extLength; ++i )
{
tmp[i] = v[lv-extLength+i];
tmp[lv+extLength+i] = v[i];
}
else
for( int i=0; i<extLength; ++i )
{
tmp[i] = 0;
tmp[lv+extLength+i] = 0;
}
}
return tmp;
}
else
{
cerr << "The extesion length should be greater zero." << endl;
return Vector<Type>(0);
}
}测试代码：
/*****************************************************************************
* utilities_test.cpp
*
* Utilities testing.
*
* Zhang Ming, 2010-01, Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <iomanip>
#include <string>
#include <utilities.h>
using namespace std;
using namespace splab;
const int N = 5;
int main()
{
Vector<int> v1(N);
for( int i=1; i<=v1.dim(); i++ )
v1(i) = i;
cout << "vector v1 : " << v1 << endl;
Vector<int> v2(N);
for( int i=1; i<=v2.dim(); ++i )
v2(i) = i+N;
cout << "vector v2 : " << v2 << endl;
int N = 11;
double a = 0;
double b = 1.0;
Vector<double> x = linspace( a, b, N );
cout << N << " points linearly spaced from 0 to 1.0"
<< x << endl;
cout << "Flipping vector v1 from left to right : " << flip(v1) << endl;
cout << "Shift vector v1 from left to right : " << shift(v1,2) << endl;
cout << "Shift vector v1 from right to left : " << shift(v1,-2) << endl;
cout << "Circle shift vector v1 from left to right : " << circshift(v1,2) << endl;
cout << "Circle shift vector v1 from right to left : " << circshift(v1,-2) << endl;
cout << "FFT shift of vector : " << fftshift(v1) << endl;
cout << "Dyadic upsampling of vector v1 by zeros at the even position : "
<< dyadUp( v1,0 ) << endl;
cout << "Dyadic upsampling of vector v1 by zeros at the odd position : "
<< dyadUp( v1,1 ) << endl;
cout << "Dyadic downsampling of vector v1 by zeros at the even position : "
<< dyadDown( v1,0 ) << endl;
cout << "Dyadic downsampling of vector v1 by zeros at the odd position : "
<< dyadDown( v1,1 ) << endl;
Vector<float> sn(N);
Vector< complex<float> > cn(N);
for( int i=0; i<N; ++i )
{
sn[i] = float(sin(i*TWOPI/N));
cn[i] = complex<float>( float(sin(i*TWOPI/N)), float(cos(i*TWOPI/N)) );
}
cout << setiosflags(ios::fixed) << setprecision(4);
cout << "real signal sn : " << sn << endl;
cout << "FFT interpolation of sn by factor fo 2 : "
<< fftInterp( sn, 2 ) << endl;
cout << "complex signal cn : " << cn << endl;
cout << "FFT interpolation of sn by factor fo 2 : "
<< fftInterp( cn, 2 ) << endl;
cout << resetiosflags(ios::fixed);
int n = 2;
string dire = "left";
string mode = "zpd";
cout << "Extending vector v1 in left direction by zeros padding : "
<< wextend( v1,n,dire,mode ) << endl;
mode = "ppd";
cout << "Extending vector v1 in left direction by periodic mode : "
<< wextend( v1,n,dire,mode ) << endl;
mode = "sym";
cout << "Extending vector v1 in left direction by symmetric mode : "
<< wextend( v1,n,dire,mode ) << endl;
dire = "right";
mode = "zpd";
cout << "Extending vector v1 in right direction by zeros padding : "
<< wextend( v1,n,dire,mode ) << endl;
mode = "ppd";
cout << "Extending vector v1 in right direction by periodic mode : "
<< wextend( v1,n,dire,mode ) << endl;
mode = "sym";
cout << "Extending vector v1 in right direction by symmetric mode : "
<< wextend( v1,n,dire,mode ) << endl;
dire = "both";
mode = "zpd";
cout << "Extending vector v1 in both direction by zeros padding : "
<< wextend( v1,n,dire,mode ) << endl;
mode = "ppd";
cout << "Extending vector v1 in both direction by periodic mode : "
<< wextend( v1,n,dire,mode ) << endl;
mode = "sym";
cout << "Extending vector v1 in both direction by symmetric mode : "
<< wextend( v1,n,dire,mode ) << endl;
cout << "Keeping the center part of vector v1 : " << wkeep( v1,3,"center" ) << endl;
cout << "Keeping the left part of vector v1 : " << wkeep( v1,3,"left" ) << endl;
cout << "Keeping the right part of vector v1 : " << wkeep( v1,3,"right" ) << endl;
cout << "Keeping the first(2) to first + L(3) elements of vector v1 : "
<< wkeep( v1,3,2 ) << endl;
cout << "The modulus of 2 divided by 5 is " << mod(2,5) << "." << endl;
cout << "The modulus of -1 divided by 5 is " << mod(-1,5) << "." << endl;
cout << endl;
cout << "The nearest integer >= 10/2 is " << ceil(10,2) << "." << endl;
cout << "The nearest integer >= 10/3 is " << ceil(10,3) << "." << endl;
cout << endl;
cout << "The numbers can be represented by the integer power of 2 "
<< "from 0 to 1000 are : " << endl;
return 0;
}运行结果：
vector v1 : size: 5 by 1
1
2
3
4
5
vector v2 : size: 5 by 1
6
7
8
9
10
11 points linearly spaced from 0 to 1.0size: 11 by 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Flipping vector v1 from left to right : size: 5 by 1
5
4
3
2
1
Shift vector v1 from left to right : size: 5 by 1
0
0
1
2
3
Shift vector v1 from right to left : size: 5 by 1
3
4
5
0
0
Circle shift vector v1 from left to right : size: 5 by 1
4
5
1
2
3
Circle shift vector v1 from right to left : size: 5 by 1
3
4
5
1
2
FFT shift of vector : size: 5 by 1
4
5
1
2
3
Dyadic upsampling of vector v1 by zeros at the even position : size: 11 by 1
0
1
0
2
0
3
0
4
0
5
0
Dyadic upsampling of vector v1 by zeros at the odd position : size: 9 by 1
1
0
2
0
3
0
4
0
5
Dyadic downsampling of vector v1 by zeros at the even position : size: 3 by 1
1
3
5
Dyadic downsampling of vector v1 by zeros at the odd position : size: 2 by 1
2
4
real signal sn : size: 11 by 1
0.0000
0.5406
0.9096
0.9898
0.7557
0.2817
-0.2817
-0.7557
-0.9898
-0.9096
-0.5406
FFT interpolation of sn by factor fo 2 : size: 22 by 1
0.0000
0.2817
0.5406
0.7557
0.9096
0.9898
0.9898
0.9096
0.7557
0.5406
0.2817
-0.0000
-0.2817
-0.5406
-0.7557
-0.9096
-0.9898
-0.9898
-0.9096
-0.7557
-0.5406
-0.2817
complex signal cn : size: 11 by 1
(0.0000,1.0000)
(0.5406,0.8413)
(0.9096,0.4154)
(0.9898,-0.1423)
(0.7557,-0.6549)
(0.2817,-0.9595)
(-0.2817,-0.9595)
(-0.7557,-0.6549)
(-0.9898,-0.1423)
(-0.9096,0.4154)
(-0.5406,0.8413)
FFT interpolation of sn by factor fo 2 : size: 22 by 1
(0.0000,1.0000)
(0.2817,0.9595)
(0.5406,0.8413)
(0.7557,0.6549)
(0.9096,0.4154)
(0.9898,0.1423)
(0.9898,-0.1423)
(0.9096,-0.4154)
(0.7557,-0.6549)
(0.5406,-0.8413)
(0.2817,-0.9595)
(0.0000,-1.0000)
(-0.2817,-0.9595)
(-0.5406,-0.8413)
(-0.7557,-0.6549)
(-0.9096,-0.4154)
(-0.9898,-0.1423)
(-0.9898,0.1423)
(-0.9096,0.4154)
(-0.7558,0.6549)
(-0.5406,0.8413)
(-0.2817,0.9595)
Extending vector v1 in left direction by zeros padding : size: 7 by 1
0
0
1
2
3
4
5
Extending vector v1 in left direction by periodic mode : size: 7 by 1
4
5
1
2
3
4
5
Extending vector v1 in left direction by symmetric mode : size: 7 by 1
2
1
1
2
3
4
5
Extending vector v1 in right direction by zeros padding : size: 7 by 1
1
2
3
4
5
0
0
Extending vector v1 in right direction by periodic mode : size: 7 by 1
1
2
3
4
5
1
2
Extending vector v1 in right direction by symmetric mode : size: 7 by 1
1
2
3
4
5
5
4
Extending vector v1 in both direction by zeros padding : size: 9 by 1
0
0
1
2
3
4
5
0
0
Extending vector v1 in both direction by periodic mode : size: 9 by 1
4
5
1
2
3
4
5
1
2
Extending vector v1 in both direction by symmetric mode : size: 9 by 1
2
1
1
2
3
4
5
5
4
Keeping the center part of vector v1 : size: 3 by 1
2
3
4
Keeping the left part of vector v1 : size: 3 by 1
1
2
3
Keeping the right part of vector v1 : size: 3 by 1
3
4
5
Keeping the first(2) to first + L(3) elements of vector v1 : size: 3 by 1
3
4
5
The modulus of 2 divided by 5 is 2.
The modulus of -1 divided by 5 is 4.
The nearest integer >= 10/2 is 5.
The nearest integer >= 10/3 is 4.
The numbers can be represented by the integer power of 2 from 0 to 1000 are :
Process returned 0 (0x0) execution time : 0.234 s
Press any key to continue.
扫一扫
366
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
