Files
radio/interpolation.c
T
jens 290fd5f5e0 Initial import
git-svn-id: http://moon:8086/svn/software/trunk/libsrc/radio@1 b431acfa-c32f-4a4a-93f1-934dc6c82436
2014-07-19 07:44:42 +00:00

831 lines
18 KiB
C
Executable File

// --------------------------------------------------------------
#include <string.h>
#include <malloc.h>
#include <math.h>
#include "interpolation.h"
// --------------------------------------------------------------
// Complex FIR filtering with downsampling
// --------------------------------------------------------------
void FirCpxInit(fir_cpx_t *pObj, uint32_t N, uint32_t M_up, uint32_t M_down)
{
pObj->pX = (cpx_t*)malloc(N*sizeof(cpx_t));
memset(pObj->pX, 0, N*sizeof(cpx_t));
pObj->N = N;
pObj->M_up = M_up;
pObj->M_down = M_down;
pObj->r = 0;
pObj->w = 0;
ClockInit(&pObj->clockUp, M_up);
ClockInit(&pObj->clockDown, M_down);
}
void FirCpxFree(fir_cpx_t *pObj)
{
free(pObj->pX);
}
uint32_t FirCpxProcessReal(fir_cpx_t *pObj, radio_float_t *pW, cpx_t *pSrc, cpx_t *pDst, uint32_t len)
{
uint32_t j, w, r, N;
uint32_t remain;
uint32_t nBytesWritten;
uint32_t nBytesRead;
const cpx_t zero = Cpx(0,0);
w = pObj->w;
r = pObj->r;
N = pObj->N;
nBytesWritten = 0;
nBytesRead = 0;
remain = len*pObj->M_up;
while(remain--)
{
if (ClockIsTick(&pObj->clockUp, 1))
{
pObj->pX[w] = pSrc[nBytesRead];
nBytesRead++;
}
else
{
pObj->pX[w] = zero;
}
if (++w >= N)
w = 0;
if (ClockIsTick(&pObj->clockDown, 1))
{
// fs/downSampleFactor
pDst[nBytesWritten] = zero;
for (j=0; j < N; j++)
{
pDst[nBytesWritten].real += pObj->pX[r].real * pW[j];
pDst[nBytesWritten].imag += pObj->pX[r].imag * pW[j];
if (--r >= N)
r = N-1;
}
nBytesWritten++;
}
if (++r >= N)
r = 0;
}
pObj->r = r;
pObj->w = w;
return nBytesWritten;
}
void FirCpxMultirateStageDownInit(fir_cpx_multirate_stage_t *pObj, radio_float_t *pW, uint32_t N, uint32_t M, uint32_t k)
{
uint32_t i;
radio_float_t *_pW;
pObj->N = 0;
for (i=k; i < N; i += M)
{
pObj->N++;
}
pObj->pW = (radio_float_t*)malloc(pObj->N*sizeof(cpx_t));
_pW = pObj->pW;
for (i=k; i < N; i += M)
{
*(_pW++) = pW[i];
}
pObj->k = k;
pObj->K = M;
pObj->pX = (cpx_t*)malloc(pObj->N*sizeof(cpx_t));
memset(pObj->pX, 0, pObj->N*sizeof(cpx_t));
}
void FirCpxMultirateStageUpInit(fir_cpx_multirate_stage_t *pObj, radio_float_t *pW, uint32_t N, uint32_t L, uint32_t k)
{
uint32_t i;
radio_float_t *_pW;
pObj->N = 0;
for (i=L-k-1; i < N; i += L)
{
pObj->N++;
}
pObj->pW = (radio_float_t*)malloc(pObj->N*sizeof(cpx_t));
_pW = pObj->pW;
if (pW)
{
for (i=L-k-1; i < N; i += L)
{
*(_pW++) = pW[i];
}
}
pObj->k = k;
pObj->K = L;
pObj->pX = (cpx_t*)malloc(pObj->N*sizeof(cpx_t));
memset(pObj->pX, 0, pObj->N*sizeof(cpx_t));
}
void FirCpxMultirateStageFree(fir_cpx_multirate_stage_t *pObj)
{
if (pObj->pX)
free(pObj->pX);
if (pObj->pW)
free(pObj->pW);
}
cpx_t FirCpxMultirateStageProcessDownAccum(fir_cpx_multirate_stage_t *pObj, cpx_t res, cpx_t src)
{
int32_t i;
for (i=pObj->N-1; i > 0; i--)
{
pObj->pX[i] = pObj->pX[i-1];
}
pObj->pX[i] = src;
return CpxMulRealAdd(res, pObj->pX, pObj->pW, pObj->N);
}
uint32_t FirCpxMultirateStageProcessUp(fir_cpx_multirate_stage_t *pObj, cpx_t *pSrc, cpx_t *pDst, uint32_t len)
{
int32_t i;
uint32_t num_written;
num_written = 0;
while(len--)
{
for (i=pObj->N-1; i > 0; i--)
{
pObj->pX[i] = pObj->pX[i-1];
}
pObj->pX[i] = *(pSrc++);
pDst[pObj->K-1 - pObj->k + num_written] = CpxMulRealAdd(Cpx(0,0), pObj->pX, pObj->pW, pObj->N);
num_written += pObj->K;
}
return num_written;
}
void FirCpxMultirateDownInit(fir_cpx_multirate_t *pObj, radio_float_t *pW, uint32_t N, uint32_t M)
{
uint32_t i;
pObj->L = 0;
pObj->M = M;
pObj->N = N;
pObj->k = M-1;
pObj->res = Cpx(0,0);
pObj->pStageUp = 0;
pObj->pStageDown = (fir_cpx_multirate_stage_t*)malloc(M*sizeof(fir_cpx_multirate_stage_t));
for (i=0; i < M; i++)
{
FirCpxMultirateStageDownInit(&pObj->pStageDown[i], pW, N, M, i);
}
}
void FirCpxMultirateUpInit(fir_cpx_multirate_t *pObj, radio_float_t *pW, uint32_t N, uint32_t L)
{
uint32_t i;
pObj->L = L;
pObj->M = 0;
pObj->N = N;
pObj->k = 0;
pObj->res = Cpx(0,0);
pObj->pStageDown = 0;
pObj->pStageUp = (fir_cpx_multirate_stage_t*)malloc(L*sizeof(fir_cpx_multirate_stage_t));
for (i=0; i < L; i++)
{
FirCpxMultirateStageUpInit(&pObj->pStageUp[i], pW, N, L, i);
}
}
void FirCpxMultirateFree(fir_cpx_multirate_t *pObj)
{
uint32_t i;
if (pObj->pStageDown)
{
for (i=0; i < pObj->M; i++)
{
FirCpxMultirateStageFree(&pObj->pStageDown[i]);
}
free(pObj->pStageDown);
}
if (pObj->pStageUp)
{
for (i=0; i < pObj->L; i++)
{
FirCpxMultirateStageFree(&pObj->pStageUp[i]);
}
free(pObj->pStageUp);
}
}
void FirCpxMultirateUpReinit(fir_cpx_multirate_t *pObj, radio_float_t *pW, uint32_t N, uint32_t L)
{
FirCpxMultirateFree(pObj);
FirCpxMultirateUpInit(pObj, pW, N, L);
}
uint32_t FirCpxDownZero(fir_cpx_multirate_t *pObj, cpx_t *pDst, uint32_t len)
{
uint32_t i;
uint32_t num_written;
num_written = 0;
for (i=0; i < len; i += pObj->M)
{
*(pDst++) = Cpx(0,0);
num_written++;
}
return num_written;
}
uint32_t FirCpxDownProcess(fir_cpx_multirate_t *pObj, cpx_t *pSrc, cpx_t *pDst, uint32_t len)
{
uint32_t num_written_soll;
uint32_t num_written_ist;
num_written_ist = 0;
num_written_soll = FirCpxDownZero(pObj, pDst, len);
while(len--)
{
pObj->res = FirCpxMultirateStageProcessDownAccum(&pObj->pStageDown[pObj->k], pObj->res, *(pSrc++));
pObj->k--;
if (pObj->k >= pObj->M)
{
pDst[num_written_ist++] = pObj->res;
pObj->res = Cpx(0,0);
pObj->k = pObj->M-1;
}
}
return num_written_ist;
}
uint32_t FirCpxUpZero(fir_cpx_multirate_t *pObj, cpx_t *pDst, uint32_t len)
{
uint32_t i;
uint32_t num_written;
num_written = 0;
for (i=0; i < pObj->L*len; i++)
{
*(pDst++) = Cpx(0,0);
num_written++;
}
return num_written;
}
uint32_t FirCpxUpProcess(fir_cpx_multirate_t *pObj, cpx_t *pSrc, cpx_t *pDst, uint32_t len)
{
uint32_t k;
uint32_t res;
uint32_t num_written_ist;
num_written_ist = 0;
for (k=0; k < pObj->L; k++)
{
res = FirCpxMultirateStageProcessUp(&pObj->pStageUp[k], pSrc, pDst, len);
if (num_written_ist < res)
num_written_ist = res;
}
return num_written_ist;
}
// --------------------------------------------------------------
// Complex FIR-based polyphase interpolation filter
// --------------------------------------------------------------
radio_float_t PolyPhaseCpxIpInit(ppip_cpx_t *pObj, radio_float_t fa, radio_float_t tsym, radio_float_t alpha, uint32_t N, uint32_t M)
{
uint32_t i, j;
fir_float_t *pCoeff;
radio_float_t k;
// Calculate polyphase filter
pCoeff = (fir_float_t*)malloc((N+1)*M*sizeof(fir_float_t));
memset(pCoeff, 0, (N+1)*M*sizeof(fir_float_t));
k = CalcFirSRRC(pCoeff, M*fa, tsym, alpha, N*M);
pObj->ppCoeff = (fir_float_t**)malloc(M*sizeof(fir_float_t*));
for(i=0; i < M; i++)
{
pObj->ppCoeff[i] = (fir_float_t*)malloc(N*sizeof(fir_float_t));
memset(pObj->ppCoeff[i], 0, N*sizeof(fir_float_t));
}
for(i=0; i < M; i++)
{
for(j=0; j < N; j++)
{
pObj->ppCoeff[i][j] = M*pCoeff[M*j+M/2+i];
}
}
// FIFO
pObj->L = N; // Influences STR-Loop delay
pObj->w = 0;
pObj->r = 0;
pObj->pFifo = (cpx_t*)malloc(pObj->L*sizeof(cpx_t));
memset(pObj->pFifo, 0, pObj->L*sizeof(cpx_t));
pObj->pX = (cpx_t*)malloc(pObj->L*sizeof(cpx_t));
memset(pObj->pX, 0, pObj->L*sizeof(cpx_t));
pObj->M = M;
pObj->N = N;
free(pCoeff);
return k/M;
}
void PolyPhaseCpxIpFree(ppip_cpx_t *pObj)
{
uint32_t i;
for(i=0; i < pObj->M; i++)
{
free(pObj->ppCoeff[i]);
}
free(pObj->ppCoeff);
pObj->ppCoeff = NULL;
free(pObj->pFifo);
pObj->pFifo = NULL;
free(pObj->pX);
pObj->pX = NULL;
}
void PolyPhaseCpxIpFeed(ppip_cpx_t *pObj, uint32_t push, cpx_t x)
{
pObj->w = (pObj->w + push) % pObj->L;
pObj->pFifo[pObj->w] = x;
}
cpx_t PolyPhaseCpxIpInterpolate(ppip_cpx_t *pObj, int32_t pop, radio_float_t mu)
{
uint32_t index, r;
int32_t j;
cpx_t y;
index = (uint32_t)dmod(dmax(0, mu*pObj->M), (radio_float_t)pObj->M);
pObj->r = (pObj->r + pop) % pObj->L;
if (pop)
{
r = pObj->r;
j = pObj->N-1;
while(r < pObj->L)
{
pObj->pX[j--] = pObj->pFifo[r++];
}
r = 0;
while(j >= 0)
{
pObj->pX[j--] = pObj->pFifo[r++];
}
}
y = CpxMulRealAdd(Cpx(0,0), pObj->pX, pObj->ppCoeff[index], pObj->N);
/*
pObj->r = (pObj->r + pop) % pObj->L;
r = pObj->r;
y = Cpx(0,0);
for (j=0; j < pObj->N; j++)
{
y.real += pObj->pFifo[r].real * pW[j];
y.imag += pObj->pFifo[r].imag * pW[j];
r = (r + pObj->L-1) % pObj->L;
}
*/
return y;
}
// --------------------------------------------------------------
// Complex FIR-based polyphase interpolation filter (Farrow structure)
// From:
// "PERFORMANCE AND DESIGN OF FARROW FILTER USED FOR ARBITRARY RESAMPLING"
// [unknown date], Fred Harris, Signal Processing ChairCommunication Systems and Signal Processing Institute
// College of Engineering, San Diego State University, San Diego, CA 92182-0190 USA
// --------------------------------------------------------------
void FarrowPPIPCpxInit(ppip_farrow_cpx_t *pObj, uint32_t M, uint32_t N)
{
uint32_t i;
pObj->ppCoeff = (fir_float_t**)malloc(M*sizeof(fir_float_t*));
for(i=0; i < M; i++)
{
pObj->ppCoeff[i] = (fir_float_t*)malloc(N*sizeof(fir_float_t));
memset(pObj->ppCoeff[i], 0, N*sizeof(fir_float_t));
}
pObj->M = M;
pObj->N = N;
pObj->pH = (cpx_t*)malloc(M*sizeof(cpx_t));
memset(pObj->pH, 0, M*sizeof(cpx_t));
pObj->pB = (cpx_t*)malloc(M*sizeof(cpx_t));
memset(pObj->pB, 0, M*sizeof(cpx_t));
pObj->L = N; // Influences STR-Loop delay
pObj->w = 0;
pObj->r = 0;
pObj->pFifo = (cpx_t*)malloc(pObj->L*sizeof(cpx_t));
memset(pObj->pFifo, 0, pObj->L*sizeof(cpx_t));
pObj->pX = (cpx_t*)malloc(pObj->L*sizeof(cpx_t));
memset(pObj->pX, 0, pObj->L*sizeof(cpx_t));
}
void FarrowPPIPCpxFree(ppip_farrow_cpx_t *pObj)
{
uint32_t i;
for(i=0; i < pObj->M; i++)
{
free(pObj->ppCoeff[i]);
}
free(pObj->ppCoeff);
pObj->ppCoeff = NULL;
free(pObj->pH);
pObj->pH = NULL;
free(pObj->pB);
pObj->pB = NULL;
free(pObj->pFifo);
pObj->pFifo = NULL;
free(pObj->pX);
pObj->pX = NULL;
}
cpx_t horner(cpx_t *a, cpx_t *b, uint32_t M, radio_float_t mu)
{
uint32_t i;
b[0] = a[0];
for (i=1; i < M; i++)
{
b[i].real = a[i].real+mu*b[i-1].real;
b[i].imag = a[i].imag+mu*b[i-1].imag;
}
return b[M-1];
}
void FarrowPPIPCpxFeed(ppip_farrow_cpx_t *pObj, uint32_t push, cpx_t x)
{
pObj->w = (pObj->w + push) % pObj->L;
pObj->pFifo[pObj->w] = x;
}
cpx_t FarrowPPIPCpxInterpolate(ppip_farrow_cpx_t *pObj, uint32_t pop, radio_float_t mu)
{
uint32_t i, r;
int32_t j;
pObj->r = (pObj->r + pop) % pObj->L;
if (pop)
{
r = pObj->r;
j = pObj->N-1;
while(r < pObj->L)
{
pObj->pX[j--] = pObj->pFifo[r++];
}
r = 0;
while(j >= 0)
{
pObj->pX[j--] = pObj->pFifo[r++];
}
}
// Partial filter responses
for (i=0; i < pObj->M; i++)
{
pObj->pH[i] = CpxMulRealAdd(Cpx(0,0), pObj->pX, pObj->ppCoeff[i], pObj->N);
}
// Combine
return horner(pObj->pH, pObj->pB, pObj->M, mu);
}
// --------------------------------------------------------------
// Real First-order interpolation filter
// --------------------------------------------------------------
void FDly_Init(fdly_t *pObj, int32_t order)
{
pObj->Nf = order + 1;
pObj->hf = (radio_float_t*)malloc(pObj->Nf*sizeof(radio_float_t));
pObj->statef = (radio_float_t*)malloc(pObj->Nf*sizeof(radio_float_t));
memset(pObj->hf, 0, pObj->Nf*sizeof(radio_float_t));
memset(pObj->statef, 0, pObj->Nf*sizeof(radio_float_t));
}
void FDly_Free(fdly_t *pObj)
{
if (pObj->hf)
free(pObj->hf);
pObj->hf = NULL;
if (pObj->statef)
free(pObj->statef);
pObj->statef = NULL;
}
radio_float_t FDly_Process(fdly_t *pObj, radio_float_t x)
{
radio_float_t y;
FIR(pObj->hf, pObj->statef, pObj->Nf, &x, &y, 1);
return y;
}
void FDly_CalcCoeff(fdly_t *pObj, radio_float_t delay)
{
int32_t n, k;
radio_float_t hn;
for (n=0; n < pObj->Nf; n++)
{
hn = 1.0;
for (k=0; k < pObj->Nf; k++)
{
if (k == n)
continue;
hn *= ((delay - (radio_float_t)k) / (radio_float_t)(n-k));
}
pObj->hf[n] = hn;
}
}
// --------------------------------------------------------------
// Real Polynomial-based Lagrange interpolation filter
// --------------------------------------------------------------
void LGIpInit(lgip_t *pObj, int32_t order, int32_t bufsize)
{
pObj->N = order + 1;
pObj->bufsize = bufsize;
RBufInit(&pObj->rbuf, bufsize);
}
void LGIpFree(lgip_t *pObj)
{
RBufFree(&pObj->rbuf);
pObj->N = 0;
}
radio_float_t LGIpProcess(lgip_t *pObj, radio_float_t x, int32_t m, radio_float_t mu)
{
int32_t n, k;
radio_float_t hn, y, delay;
RBufPut(&pObj->rbuf, x);
y = 0;
delay = 1-mu;
for (n=0; n < pObj->N; n++)
{
hn = 1.0;
for (k=0; k < pObj->N; k++)
{
if (k == n)
continue;
hn *= ((delay - (radio_float_t)k) / (radio_float_t)(n-k));
}
y += hn*RBufGetAt(&pObj->rbuf, pObj->bufsize/2 + m-n);
}
return y;
}
// --------------------------------------------------------------
// Complex Polynomial-based Lagrange interpolation filter
// --------------------------------------------------------------
void LGCpxIpInit(lgip_cpx_t *pObj, int32_t order)
{
// FIFO
pObj->N = order + 1;
pObj->w = 0;
pObj->r = 0;
pObj->L = 64*pObj->N; // Influences STR-Loop delay
pObj->pFifo = (cpx_t*)malloc(pObj->L*sizeof(cpx_t));
memset(pObj->pFifo, 0, pObj->L*sizeof(cpx_t));
}
void LGCpxIpFree(lgip_cpx_t *pObj)
{
free(pObj->pFifo);
pObj->pFifo = NULL;
pObj->N = 0;
}
cpx_t LGCpxIpProcess(lgip_cpx_t *pObj, cpx_t x, int32_t m, radio_float_t mu)
{
uint32_t i, j, k, N, L, r;
radio_float_t h;
cpx_t y;
L = pObj->L;
N = pObj->N;
pObj->pFifo[pObj->w++] = x;
if (pObj->w == L)
pObj->w = 0;
r = (m-N) % L;
k = 0;
y.real = 0.0;
y.imag = 0.0;
for (i=0; i < pObj->N; i++)
{
h = 1.0;
for (j=0; j < pObj->N; j++)
{
if (i == j)
continue;
h *= (((radio_float_t)pObj->N/2 + mu - (radio_float_t)j - (radio_float_t)0.5) / (radio_float_t)(i-j));
}
y.real += pObj->pFifo[r].real * h;
y.imag += pObj->pFifo[r].imag * h;
r = (r+1) % L;
}
return y;
}
// --------------------------------------------------------------
// Real FIR-based polyphase interpolation SRRC (Square Root Raised Cosine) matched filter #1
// --------------------------------------------------------------
radio_float_t PolyPhaseSRRCInit(fir_pp_t *pObj, radio_float_t fa, radio_float_t tsym, radio_float_t alpha, uint32_t N, uint32_t M)
{
uint32_t i, j;
fir_float_t *pCoeff;
radio_float_t k;
pCoeff = (fir_float_t*)malloc((N+1)*M*sizeof(fir_float_t));
memset(pCoeff, 0, M*(N+1)*sizeof(fir_float_t));
k = CalcFirSRRC(pCoeff, M*fa, tsym, alpha, M*N);
pObj->ppCoeff = (fir_float_t**)malloc(M*sizeof(fir_float_t*));
for(i=0; i < M; i++)
{
pObj->ppCoeff[i] = (fir_float_t*)malloc(N*sizeof(fir_float_t));
memset(pObj->ppCoeff[i], 0, N*sizeof(fir_float_t));
}
pObj->pState = (fir_float_t*)malloc(N*sizeof(fir_float_t));
memset(pObj->pState, 0, N*sizeof(fir_float_t));
for(i=0; i < M; i++)
{
for(j=0; j < N; j++)
{
pObj->ppCoeff[i][j] = M*pCoeff[M*j+M/2+i];
}
}
pObj->M = M;
pObj->N = N;
free(pCoeff);
return k/M;
}
void PolyPhaseSRRCFree(fir_pp_t *pObj)
{
uint32_t i;
for(i=0; i < pObj->M; i++)
{
free(pObj->ppCoeff[i]);
}
free(pObj->ppCoeff);
pObj->ppCoeff = NULL;
free(pObj->pState);
pObj->pState = NULL;
}
uint32_t PolyPhaseSRRCProcess(fir_pp_t *pObj, radio_float_t mu, radio_float_t *pX, radio_float_t *pY, uint32_t len)
{
uint32_t m;
// m = (pObj->M/2 + (int32_t)(mu * pObj->M/2)) % pObj->M;
m = (int32_t)dmod(mu*pObj->M, (radio_float_t)pObj->M);
FIR(pObj->ppCoeff[m], pObj->pState, pObj->N, pX, pY, len);
return len;
}
// --------------------------------------------------------------
// Real FIR-based polyphase interpolation SRRC (Square Root Raised Cosine) matched filter #2
// --------------------------------------------------------------
radio_float_t PolyPhaseIpInit(ppip_t *pObj, radio_float_t fa, radio_float_t tsym, radio_float_t alpha, uint32_t N, uint32_t M)
{
uint32_t i, j;
fir_float_t *pCoeff;
radio_float_t k;
pCoeff = (fir_float_t*)malloc((N+1)*M*sizeof(fir_float_t));
memset(pCoeff, 0, (N+1)*M*sizeof(fir_float_t));
k = CalcFirSRRC(pCoeff, M*fa, tsym, alpha, N*M);
pObj->ppCoeff = (fir_float_t**)malloc(M*sizeof(fir_float_t*));
for(i=0; i < M; i++)
{
pObj->ppCoeff[i] = (fir_float_t*)malloc(N*sizeof(fir_float_t));
memset(pObj->ppCoeff[i], 0, N*sizeof(fir_float_t));
}
for(i=0; i < M; i++)
{
for(j=0; j < N; j++)
{
pObj->ppCoeff[i][j] = M*pCoeff[M*j+M/2+i];
}
}
pObj->pState = (fir_float_t*)malloc(N*sizeof(fir_float_t));
memset(pObj->pState, 0, N*sizeof(fir_float_t));
pObj->M = M;
pObj->N = N;
RBufInit(&pObj->rbuf, 2*N); // Influences STR-Loop delay
free(pCoeff);
return k/M;
}
void PolyPhaseIpFree(ppip_t *pObj)
{
uint32_t i;
for(i=0; i < pObj->M; i++)
{
free(pObj->ppCoeff[i]);
}
free(pObj->ppCoeff);
pObj->ppCoeff = NULL;
free(pObj->pState);
pObj->pState = NULL;
RBufFree(&pObj->rbuf);
}
radio_float_t PolyPhaseIpProcess(ppip_t *pObj, radio_float_t x, int32_t m, radio_float_t mu)
{
uint32_t index, i;
radio_float_t y;
index = (int32_t)dmod(mu*pObj->M, (radio_float_t)pObj->M);
for(i=0; i < pObj->N; i++)
{
pObj->pState[i] = RBufGetAt(&pObj->rbuf, m-i);
}
FIR(pObj->ppCoeff[index], pObj->pState, pObj->N, &x, &y, 1);
RBufPut(&pObj->rbuf, x);
return y;
}