1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
|
#include <iostream>
#include <complex>
#include <cmath>
#include <iterator>
using namespace std;
unsigned int bitReverse(unsigned int x, int log2n)
{
int n = 0;
for (int i=0; i < log2n; i++)
{
n <<= 1;
n |= (x & 1);
x >>= 1;
}
return n;
}
const double PI = 3.1415926536;
template<class Iter_T>
void fft(Iter_T a, Iter_T b, int log2n)
{
const complex<double> J(0, 1);
int n = 1 << log2n;
for (unsigned int i=0; i < n; ++i)
{
b[bitReverse(i, log2n)] = a[i];
}
for (int s = 1; s <= log2n; ++s)
{
int m = 1 << s;
int m2 = m >> 1;
complex<double>w(1, 0);
complex<double> wm = exp(-J * (PI / m2));
for (int j=0; j < m2; ++j)
{
for (int k=j; k < n; k += m)
{
complex <double> t = w * b[k + m2];
complex <double> u = b[k];
b[k] = u + t;
b[k + m2] = u - t;
}
w *= wm;
}
}
}
template <typename Type>
Vector<Type> yulewalkerPSE( const Vector<Type> &xn, int p, Type &sigma2 )
{
int N = xn.size();
assert( p <= N );
Vector<Type> rn(p+1);
for( int i=0; i<=p; ++i )
for( int k=0; k<N-i; ++k )
rn[i] += xn[k+i]*xn[k];
rn /= Type(N);
return levinson( rn, sigma2 );
}
template <typename Type>
Vector<Type> levinson( const Vector<Type> &rn, Type &sigma2 )
{
int p = rn.size()-1;
Type tmp;
Vector<Type> ak(p+1), akPrev(p+1);
ak[0] = Type(1.0);
sigma2 = rn[0];
ak[1] = -rn[1]/sigma2;
sigma2 *= 1 - ak[1]*ak[1];
for( int k=2; k<=p; ++k )
{
tmp = 0;
for( int i=0; i<k; ++i )
tmp += ak[i]*rn[k-i];
ak[k] = -tmp/sigma2;
for( int i=1; i<k; ++i )
akPrev[i] = ak[i] + ak[k]*ak[k-i];
for( int i=1; i<k; ++i )
ak[i] = akPrev[i];
sigma2 *= 1 - ak[k]*ak[k];
}
return ak;
}
typedef double Type;
int main( )
{
typedef complex <double> cx;
Type sigma2;
cx a[127];
for(int i=0;i<127;i++)
{
a[i]=sin(2*PI*10*i*0.005);
}
cx b[127];
fft(a, b, 7);
for (int i=0; i<128; ++i)
{
cout << b[i] << "\n";
}
int y=yulewalkerPSE(b,4,sigma2);
cout<<y<<"\n";
getchar();
return 0;
}
|