jens/avpflms
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
master
avpflms/pfft.c
T
jens e93e7927cb Initial import 
git-svn-id: http://moon:8086/svn/software/trunk/libsrc/avpflms@1 b431acfa-c32f-4a4a-93f1-934dc6c82436
2014-07-19 07:44:42 +00:00

357 lines
12 KiB
C
Executable File
Raw Permalink Blame History

/*-------------------------------------------------------------------------
* pfft.c Partioned Fast Fourier Transform
* $Id: $
*-------------------------------------------------------------------------
*
* Copyright (C) 2001 Algo Vision Systems GmbH
*
*-------------------------------------------------------------------------
*/
#define AVNEEDFLOAT
#ifdef _DEBUG
#include <stdio.h>
#endif
#include <math.h>
#include "avtypes.h"
#include "averror.h"
#include "avrtl.h"
#include "avfft.h"
#include "pfft.h"
/*************************************************************************/
AVERR PfftInit(
PFFT *pObj, /* Zeiger auf Objekt */
UINT32 N, /* Filterlaenge N */
UINT32 P, /* Anzahl Filterpartitionen */
UINT32 S, /* Anzahl Filtersegmente pro Partition */
UINT32 L, /* Laenge der Eingangsdaten (Blocklaenge */
UINT32 C /* FFT-Laenge */
)
{
UINT32 i, s, NP, C_L;
CMPXBUF *pX; /* Temporaerer Zeiger fuer Bufferinit.*/
/* Auto-Param */
if ((N*L)==0)
return AV_E_FAIL;
if (P==0)
{
if ((N%L) != 0)
return AV_E_FAIL;
P = N/L;
}
/* Auto-guess FFT-Groesse C */
if (C==0)
pObj->C = (UINT32)pow(2,ceil(log(L+N/P-1)/log(2.0)));
else
pObj->C = C;
/* Objekt initialisieren */
pObj->L = L;
pObj->N = N;
pObj->P = P;
pObj->S = S;
NP = pObj->S * pObj->L;
C_L = pObj->C - pObj->L;
/*-------------------------------------------------------*/
/* Speicher allokieren */
/*-------------------------------------------------------*/
/* Speicher fuer X[P*S][C] (komplex) */
pObj->pBufX = (CMPXBUF*)AvMemAlloc(pObj->P*pObj->S*sizeof(CMPXBUF));
for (i=0; i < pObj->P*pObj->S; i++)
{
pObj->pBufX[i].cmpxData.pReal = (avfloat_t*)AvMemAlloc(pObj->C*sizeof(avfloat_t));
pObj->pBufX[i].cmpxData.pImag = (avfloat_t*)AvMemAlloc(pObj->C*sizeof(avfloat_t));
AvZeroMem(pObj->pBufX[i].cmpxData.pReal,pObj->C*sizeof(avfloat_t));
AvZeroMem(pObj->pBufX[i].cmpxData.pImag,pObj->C*sizeof(avfloat_t));
pObj->pBufX[i].pNext = &pObj->pBufX[i+1];
pObj->pBufX[i].pLast = &pObj->pBufX[i-1];
pObj->pBufX[i].pLastPS = NULL;
pObj->pBufX[i].user = i;
}
pObj->pBufX[i-1].pNext = &pObj->pBufX[0];
pObj->pBufX[0].pLast = &pObj->pBufX[i-1];
/* Zeiger auf X[k-p*S] */
for (i=0; i < pObj->P*pObj->S; i++)
{
pX = pObj->pBufX[i].pLast;
for (s=1; s < pObj->S; s++)
pX = pX->pLast;
pObj->pBufX[i].pLastPS = pX;
}
/* Ergebnis 'Y' der Faltung (komplex) */
pObj->pY = (COMPLEX*)AvMemAlloc(sizeof(COMPLEX));
pObj->pY->pReal = (avfloat_t*)AvMemAlloc(pObj->C*sizeof(avfloat_t));
pObj->pY->pImag = (avfloat_t*)AvMemAlloc(pObj->C*sizeof(avfloat_t));
AvZeroMem(pObj->pY->pReal,pObj->C*sizeof(avfloat_t));
AvZeroMem(pObj->pY->pImag,pObj->C*sizeof(avfloat_t));
/* Overlap-Save 'xs' (reell) */
pObj->BufXsave.pReal = (avfloat_t*)AvMemAlloc(C_L*sizeof(avfloat_t));
pObj->BufXsave.pImag = (avfloat_t*)AvMemAlloc(C_L*sizeof(avfloat_t));
AvZeroMem(pObj->BufXsave.pReal, C_L*sizeof(avfloat_t));
AvZeroMem(pObj->BufXsave.pImag, C_L*sizeof(avfloat_t));
/* FFT initialisieren */
pObj->pFFT = (FFT*)AvMemAlloc(sizeof(FFT));
FFTinit(pObj->pFFT, pObj->C);
/* Arbeitszeiger initialisieren */
pObj->pX = pObj->pBufX;
return AV_E_OK;
}
/*-----------------------------------------------------------------------*/
/* Partitioned FLMS Filterinitialisierung
/* 1. Partitionierung der Filterstartwerte in P-Partitionen
/* 2. Transformation der P Teil-Filter in den Frequenzbereich
/* Element pBufWS wird veraendert
/*-----------------------------------------------------------------------*/
AVERR PfftFilterInit(
PFFT *pObj, /* Zeiger auf Objekt */
COMPLEX WTD,
CMPXBUF *pWS) /* N Filterkoeffizienten im Zeitbereich */
{
UINT32 p, NP, C_SL;
NP = pObj->S * pObj->L;
C_SL = pObj->C - NP;
for (p=0; p < pObj->P; p++)
{
/* Arbeitspuffer WS auffuellen
Re{WS[p][0..C-L-1]}=wi[p*N/P..(p+1)*N/P-1] */
AvMemCpy(pWS[p].cmpxData.pReal,
&WTD.pReal[p*NP], NP*sizeof(avfloat_t));
/* Arbeitspuffer WS auffuellen
Im{WS[p][0..C-L-1]}=wi[p*N/P..(p+1)*N/P-1] */
AvMemCpy(pWS[p].cmpxData.pImag,
&WTD.pImag[p*NP], NP*sizeof(avfloat_t));
/* Re{WS[C-L..C]} = {0} (vorsichtshalber) */
AvZeroMem(&pWS[p].cmpxData.pReal[NP], C_SL *sizeof(avfloat_t));
/* Im{WS[C-L..C]} = {0} (vorsichtshalber) */
AvZeroMem(&pWS[p].cmpxData.pImag[NP], C_SL *sizeof(avfloat_t));
/* In den Frequenzbereich transformieren
WS[p][0..C-1] = 1/C*FFT{wi[p*N/P..(p+1)*N/P-1]} */
ffts(pObj->pFFT, pWS[p].cmpxData.pReal, pWS[p].cmpxData.pImag);
}
return AV_E_OK;
}
avfloat_t *PfftGetBufInRe(PFFT *pObj)
{
UINT32 C_L;
C_L = pObj->C - pObj->L; /* C-L */
return &pObj->pX->cmpxData.pReal[C_L];
}
avfloat_t *PfftGetBufInIm(PFFT *pObj)
{
UINT32 C_L;
C_L = pObj->C - pObj->L; /* C-L */
return &pObj->pX->cmpxData.pImag[C_L];
}
avfloat_t *PfftGetBufOutRe(PFFT *pObj)
{
UINT32 C_L;
C_L = pObj->C - pObj->L; /* C-L */
return &pObj->pY->pReal[C_L];
}
avfloat_t *PfftGetBufOutIm(PFFT *pObj)
{
UINT32 C_L;
C_L = pObj->C - pObj->L; /* C-L */
return &pObj->pY->pImag[C_L];
}
/*---------------------------------------------------------------*/
/* Speicher zuweisen fuer Koeffizienten H
/*---------------------------------------------------------------*/
AVERR PfftFilterAlloc(PFFT *pObj, CMPXBUF **ppH)
{
UINT32 i;
/* Speicher fuer WS[P][C] (komplex) */
*ppH = (CMPXBUF*)AvMemAlloc(pObj->P*pObj->S*sizeof(CMPXBUF));
for (i=0; i < pObj->P*pObj->S; i++)
{
(*ppH)[i].cmpxData.pReal = (avfloat_t*)AvMemAlloc(pObj->C*sizeof(avfloat_t));
(*ppH)[i].cmpxData.pImag = (avfloat_t*)AvMemAlloc(pObj->C*sizeof(avfloat_t));
AvZeroMem((*ppH)[i].cmpxData.pReal,pObj->C*sizeof(avfloat_t));
AvZeroMem((*ppH)[i].cmpxData.pImag,pObj->C*sizeof(avfloat_t));
(*ppH)[i].pNext = &(*ppH)[i+1];
(*ppH)[i].pLast = &(*ppH)[i-1];
(*ppH)[i].pLastPS = NULL;
(*ppH)[i].user = i;
}
(*ppH)[i-1].pNext = &(*ppH)[0];
(*ppH)[0].pLast = &(*ppH)[i-1];
return AV_E_OK;
}
/*---------------------------------------------------------------*/
/* Partitioned FFT-Filterung
/*---------------------------------------------------------------*/
AVERR PfftFilter(PFFT *pObj, CMPXBUF *pH, COMPLEX x, COMPLEX y)
{
UINT32 p, NP, C_L;
CMPXBUF *pX;
COMPLEX *pY;
/* Arbeitszeiger */
pX = pObj->pX; /* Aktueller Zeiger X[k] */
pY = pObj->pY; /* 'Y' wird nach der Filterung als
temporaerer Speicher benutzt */
NP = pObj->S*pObj->L; /* N/P = S*L */
C_L = pObj->C - pObj->L; /* C-L */
/* Arbeitspuffer 'X' auffuellen Re{X[S*L..C-1]} = x[0..L-1] */
AvMemCpy(&pX->cmpxData.pReal[C_L],x.pReal,
pObj->L*sizeof(avfloat_t));
/* Saveblock 'xs' anfuegen Re{X[0..C-L-1]} = xs[0..C-L-1] */
AvMemCpy(pX->cmpxData.pReal,pObj->BufXsave.pReal,
C_L*sizeof(avfloat_t));
/* Saveblock aktualisieren xs[0..C-L-1] = x[L..C-1] */
AvMemCpy(pObj->BufXsave.pReal, &pX->cmpxData.pReal[pObj->L],
C_L*sizeof(avfloat_t));
/* Arbeitspuffer 'X' auffuellen Im{X[S*L..C-1]} = x[0..L-1] */
AvMemCpy(&pX->cmpxData.pImag[C_L],x.pImag,
pObj->L*sizeof(avfloat_t));
/* Saveblock 'xs' anfuegen Im{X[0..C-L-1]} = xs[0..C-L-1] */
AvMemCpy(pX->cmpxData.pImag,pObj->BufXsave.pImag,
C_L*sizeof(avfloat_t));
/* Saveblock aktualisieren xs[0..C-L-1] = x[L..C-1] */
AvMemCpy(pObj->BufXsave.pImag, &pX->cmpxData.pImag[pObj->L],
C_L*sizeof(avfloat_t));
/* X = FFT{x} */
fft(pObj->pFFT, pX->cmpxData.pReal, pX->cmpxData.pImag);
/* 1. Partition Multiplikation im Frequenzbereich
Y = X[k][0..C-1] * H[0][0..C-1] */
CmpxVectMul(&pX->cmpxData, &pH->cmpxData,
pY, pObj->C);
/* 2. Partition bis P-te Partition */
for (p=1; p < pObj->P; p++)
{
pH = pH->pNext;
/* X[k-p*S] suchen */
pX = pX->pLastPS;
/* Y = X[k-p*S][0..C-1] * H[p][0..C-1] */
CmpxVectMac(&pX->cmpxData, &pH->cmpxData,
pY, pObj->C);
}
/* In den Zeitbereich transformieren, y = IFFT{Y} */
ifft(pObj->pFFT, pY->pReal, pY->pImag);
/* Abspeichern der letzten L Daten ys[0..L-1] = Re{Y[C-L..C-1] */
AvMemCpy(y.pReal, &pY->pReal[C_L], pObj->L*sizeof(avfloat_t));
/* Abspeichern der letzten L Daten ys[0..L-1] = Im{Y[C-L..C-1] */
AvMemCpy(y.pImag, &pY->pImag[C_L], pObj->L*sizeof(avfloat_t));
/* Fuer den naechsten Aufruf Zeiger aktualisieren */
/* Naechstes X[k] ist: */
pObj->pX = pObj->pX->pNext;
return AV_E_OK;
}
AVERR PfftFilterFast(PFFT *pObj, CMPXBUF *pH)
{
UINT32 p, NP, C_L;
CMPXBUF *pX;
COMPLEX *pY;
/* Arbeitszeiger */
pX = pObj->pX; /* Aktueller Zeiger X[k] */
pY = pObj->pY; /* 'Y' wird nach der Filterung als
temporaerer Speicher benutzt */
NP = pObj->S*pObj->L; /* N/P = S*L */
C_L = pObj->C - pObj->L; /* C-L */
/* Saveblock 'xs' anfuegen Re{X[0..C-L-1]} = xs[0..C-L-1] */
AvMemCpy(pX->cmpxData.pReal,pObj->BufXsave.pReal,
C_L*sizeof(avfloat_t));
/* Saveblock aktualisieren xs[0..C-L-1] = x[L..C-1] */
AvMemCpy(pObj->BufXsave.pReal, &pX->cmpxData.pReal[pObj->L],
C_L*sizeof(avfloat_t));
/* Saveblock 'xs' anfuegen Im{X[0..C-L-1]} = xs[0..C-L-1] */
AvMemCpy(pX->cmpxData.pImag,pObj->BufXsave.pImag,
C_L*sizeof(avfloat_t));
/* Saveblock aktualisieren xs[0..C-L-1] = x[L..C-1] */
AvMemCpy(pObj->BufXsave.pImag, &pX->cmpxData.pImag[pObj->L],
C_L*sizeof(avfloat_t));
/* X = FFT{x} */
fft(pObj->pFFT, pX->cmpxData.pReal, pX->cmpxData.pImag);
/* 1. Partition Multiplikation im Frequenzbereich
Y = X[k][0..C-1] * H[0][0..C-1] */
CmpxVectMul(&pX->cmpxData, &pH->cmpxData,
pY, pObj->C);
/* 2. Partition bis P-te Partition */
for (p=1; p < pObj->P; p++)
{
pH = pH->pNext;
/* X[k-p*S] suchen */
pX = pX->pLastPS;
/* Y = X[k-p*S][0..C-1] * H[p][0..C-1] */
CmpxVectMac(&pX->cmpxData, &pH->cmpxData,
pY, pObj->C);
}
/* In den Zeitbereich transformieren, y = IFFT{Y} */
ifft(pObj->pFFT, pY->pReal, pY->pImag);
/* Fuer den naechsten Aufruf Zeiger aktualisieren */
/* Naechstes X[k] ist: */
pObj->pX = pObj->pX->pNext;
return AV_E_OK;
}
/*------------------------------------------------------------------*/
Reference in New Issue View Git Blame Copy Permalink