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

627 lines
17 KiB
C
Executable File
Raw Blame History

/***************************************************************************/
/* FFT.C
/* Fast-Fourier-Transformation
/* Author: Jens Ahrensfeld
/* Datum : 24.06.1999
/* letzte Änderung: 09.06.2000
/***************************************************************************/
#define AVNEEDFLOAT
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.1415926535897932384626433832795
/***************************************************************************/
/* Konstruktor FFT-Objekt */
/* */
/***************************************************************************/
void FFTinit(fft_t *pFFT, UINT32 N)
{
pFFT->m_numPoints = 0;
pFFT->m_numStages = 0;
pFFT->pTwfRe = NULL;
pFFT->pTwfIm = NULL;
if(IsPowerOfTwo(N) == FFT_ERROR)
{
printf("\n%u is not a power of 2!\n",N);
exit(1);
}
pFFT->m_numStages = (UINT32) (1e-06 + log((fft_float_t)N)/log(2.0));
pFFT->m_numPoints = (UINT32 )N;
pFFT->pTwfRe = (fft_float_t*)malloc(pFFT->m_numPoints * sizeof(fft_float_t) /2);
pFFT->pTwfIm = (fft_float_t*)malloc(pFFT->m_numPoints * sizeof(fft_float_t) /2);
FFTCalcTwiddleTable(pFFT->pTwfRe, pFFT->pTwfIm, pFFT->m_numPoints);
#ifdef FFTMsg
printf("\n%d-Point FFT object created.\n",pFFT->m_numPoints);
#endif
}
/***************************************************************************/
/* Destruktor FFT-Objekt */
/* */
/***************************************************************************/
void FFTfree(fft_t *pFFT)
{
if (pFFT->pTwfRe != NULL)
{
free(pFFT->pTwfRe);
pFFT->pTwfRe = NULL;
}
if (pFFT->pTwfIm != NULL)
{
free(pFFT->pTwfIm);
pFFT->pTwfIm = NULL;
}
pFFT->m_numPoints = 0;
pFFT->m_numStages = 0;
#ifdef FFTMsg
printf("\n%d-Point FFT object deleted.\n",pFFT->m_numPoints);
#endif
}
/***************************************************************************/
/* FFT */
/* FAST-FOURIER-TRANSFORMATION s(t) -> S(f) */
/***************************************************************************/
void fft(fft_t *pFFT, fft_float_t *in_re, fft_float_t *in_im)
{
register UINT32 stageCnt, k, kk, twf, numNodesPerStage, nodeCnt, opsPerNode, opCnt, i, i2, j;
fft_float_t tempr, tempi, s, c;
#ifdef FFTMsg
UINT32 numMul, numAdd;
printf ("\nStart %d-Point FFT.\n\n",pFFT->m_numPoints);
#endif
// Do the bit reversal
i2 = pFFT->m_numPoints >> 1;
j = 0;
for (i=0; i <pFFT->m_numPoints-1;i++) {
if (i < j) {
tempr = in_re[i];
tempi = in_im[i];
in_re[i] = in_re[j];
in_im[i] = in_im[j];
in_re[j] = tempr;
in_im[j] = tempi;
}
k = i2;
while (k <= j) {
j = j-k;
k >>= 1;
}
j = j+k;
}
// Calculate the FFT
#ifdef FFTMsg
numMul =0;
numAdd =0;
#endif
numNodesPerStage = pFFT->m_numPoints;
for (stageCnt=1; stageCnt <= pFFT->m_numStages; stageCnt++) {
k=0;
numNodesPerStage = numNodesPerStage/2;
opsPerNode = pFFT->m_numPoints/(2*numNodesPerStage);
#ifdef FFTMsg
printf ("Processing STAGE # %2d",stageCnt);
#endif
for (nodeCnt=1; nodeCnt <= numNodesPerStage; nodeCnt++) {
twf = 0;
opCnt = opsPerNode;
while (opCnt) {
#ifdef FFTMsg
numMul++;
numAdd +=2;
#endif
c = pFFT->pTwfRe[twf];
s = - pFFT->pTwfIm[twf];
kk = k + opsPerNode;
tempr = (c * in_re[kk]) - (s * in_im[kk]);
tempi = (s * in_re[kk]) + (c * in_im[kk]);
in_re[kk] = in_re[k] - tempr;
in_im[kk] = in_im[k] - tempi;
in_re[k] = in_re[k] + tempr;
in_im[k] = in_im[k] + tempi;
k++;
opCnt--;
twf += numNodesPerStage;
}
k += opsPerNode;
}
#ifdef FFTMsg
printf (" ... finished.\n");
#endif
}
#ifdef FFTMsg
printf ("\n%d-Point FFT is completed.\n",pFFT->m_numPoints);
printf ("Total number of complex additions = %u\n",numAdd);
printf ("Total number of complex multiplications = %u\n",numMul);
#endif
}
/***************************************************************************/
/* IFFT */
/* INVERSE-FAST-FOURIER-TRANSFORMATION S(f) -> s(t) */
/***************************************************************************/
void ifft(fft_t *pFFT, fft_float_t *in_re, fft_float_t *in_im)
{
register UINT32 stageCnt, k, kk, twf, numNodesPerStage, nodeCnt, opsPerNode, opCnt, i, i2, j;
fft_float_t tempr, tempi, s, c;
#ifdef FFTMsg
printf ("\nStart %d-Point IFFT.\n\n",pFFT->m_numPoints);
#endif
// Do the bit reversal
i2 = pFFT->m_numPoints >> 1;
j = 0;
for (i=0; i <pFFT->m_numPoints-1;i++) {
if (i < j) {
tempr = in_re[i];
tempi = in_im[i];
in_re[i] = in_re[j];
in_im[i] = in_im[j];
in_re[j] = tempr;
in_im[j] = tempi;
}
k = i2;
while (k <= j) {
j = j-k;
k >>= 1;
}
j = j+k;
}
// Calculate the IFFT
numNodesPerStage = pFFT->m_numPoints;
for (stageCnt=1; stageCnt <= pFFT->m_numStages; stageCnt++) {
k=0;
numNodesPerStage = numNodesPerStage/2;
opsPerNode = pFFT->m_numPoints/(2*numNodesPerStage);
#ifdef FFTMsg
printf ("Processing STAGE # %2d",stageCnt);
#endif
for (nodeCnt=1; nodeCnt <= numNodesPerStage; nodeCnt++) {
twf = 0;
opCnt = opsPerNode;
while (opCnt) {
c = pFFT->pTwfRe[twf];
s = pFFT->pTwfIm[twf];
kk = k + opsPerNode;
tempr = (c * in_re[kk]) - (s * in_im[kk]);
tempi = (s * in_re[kk]) + (c * in_im[kk]);
in_re[kk] = in_re[k] - tempr;
in_im[kk] = in_im[k] - tempi;
in_re[k] = in_re[k] + tempr;
in_im[k] = in_im[k] + tempi;
k++;
opCnt--;
twf += numNodesPerStage;
}
k += opsPerNode;
}
#ifdef FFTMsg
printf (" ... finished.\n");
#endif
}
#ifdef FFTMsg
printf ("\n%d-Point IFFT is completed.\n",pFFT->m_numPoints);
#endif
}
/***************************************************************************/
/* FFTS Autoscale Version X = 1/N * FFT{x}
/* FAST-FOURIER-TRANSFORMATION s(t) -> S(f)
/***************************************************************************/
void ffts(fft_t *pFFT, fft_float_t *in_re, fft_float_t *in_im)
{
register UINT32 stageCnt, k, kk, twf, numNodesPerStage, nodeCnt, opsPerNode, opCnt, i, i2, j;
fft_float_t tempr, tempi, s, c;
#ifdef FFTMsg
UINT32 numMul, numAdd;
printf ("\nStart %d-Point FFT.\n\n",pFFT->m_numPoints);
#endif
// Do the bit reversal
i2 = pFFT->m_numPoints >> 1;
j = 0;
for (i=0; i <pFFT->m_numPoints-1;i++) {
if (i < j) {
tempr = in_re[i];
tempi = in_im[i];
in_re[i] = in_re[j];
in_im[i] = in_im[j];
in_re[j] = tempr;
in_im[j] = tempi;
}
k = i2;
while (k <= j) {
j = j-k;
k >>= 1;
}
j = j+k;
}
// Calculate the FFT
#ifdef FFTMsg
numMul =0;
numAdd =0;
#endif
numNodesPerStage = pFFT->m_numPoints;
for (stageCnt=1; stageCnt <= pFFT->m_numStages; stageCnt++) {
k=0;
numNodesPerStage = numNodesPerStage/2;
opsPerNode = pFFT->m_numPoints/(2*numNodesPerStage);
#ifdef FFTMsg
printf ("Processing STAGE # %2d",stageCnt);
#endif
for (nodeCnt=1; nodeCnt <= numNodesPerStage; nodeCnt++) {
twf = 0;
opCnt = opsPerNode;
while (opCnt) {
#ifdef FFTMsg
numMul++;
numAdd +=2;
#endif
c = pFFT->pTwfRe[twf];
s = - pFFT->pTwfIm[twf];
kk = k + opsPerNode;
tempr = (c * in_re[kk]) - (s * in_im[kk]);
tempi = (s * in_re[kk]) + (c * in_im[kk]);
in_re[kk] = (in_re[k] - tempr)/2;
in_im[kk] = (in_im[k] - tempi)/2;
in_re[k] = (in_re[k] + tempr)/2;
in_im[k] = (in_im[k] + tempi)/2;
k++;
opCnt--;
twf += numNodesPerStage;
}
k += opsPerNode;
}
#ifdef FFTMsg
printf (" ... finished.\n");
#endif
}
#ifdef FFTMsg
printf ("\n%d-Point FFT is completed.\n",pFFT->m_numPoints);
printf ("Total number of complex additions = %u\n",numAdd);
printf ("Total number of complex multiplications = %u\n",numMul);
#endif
}
/***************************************************************************/
/* TWIDDLE-FAKTOR-TABLE */
/* Erstellt Twiddle-Faktor-Tabelle von WN^0 bis WN^N/2 */
/***************************************************************************/
void FFTCalcTwiddleTable (fft_float_t *pRealData, fft_float_t *pImagData, UINT32 numPoints)
{
UINT32 i, size;
fft_float_t arg1;
size = numPoints/2;
arg1 = (fft_float_t)(2.0 * PI / (fft_float_t)numPoints );
#ifdef FFTMsg
printf ("Calculating and initializing Twiddle-Table (%d kB)...",size*sizeof(COMPLEX));
#endif
for (i = 0; i < size; i++) {
pRealData[i] = (fft_float_t)cos(arg1* (fft_float_t)i);
pImagData[i] = (fft_float_t)sin(arg1* (fft_float_t)i);
}
#ifdef FFTMsg
printf(" completed.\n");
#endif
}
/***************************************************************************/
/* MODULUS */
/* Berechnet den Betrag einer komplexen Zahl */
/***************************************************************************/
void Modulus(fft_float_t *pRealData, fft_float_t *pImagData, UINT32 N)
{
UINT32 i;
fft_float_t temp;
for (i = 0; i < N; i++) {
temp = (fft_float_t)sqrt( pRealData[i] * pRealData[i] + pImagData[i] * pImagData[i]);
pImagData[i] = (fft_float_t)atan2(pImagData[i],pRealData[i]);
pRealData[i] = temp;
}
}
/***************************************************************************/
/* Hanning */
/* Legt das Hanningfenster auf die Abtastwerte im Zeitbereich der Groesse N */
/***************************************************************************/
void Hanning(fft_float_t *pRealData, fft_float_t *pImagData, UINT32 N, UINT32 maximum)
{
UINT32 n;
fft_float_t arg;
for (n=0; n < N; n++)
{
arg = (fft_float_t)(2*PI*n /(N-1) - 2*PI*maximum/N);
pRealData[n] = pRealData[n] * (fft_float_t)(1 + cos(arg)) /2;
pImagData[n] = pImagData[n] * (fft_float_t)(1 + cos(arg)) /2;
}
}
/***************************************************************************/
/* HANNING_K */
/* Gibt einen Faktor k in Abhängigkeit von n bezogen auf N zurück. */
/***************************************************************************/
fft_float_t hanning_k(UINT32 n, UINT32 N)
{
fft_float_t arg;
arg = (fft_float_t)(2 * PI / (fft_float_t)N);
return ((1-(fft_float_t)cos(arg*(fft_float_t)n))/2);
}
/***************************************************************************/
/* GAUSS_K */
/* Gibt einen Faktor k in Abhängigkeit von n bezogen auf N zurück. */
/***************************************************************************/
fft_float_t gauss_k(UINT32 n, UINT32 m, UINT32 s)
{
fft_float_t arg;
arg = (fft_float_t)((n-m)*(n-m)/(2*s*s));
return (fft_float_t)(exp(-arg));
}
/***************************************************************************/
/* Normalize) */
/* */
/***************************************************************************/
void Scale(fft_float_t *pRealData, fft_float_t *pImagData, fft_float_t scaleFactor, UINT32 N)
{
UINT32 k;
// Scaling Data
for (k=0; k < N; k++)
{
pRealData[k] *= (fft_float_t)scaleFactor;
pImagData[k] *= (fft_float_t)scaleFactor;
}
}
/***************************************************************************/
/* BiPower()
/* checks if N is power of 2
/***************************************************************************/
UINT32 IsPowerOfTwo(UINT32 N)
{
UINT32 i, iterations;
iterations = sizeof(UINT32)*8;
if (N == 0)
return AV_E_FALSE;
for (i=0; i <iterations; i++)
{
if (N == (UINT32)(2 << i))
return AV_E_OK;
}
return FFT_ERROR;
}
/***************************************************************************/
/* FFT2Dinit()
/*
/***************************************************************************/
void FFT2Dinit(fft2_t *pFFT, UINT32 Nx, UINT32 Ny)
{
pFFT->pXfft = NULL;
pFFT->pYfft = NULL;
pFFT->m_Nx = Nx;
pFFT->m_Ny = Ny;
pFFT->pXfft = (fft_t*)malloc(sizeof(pFFT->pXfft));
FFTinit(pFFT->pXfft, Nx);
if (Nx == Ny)
pFFT->pYfft = pFFT->pXfft;
else
{
pFFT->pYfft = (fft_t*)malloc(sizeof(pFFT->pYfft));
FFTinit(pFFT->pYfft, Ny);
}
}
/***************************************************************************/
/* FFT2Dfree()
/*
/***************************************************************************/
void FFT2Dfree(fft2_t *pFFT)
{
if (pFFT->pXfft != NULL)
{
FFTfree(pFFT->pXfft);
pFFT->pXfft = NULL;
}
if (pFFT->pYfft != NULL)
{
FFTfree(pFFT->pYfft);
pFFT->pYfft = NULL;
}
pFFT->m_Nx = 0;
pFFT->m_Ny = 0;
}
/***************************************************************************/
/* fft2d()
/*
/***************************************************************************/
void fft2d(fft2_t *pFFT, fft_float_t **ppReal, fft_float_t **ppImag)
{
fft_float_t *pTempr, *pTempi;
UINT32 i, row;
pTempr = (fft_float_t*)malloc(pFFT->m_Ny*sizeof(fft_float_t));
pTempi = (fft_float_t*)malloc(pFFT->m_Ny*sizeof(fft_float_t));
/* Transform ROWs */
for (row=0; row <pFFT->m_Ny; row++)
fft(pFFT->pXfft, ppReal[row], ppImag[row]);
/* Transform COLs */
for (row=0; row <pFFT->m_Nx; row++)
{
for(i=0; i <pFFT->m_Ny; i++)
{
pTempr[i] = ppReal[i][row];
pTempi[i] = ppImag[i][row];
}
fft(pFFT->pYfft, pTempr, pTempi);
for(i=0; i <pFFT->m_Ny; i++)
{
ppReal[i][row] = pTempr[i];
ppImag[i][row] = pTempi[i];
}
}
free(pTempr);
free(pTempi);
}
/***************************************************************************/
/* ifft2d()
/*
/***************************************************************************/
void ifft2d(fft2_t *pFFT, fft_float_t **ppReal, fft_float_t **ppImag)
{
fft_float_t *pTempr, *pTempi;
UINT32 i, row;
pTempr = (fft_float_t*)malloc(pFFT->m_Ny*sizeof(fft_float_t));
pTempi = (fft_float_t*)malloc(pFFT->m_Ny*sizeof(fft_float_t));
/* Transform COLs */
for (row=0; row <pFFT->m_Nx; row++)
{
for(i=0; i <pFFT->m_Ny; i++)
{
pTempr[i] = ppReal[i][row];
pTempi[i] = ppImag[i][row];
}
ifft(pFFT->pYfft, pTempr, pTempi);
for(i=0; i <pFFT->m_Ny; i++)
{
ppReal[i][row] = pTempr[i];
ppImag[i][row] = pTempi[i];
}
}
free(pTempr);
free(pTempi);
/* Transform ROWs */
for (row=0; row <pFFT->m_Ny; row++)
ifft(pFFT->pXfft, ppReal[row], ppImag[row]);
}
/***************************************************************************/
/* DFT()
/* DISKRETE FOURIER-TRANSFORMATION S(f) -> s(t) */
/***************************************************************************/
void DFT(fft_float_t *pDataR, fft_float_t *pDataI, UINT32 N)
{
fft_float_t *a, *b;
fft_float_t phi, c, s;
UINT32 i, j;
a = (fft_float_t*)malloc(N * sizeof(fft_float_t));
b = (fft_float_t*)malloc(N * sizeof(fft_float_t));
memcpy((fft_float_t*)a,(fft_float_t*)pDataR,N * sizeof(fft_float_t));
memcpy((fft_float_t*)b,(fft_float_t*)pDataI,N * sizeof(fft_float_t));
for (i = 0; i < N; i++)
{
pDataR[i] = 0;
pDataI[i] = 0;
for (j = 0; j < N; j++)
{
phi = (fft_float_t)(2*i*j*PI/N);
c = (fft_float_t)(0.5*cos(phi));
s = (fft_float_t)(0.5*sin(phi));
pDataR[i] += c*a[j] - s*b[j];
pDataI[i] += s*a[j] + c*b[j];
}
}
free(a);
free(b);
}
/***************************************************************************/
/* IDFT()
/* INVERSE DISKRETE FOURIER-TRANSFORMATION S(f) -> s(t) */
/***************************************************************************/
void IDFT(fft_float_t *pDataR, fft_float_t *pDataI, UINT32 N)
{
fft_float_t *a, *b;
fft_float_t phi, c, s;
UINT32 i, j;
a = (fft_float_t*)malloc(N * sizeof(fft_float_t));
b = (fft_float_t*)malloc(N * sizeof(fft_float_t));
memcpy((fft_float_t*)a,(fft_float_t*)pDataR,N * sizeof(fft_float_t));
memcpy((fft_float_t*)b,(fft_float_t*)pDataI,N * sizeof(fft_float_t));
for (i = 0; i < N; i++)
{
pDataR[i] = 0;
pDataI[i] = 0;
for (j = 0; j < N; j++)
{
phi = (fft_float_t)(2*i*j*PI/N);
c = (fft_float_t) (0.5*(fft_float_t)cos(phi));
s = (fft_float_t)(-0.5*(fft_float_t)sin(phi));
pDataR[i] += c*a[j] - s*b[j];
pDataI[i] += s*a[j] + c*b[j];
}
}
free(a);
free(b);
}
Reference in New Issue View Git Blame Copy Permalink