2019獨(dú)角獸企業(yè)重金招聘Python工程師標(biāo)準(zhǔn)>>>

頭文件：

/** 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*//****************************************************************************** wft.h** Windowed Fourier Transform.** These routines are designed for calculating WFT and IWFT of 1D signals.* The windows for forward and backword transform should be the same.** In order to eliminate border effect, the input signal is extended by* three forms: zeros padded("zpd"), periodized extension("ppd") and* symetric extension("sym").** Zhang Ming, 2010-03, Xi'an Jiaotong University.*****************************************************************************/#ifndef WFT_H #define WFT_H#include <string> #include <complex> #include <fft.h> #include <matrix.h> #include <utilities.h>namespace splab {template<typename Type>Matrix< complex<Type> > wft( const Vector<Type>&, const Vector<Type>&,const string &mode = "zpd" );template<typename Type>Vector<Type> iwft( const Matrix< complex<Type> >&, const Vector<Type>& );#include <wft-impl.h>} // namespace splab#endif // WFT_H

實(shí)現(xiàn)文件：

/** 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*//****************************************************************************** wft-impl.h** Implementation for Windowed Fourier Transform.** Zhang Ming, 2010-03, Xi'an Jiaotong University.*****************************************************************************//*** Compute WFT of 1D signal("xn") using "wn" as window. The time-frequency* coeffitions are stored in "coefs". The column represents time, and row* represents frequency.*/ template <typename Type> Matrix< complex<Type> > wft( const Vector<Type> &xn, const Vector<Type> &wn,const string &mode ) {int Lx = xn.size(),Lw = wn.size();// extends the input signalVector<Type> tmp = wextend( xn, Lw/2, "both", mode );Matrix< complex<Type> > coefs( Lw, Lx );Vector<Type> sn( Lw );Vector< complex<Type> > Sk( Lw );for( int i=0; i<Lx; ++i ){// intercept xn by wn functionfor( int j=0; j<Lw; ++j )sn[j] = tmp[i+j] * wn[j];Sk = fft(sn);// compute the Foureier transformcoefs.setColumn( Sk, i );}return coefs; }/*** Compute the inverse windowed Fourier transform of "coefs". The window* "wn" should be the same as forward transform. The reconstruction signal* is stored in "xn".*/ template <typename Type> Vector<Type> iwft( const Matrix< complex<Type> > &coefs,const Vector<Type> &wn ) {int Lw = wn.size(),Lx = coefs.cols();Vector<Type> xn(Lx);Matrix<Type> tmp( Lw, Lx );Vector< complex<Type> > Sk( Lw );Vector<Type> sn( Lw );// compute the inverse Fourier transform of coefsfor( int i=0; i<Lx; ++i ){Sk = coefs.getColumn(i);sn = ifftc2r(Sk);tmp.setColumn( sn, i );}int mid = Lw / 2;for( int i=0; i<Lx; ++i )xn[i] = tmp[mid][i] / wn[mid];return xn; }

測試代碼：

/****************************************************************************** wft_test.cpp** Windowed Fourier transform testing.** Zhang Ming, 2010-03, Xi'an Jiaotong University.*****************************************************************************/#define BOUNDS_CHECK#include <iostream> #include <vectormath.h> #include <timing.h> #include <wft.h>using namespace std; using namespace splab;typedef long double Type; const int Lg = 100; const int Ls = 1000; const int Fs = 1000;int main() {/******************************* [ signal ] ******************************/Type a = 0;Type b = Ls-1;Vector<Type> t = linspace(a,b,Ls) / Type(Fs);Vector<Type> s = sin( Type(400*PI) * pow(t,Type(2.0)) );/******************************** [ widow ] ******************************/a = 0;b = Type(Lg-1);Type u = (Lg-1)/Type(2);Type r = Lg/Type(8);t = linspace(a,b,Lg);Vector<Type> g = gauss(t,u,r);g = g/norm(g);/********************************* [ WFT ] *******************************/Type runtime = 0;Timing cnt;cout << "Taking windowed Fourier transform." << endl;cnt.start();Matrix< complex<Type> > coefs = wft( s, g );cnt.stop();runtime = cnt.read();cout << "The running time = " << runtime << " (ms)" << endl << endl;/******************************** [ IWFT ] *******************************/cout << "Taking inverse windowed Fourier transform." << endl;cnt.start();Vector<Type> x = iwft( coefs, g );cnt.stop();runtime = cnt.read();cout << "The running time = " << runtime << " (ms)" << endl << endl;cout << "The relative error is : " << "norm(s-x) / norm(s) = "<< norm(s-x)/norm(s) << endl << endl;return 0; }

運(yùn)行結(jié)果：

Taking windowed Fourier transform. The running time = 0.031 (ms)Taking inverse windowed Fourier transform. The running time = 0.047 (ms)The relative error is : norm(s-x) / norm(s) = 7.013e-020Process returned 0 (0x0) execution time : 0.156 s Press any key to continue.

轉(zhuǎn)載于:https://my.oschina.net/zmjerry/blog/4261

總結(jié)

以上是生活随笔為你收集整理的加窗Fourier变换算法的C++实现的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。