窓付きFourier変換アルゴリズムのC++実装
9416 ワード
ヘッダファイル:
実装ファイル:
テストコード:
実行結果:
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), [email protected]
*
* 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
実装ファイル:
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), [email protected]
*
* 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 signal
Vector<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 function
for( int j=0; j<Lw; ++j )
sn[j] = tmp[i+j] * wn[j];
Sk = fft(sn);
// compute the Foureier transform
coefs.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 coefs
for( 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;
}
実行結果:
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-020
Process returned 0 (0x0) execution time : 0.156 s
Press any key to continue.