最近在研究信号处理的C语言算法,突然发现一个西安交大的师兄之前已经写出来了一个完整的库,同样是研究生三年,差别怎么这样大呢。
先从用轮子开始吧。
一、SP++3.0安装及测试
官网下载地址:
https://code.google.com/archive/p/tspl/downloads
将 SP++3.0 压缩包解压到某一路径下,如 D:\Program Files\SP++3.0
建立VS2010 项目,在Project->Propertiesr的VC++Directories->Include Directories
中加入D:\Program Files\SP++3.0\include
然后添加一个fft_test.cpp,进行测试
/*****************************************************************************
* fft_test.cpp
*
* FFT test.
*
* Zhang Ming, 2010-09, Xi'an Jiaotong University.
*****************************************************************************/
#include <iostream>
#include <cstdlib>
#include <vectormath.h>
#include <fft.h>
using namespace std;
using namespace splab;
typedef double Type;
const int MINLEN = 1;
const int MAXLEN = 1000;
const int STEP = 10;
int main()
{
Vector< complex<Type> > sn, Rk, Sk, xn;
Vector<Type> rn, tn;
cout << "forward transform: complex to complex." << endl;
cout << "inverse transform: complex to complex." << endl << endl;
cout << "signal length" << "\t" << "mean(abs((sn-xn))" << endl;
for( int len=MINLEN; len<MAXLEN; len+=STEP )
{
sn.resize(len);
for( int i=0; i<len; ++i )
sn[i] = complex<Type>( rand()%10, rand()%10 );
Sk = fftc2c( sn );
xn = ifftc2c( Sk );
// Sk = fft( sn );
// xn = ifft( Sk );
cout << " " << len << "\t\t" << " " << sum(abs(sn-xn))/len << endl;
}
cout << endl << endl;
cout << "forward transform: real to complex ." << endl;
cout << "inverse transform: complex to real." << endl << endl;
cout << "signal length" << "\t" << "mean(abs((rn-tn))" << endl;
for( int len=MINLEN; len<MAXLEN; len+=STEP )
{
rn.resize(len);
for( int i=0; i<len; ++i )
rn[i] = rand()%10;
Rk = fftr2c( rn );
tn = ifftc2r( Rk );
// Rk = fft( rn );
// tn = real( ifft(Rk) );
cout << " " << len << "\t\t" << " " << sum(abs(rn-tn))/len << endl;
}
cout << endl;
return 0;
}
编译会产生错误:min,max出现在未命名空间,直接注释即可(也可在主文件添加#inlucde<algorthm>)。
之后编译成功。
二、SP++3.0测试含matlab plot代码
注意:
此处必须先加载SP++3.0\include,然后再加载Matalb\extern\inlcude,因为SP++3.0 与Matlab中都有matrix.h头文件,如果加载顺序相反,则无法通过编译;
此处配置Matlab环境,参考文章:https://www.cnblogs.com/shuqingstudy/p/10134254.html
配置完成后,测试cppmatlab_test.cpp
/*****************************************************************************
* CppMatlab_test.cpp
*
* C++ and Matlab mixed programming testing.
*
* Zhang Ming, 2010-10, Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <cstring>
#include <vectormath.h>
#include <matrixmath.h>
#include <wft.h>
#include "engine.h"
using namespace std;
using namespace splab;
typedef double Type;
const int Lg = 128;
const int Ls = 1000;
const Type Fs = 1000;
int main()
{
/******************************* [ signal ] ******************************/
Vector<Type> t = linspace( Type(0), Type(Ls-1), Ls ) / Type(Fs);
Vector<Type> s = sin( Type(400*PI) * pow(t,Type(2.0)) );
/******************************** [ widow ] ******************************/
t = linspace(Type(0),Type(Lg-1),Lg);
Vector<Type> g = gauss( t, (Lg-1)/Type(2), Lg/Type(8) );
/********************************* [ WFT ] *******************************/
cout << "Taking windowed Fourier transform." << endl;
Matrix< complex<Type> > coefs = wft( s, g );
/******************************** [ IWFT ] *******************************/
cout << "Taking inverse windowed Fourier transform." << endl;
Vector<Type> x = iwft( coefs, g );
cout << "The relative error is : " << "norm(s-x) / norm(s) = "
<< norm(s-x)/norm(s) << endl << endl;
/******************************** [ PLOT ] *******************************/
Engine *ep = engOpen( NULL );
if( !ep )
{
cerr << "Cannot open Matlab Engine!" << endl;
exit(1);
}
Matrix<Type> C = trT( abs(coefs) );
int M = C.rows(), N = C.cols();
// define mxArray as 1-by-1 Real Scalar
mxArray *mFs = mxCreateDoubleMatrix( 1, 1, mxREAL );
// define mxArray as N-by-1 Real Vector
mxArray *ms = mxCreateDoubleMatrix( Ls, 1, mxREAL );
mxArray *mx = mxCreateDoubleMatrix( Ls, 1, mxREAL );
// define mxArray as N-by-M Real Matrix, BECAUSE matalb is ROW MAJOR
// and C/C++ is COLUMN MAJOR, the row of SP++ matrix is copied as the
// column of Matlab matrix
mxArray *mC = mxCreateDoubleMatrix( N, M, mxREAL );
// array copy from Scalar to mxArray.
memcpy( mxGetPr(mFs), &Fs, sizeof(Type) );
// array copy from Vectors to mxArray.
memcpy( mxGetPr(ms), s, Ls*sizeof(Type) );
memcpy( mxGetPr(mx), x, Ls*sizeof(Type) );
// array copy from Matrix to mxArray.
memcpy( mxGetPr(mC), C, C.size()*sizeof(Type) );
// send command to Matlab engine
engPutVariable( ep, "fs", mFs );
engPutVariable( ep, "s", ms );
engPutVariable( ep, "x", mx );
engPutVariable( ep, "C", mC );
// Matlab commands
const char *mCmd = " \
figure('name','C++ and Matlab Mixed Programming Testing'); \
hFN = floor(size(C,1)/2); tN = size(C,2); \
subplot(2,2,1); plot((0:tN-1), s); \
axis([0,tN,min(s),max(s)]); \
xlabel('Time (ms)', 'fontsize',12); ylabel('Amplitude', 'fontsize',12); \
title('(a)'); \
subplot(2,2,2); pcolor((0:tN-1),(0:hFN)'/hFN, C(1:hFN+1,:)); \
shading interp; \
axis([0,tN, 0,1]); \
yt = 0 : 0.2 : 1; set(gca, 'YTick',yt); set(gca, 'YTickLabel',fs*yt/2); \
xlabel('Time (ms)','fontsize',12); ylabel('Frequency (Hz)','fontsize',12); \
title('(b)'); \
subplot(2,2,3); plot((0:tN-1),x); \
axis([0,tN,min(x),max(x)]); \
xlabel('Time (ms)','fontsize',12); \
ylabel('Amplitude', 'fontsize',12); \
title('(c)'); \
subplot(2,2,4); e=s-x; plot((0:tN-1),e); \
axis([0,tN,min(e),max(e)]); \
xlabel('Time (ms)','fontsize',12); \
ylabel('Amplitude', 'fontsize',12); \
title('(d)'); \
";
// send command to Matlab engine
engEvalString( ep, mCmd );
// delete mxArray
mxDestroyArray( mFs );
mxDestroyArray( ms );
mxDestroyArray( mx );
mxDestroyArray( mC );
system( "pause" );
engClose(ep);
return 0;
}
输出:
项目代码缺陷:
目前VS2015环境仅仅配置了SP++库环境和Matlab环境,SP++内部使用fftw并没有配置,因此后面如果需要再进行配置。