// -------------------------------------------------------------- #include #include #include #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; }