📄 rec_signal.cpp
字号:
int count, i, checksum, o1, o2;
CSignal Aux(1023);
switch (PRNNo)
{
case 1:
tap1 = 2;
tap2 = 6;
break;
case 2:
tap1 = 3;
tap2 = 7;
break;
case 3:
tap1 = 4;
tap2 = 8;
break;
case 4:
tap1 = 5;
tap2 = 9;
break;
case 5:
tap1 = 1;
tap2 = 9;
break;
case 6:
tap1 = 2;
tap2 = 10;
break;
case 7:
tap1 = 1;
tap2 = 8;
break;
case 8:
tap1 = 2;
tap2 = 9;
break;
case 9:
tap1 = 3;
tap2 = 10;
break;
case 10:
tap1 = 2;
tap2 = 3;
break;
case 11:
tap1 = 3;
tap2 = 4;
break;
case 12:
tap1 = 5;
tap2 = 6;
break;
case 13:
tap1 = 6;
tap2 = 7;
break;
case 14:
tap1 = 7;
tap2 = 8;
break;
case 15:
tap1 = 8;
tap2 = 9;
break;
case 16:
tap1 = 9;
tap2 = 10;
break;
case 17:
tap1 = 1;
tap2 = 4;
break;
case 18:
tap1 = 2;
tap2 = 5;
break;
case 19:
tap1 = 3;
tap2 = 6;
break;
case 20:
tap1 = 4;
tap2 = 7;
break;
case 21:
tap1 = 5;
tap2 = 8;
break;
case 22:
tap1 = 6;
tap2 = 9;
break;
case 23:
tap1 = 1;
tap2 = 3;
break;
case 24:
tap1 = 4;
tap2 = 6;
break;
case 25:
tap1 = 5;
tap2 = 7;
break;
case 26:
tap1 = 6;
tap2 = 8;
break;
case 27:
tap1 = 7;
tap2 = 9;
break;
case 28:
tap1 = 8;
tap2 = 10;
break;
case 29:
tap1 = 1;
tap2 = 6;
break;
case 30:
tap1 = 2;
tap2 = 7;
break;
case 31:
tap1 = 3;
tap2 = 8;
break;
case 32:
tap1 = 4;
tap2 = 9;
break;
default:
cerr << " wrong prn No " << endl;
exit(1);
}
// initialize variables
count = 0;
for (i = 0; i < 10; i++ )
{
reg1[i] = -1;
reg2[i] = -1;
check[i] = 1;
}
checksum = 10;
while ( checksum != -10 )
{
o1 = reg1[2] * reg1[9];
o2 = reg2[1] * reg2[2] * reg2[5] * reg2[7] * reg2[8] * reg2[9];
Aux.data[count] = -reg1[9] * reg2[tap1-1] * reg2[tap2-1];
count++;
for (i = 8; i>=0; i--)
{
reg1[i+1] = reg1[i];
reg2[i+1] = reg2[i];
}
reg1[0] = o1;
reg2[0] = o2;
checksum = 0;
for ( i=0; i < 10; i++)
{
checksum += reg2[i];
}
} // end while (checksum != -10 )
return(Aux);
}
//----------------------------------------------------------------------------------------//
// class CFixedFFT implementation //
// //
// Purpose: //
// //
// to implement fixed fft //
//----------------------------------------------------------------------------------------//
// Constructor, destructor
CFixedFFT::CFixedFFT(unsigned int _n):m_n(_n)
{
if(_n<0 || (_n&(_n-1))!=0 ){
cerr<<"Error: illegal parameter in CFixedFFT, exit(0)"<<endl;
exit(1);
}
unsigned int i,j,l, L;
unsigned int n;
unsigned int LE; /* Number of points in sub DFT at stage L
and offset to next DFT in stage */
unsigned int LE1; /* Number of butterflies in one DFT at
stage L. Also is offset to lower point
in butterfly at stage L */
EXP = 0; // Log2(m_n)
while((1<<EXP++) != m_n);
EXP--;
n = m_n - (m_n>>2);
Sinewave = new int [n];
for (L=0; L<n; L++) /* Create twiddle factor table */
{
Sinewave[L] = (int)((0x7fff*sin(pi2*L/m_n))+0.5);
}
}
CFixedFFT::~CFixedFFT()
{
delete [] Sinewave;
}
inline int FIX_MPY(int a, int b)
{
/* shift right one less bit (i.e. 15-1) */
int c = ((int)a * (int)b) >> 14;
/* last bit shifted out = rounding-bit */
b = c & 0x01;
/* last shift + rounding bit */
a = (c >> 1) + b;
return a;
}
//----------------------------------------------------------------
// FFTForward
// realize FFT transform
//
// IN:
// _re ------------ real part of input signal
// _im ------------ image part of input signal
// inverse == 0, forward fft, otherwise backward fft
// Out:
// _re ------------ real part of output signal
// _im ------------ image part of output signal
//
int CFixedFFT::FFT(CSignal& _re, CSignal& _im, short inverse)
{
accuracy *fr = _re.data;
accuracy *fi = _im.data;
int m;
int N_WAVE = m_n;
int mr, nn, i, j, l, k, istep, scale, shift;
int qr, qi, tr, ti, wr, wi;
int n = _re.m_n;
/* max FFT size = N_WAVE */
if (n > N_WAVE)
return -1;
mr = 0;
nn = n - 1;
scale = 0;
/* decimation in time - re-order data */
for (m=1; m<=nn; ++m) {
l = n;
do {
l >>= 1;
} while (mr+l > nn);
mr = (mr & (l-1)) + l;
if (mr <= m)
continue;
tr = fr[m];
fr[m] = fr[mr];
fr[mr] = tr;
ti = fi[m];
fi[m] = fi[mr];
fi[mr] = ti;
}
l = 1;
k = EXP-1;
while (l < n) {
if (inverse) {
/* variable scaling, depending upon data */
shift = 0;
for (i=0; i<n; ++i) {
j = fr[i];
if (j < 0)
j = -j;
m = fi[i];
if (m < 0)
m = -m;
if (j > 16383 || m > 16383) {
shift = 1;
break;
}
}
if (shift)
++scale;
} else {
/*
fixed scaling, for proper normalization --
there will be log2(n) passes, so this results
in an overall factor of 1/n, distributed to
maximize arithmetic accuracy.
*/
shift = 1;
}
/*
it may not be obvious, but the shift will be
performed on each data point exactly once,
during this pass.
*/
istep = l << 1;
for (m=0; m<l; ++m) {
j = m << k;
/* 0 <= j < N_WAVE/2 */
wr = Sinewave[j+N_WAVE/4];
wi = -Sinewave[j];
if (inverse)
wi = -wi;
if (shift) {
wr >>= 1;
wi >>= 1;
}
for (i=m; i<n; i+=istep) {
j = i + l;
tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
qr = fr[i];
qi = fi[i];
if (shift) {
qr >>= 1;
qi >>= 1;
}
fr[j] = qr - tr;
fi[j] = qi - ti;
fr[i] = qr + tr;
fi[i] = qi + ti;
}
}
--k;
l = istep;
}
for( i = 0; i<n; i++ )
{
fr[i] = fr[i]<<scale;
fi[i] = fi[i]<<scale;
}
return scale;
}
CSignal sig_sqrt(const CSignal& _sig)
{
CSignal Aux(_sig.m_n);
int i;
for( i=0; i<_sig.m_n; i++ )
{
Aux.data[i] = fix_sqrt(_sig.data[i]);
}
return Aux;
}
int fix_sqrt (int x)
{
if( x<0)
{
cerr<<"Error: argument less than zero in fix_sqrt(), exit"<<endl;
exit(1);
}
int xt, scr;
int i;
i = 0;
xt = x;
do {
xt = xt >> 1;
i++;
}while (xt > 0);
i = (i >> 1) + 1;
xt = x >> i;
do {
scr = xt * xt;
scr = x - scr;
scr = scr >> 1;
scr = scr / xt;
xt = scr + xt;
}while (scr != 0);
if(xt*xt-x >= xt && xt != 0 ) xt--;
//xt = xt << 7;
return (xt);
}
/******************************************************************************
FUNCTION rss(long a, long b)
RETURNS long integer
PARAMETERS
a long integer
b long integer
PURPOSE
This function finds the fixed point magnitude of a 2 dimensional vector
WRITTEN BY
Clifford Kelley
******************************************************************************/
int rss (int a, int b)
{
int result, c, d;
c = labs (a);
d = labs (b);
if (c == 0 && d == 0)
result = 0;
else {
if (c > d)
result = (d >> 1) + c;
else
result = (c >> 1) + d;
}
return (result);
}
/******************************************************************************
FUNCTION fix_atan2(long y,long x)
RETURNS long integer
PARAMETERS
x long in-phase fixed point value
y long quadrature fixed point value
PURPOSE
This function computes the fixed point arctangent represented by
x and y in the parameter list
1 radian = 16384
based on the power series f-f^3*2/9
WRITTEN BY
Clifford Kelley
Fixed for y==x added special code for x==0 suggested by Joel
Barnes, UNSW
******************************************************************************/
int fix_atan2 (int y, int x)
{
static int const SCALED_PI_ON_2 = 25736L;
static int const SCALED_PI = 51472L;
long result = 0, n, n3;
if ((x == 0) && (y == 0))
return (0); /* invalid case */
if (x > 0 && x >= labs (y)) {
n = (y << 14) / x;
n3 = ((((n * n) >> 14) * n) >> 13) / 9;
result = n - n3;
}
else if (x <= 0 && -x >= labs (y)) {
n = (y << 14) / x;
n3 = ((((n * n) >> 14) * n) >> 13) / 9;
if (y > 0)
result = n - n3 + SCALED_PI;
else if (y <= 0)
result = n - n3 - SCALED_PI;
}
else if (y > 0 && y > labs (x)) {
n = (x << 14) / y;
n3 = ((((n * n) >> 14) * n) >> 13) / 9;
result = SCALED_PI_ON_2 - n + n3;
}
else if (y < 0 && -y > labs (x)) {
n = (x << 14) / y;
n3 = ((((n * n) >> 14) * n) >> 13) / 9;
result = -n + n3 - SCALED_PI_ON_2;
}
return (result);
}
///////////////////////////////////////////////////////////////////////////
// fft(double reel[], double imag[], int n, int sign) ///
// input: ///
// *reel -- real part of the signal ///
// *imag -- imag part of the signal ///
// mm -- length of the signal, max be the pow of two ///
// sign-- sign of the fft, sign= 1 for forward transform ///
// sign= -1 for backward transform ///
// output: ///
// *reel -- real part of the fft ///
// *imag -- image part of the fft ///
// Author: YI Weiyong
//template <class type1>
void fft(float *reel, float *imag, int mm, int sign)
{
int n, m, m2, i, j, k, l, log2n;
float c1, c2, norm, norm2, cphi, sphi, temp;
if ( mm &(mm - 1) )
{
cerr << "Error: the length of the data in fft is not the power of 2, exit(1) " << endl;
exit(1);
}
m = 1;
for (i = 1; i < mm; i++ )
{
m <<= 1;
if ( m == mm )
{
log2n = i;
break;
}
}
n = 1<<log2n;
/* Inversement des bits */
for(i=0; i<n; i++)
{
for(j=log2n-1, m=0, k=i; j>=0; j--, k>>=1)
m += (k&1)<<j;
if(m > i)
{
temp = reel[i];
reel[i] = reel[m];
reel[m] = temp;
temp = imag[i];
imag[i] = imag[m];
imag[m] = temp;
}
}
/* normalisation */
/* norm = 1.0/sqrt((double)n);
for(i=0; i<n ;i++) {
reel[i] *= norm;
imag[i] *= norm;
}*/
/* calculate fft */
for( j = 0; j < log2n; j++ ) {
m = 1<<j;
m2 = 2*m;
c1 = 1.0;
c2 = 0.0;
cphi = (float)cos(sign*2.0*pi/((double)m2));
sphi = (float)sin(- sign*2.0*pi/((double)m2));
for(k=0; k<m; k++)
{
for(i=k; i<n; i+=m2)
{
l = i + m;
norm = c1*reel[l] - c2*imag[l];
norm2 = c1*imag[l] + c2*reel[l];
reel[l] = reel[i] - norm;
imag[l] = imag[i] - norm2;
reel[i] += norm;
imag[i] += norm2;
}
norm = c1*cphi - c2*sphi;
norm2 = c1*sphi + c2*cphi;
c1 = norm; c2 = norm2;
}
}
if ( sign == -1 )
{
for (i = 0; i < n; i++ )
{
reel[i] = reel[i]/n;
imag[i] = imag[i]/n;
}
}
}
unsigned int bit_reverse(unsigned int x, unsigned int m)
{
if ( x >= (1<<m) )
{
cerr << " Warning: There might be some error in bit_reverse fuction, x > 2^m" <<endl;
}
int n = 0;
for (int i=0; i < m; i++)
{
n <<= 1;
n |= (x & 1);
x >>= 1;
}
return n;
}
///////////////////////////////////////////////////////////////////////////
// fft(CSignal& _real, CSignal& _imag, int sign) ///
// input: ///
// _real -- real part of the signal ///
// _imag -- imag part of the signal ///
// sign-- sign of the fft, sign= 1 for forward transform ///
// sign= -1 for backward transform ///
// output: ///
// _real -- real part of the fft ///
// _imag -- image part of the fft ///
// Author: YI Weiyong
void fft(CSignal& _real, CSignal& _imag, int sign)
{
float *real, *imag;
int n = (int)_real.m_n;
int i,j;
real = new float [n];
imag = new float [n];
for(i=0; i<n; i++)
{
real[i] = (float) _real.data[i];
imag[i] = (float) _imag.data[i];
}
fft(real, imag, n, sign);
for(i=0; i<n; i++)
{
_real.data[i] = (int)real[i];
_imag.data[i] = (int)imag[i];
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -