changeset 772:59df354b99cc

Add some utility functions in the FFT module.
author Matti Hamalainen <ccr@tnsp.org>
date Fri, 10 May 2013 12:34:45 +0300
parents 71af31139e46
children 95d9bd557a54
files dmfft.c dmfft.h
diffstat 2 files changed, 81 insertions(+), 10 deletions(-) [+]
line wrap: on
line diff
--- a/dmfft.c	Thu May 09 19:58:06 2013 +0300
+++ b/dmfft.c	Fri May 10 12:34:45 2013 +0300
@@ -111,8 +111,8 @@
     /*
      *  Massage output to get the output for a real input sequence.
      */
-    int *br1 = ctx->breversed + 1;
-    int *br2 = ctx->breversed + ctx->npoints - 1;
+    const int *br1 = ctx->breversed + 1;
+    const int *br2 = ctx->breversed + ctx->npoints - 1;
 
     while (br1 <= br2)
     {
@@ -151,17 +151,20 @@
 }
 
 
-int dmConvertFFTtoFreqDomain(DMFFTContext *ctx, DMFFTType *buffer, DMFFTType *real, DMFFTType *imag)
+int dmConvertFFTtoFreqDomain(DMFFTContext *ctx, DMFFTType *buffer,
+    DMFFTType *real, DMFFTType *imag)
 {
     int i;
     if (ctx == NULL || buffer == NULL || real == NULL || imag == NULL)
         return DMERR_NULLPTR;
 
     // Convert from bitreversed form to normal frequency domain
+    const int *breversed = ctx->breversed;
+
     for (i = 1; i < ctx->npoints; i++)
     {
-        real[i] = buffer[ctx->breversed[i]  ];
-        imag[i] = buffer[ctx->breversed[i]+1];
+        real[i] = buffer[breversed[i]  ];
+        imag[i] = buffer[breversed[i]+1];
     }
 
     // Make the ends meet
@@ -174,16 +177,79 @@
 }
 
 
+#define DM_PW(a, b, c) sqrt(a * a + b * b) * c
+
+
+int dmConvertFFTtoFreqAndPower(DMFFTContext *ctx, DMFFTType *buffer,
+    DMFFTType *real, DMFFTType *imag, DMFFTType *power, const DMFFTType scale)
+{
+    int i;
+    if (ctx == NULL || buffer == NULL || real == NULL || imag == NULL || power == NULL)
+        return DMERR_NULLPTR;
+
+    // Convert from bitreversed form to normal frequency domain
+    const int *breversed = ctx->breversed;
+
+    for (i = 1; i < ctx->npoints; i++)
+    {
+        DMFFTType fre = real[i] = buffer[breversed[i]  ],
+                  fim = imag[i] = buffer[breversed[i]+1];
+
+        power[i] = DM_PW(fre, fim, scale);
+    }
+
+    // Make the ends meet
+    real[0]             = buffer[0];
+    imag[0]             = 0;
+    power[0]            = DM_PW(buffer[0], 0, scale);
+    
+    real[ctx->npoints]  = buffer[1];
+    imag[ctx->npoints]  = 0;
+    power[ctx->npoints] = DM_PW(buffer[1], 0, scale);
+
+    return DMERR_OK;
+}
+
+
+int dmConvertFFTtoPowerAndSum(DMFFTContext *ctx, DMFFTType *buffer,
+    DMFFTType *power, const DMFFTType pscale, DMFFTType *sum, const DMFFTType sscale)
+{
+    int i;
+    if (ctx == NULL || buffer == NULL || power == NULL || sum == NULL)
+        return DMERR_NULLPTR;
+
+    // Convert from bitreversed form to normal frequency domain
+    const int *breversed = ctx->breversed;
+    DMFFTType psum = 0;
+
+    for (i = 1; i < ctx->npoints; i++)
+    {
+        DMFFTType fre = buffer[breversed[i]  ],
+                  fim = buffer[breversed[i]+1],
+                  fpw = sqrt(fre * fre + fim * fim);
+
+        power[i]  = fpw * pscale;
+        psum     += fpw * sscale;
+    }
+
+    // Make the ends meet
+    power[0]            = DM_PW(buffer[0], 0, pscale);
+    power[ctx->npoints] = DM_PW(buffer[1], 0, pscale);
+
+    return DMERR_OK;
+}
+
+
 int dmConvertFFTtoTimeDomain(DMFFTContext *ctx, DMFFTType *buffer, DMFFTType *tdom)
 {
-    int i;
+    int i, n;
     if (ctx == NULL || buffer == NULL || tdom == NULL)
         return DMERR_NULLPTR;
 
-    for (i = 1; i < ctx->npoints; i++)
+    for (n = i = 0; i < ctx->npoints; i++, n += 2)
     {
-        tdom[i*2  ] = buffer[ctx->breversed[i]  ];
-        tdom[i*2+1] = buffer[ctx->breversed[i]+1];
+        tdom[n  ] = buffer[ctx->breversed[i]  ];
+        tdom[n+1] = buffer[ctx->breversed[i]+1];
     }
 
     return DMERR_OK;
--- a/dmfft.h	Thu May 09 19:58:06 2013 +0300
+++ b/dmfft.h	Fri May 10 12:34:45 2013 +0300
@@ -18,7 +18,12 @@
 int   dmInitializeFFT(DMFFTContext *, int);
 void  dmEndFFT(DMFFTContext *);
 int   dmRealFFT(DMFFTContext *, DMFFTType *);
-int   dmConvertFFTtoFreqDomain(DMFFTContext *ctx, DMFFTType *buffer, DMFFTType *real, DMFFTType *imag);
+int   dmConvertFFTtoFreqDomain(DMFFTContext *ctx, DMFFTType *buffer,
+      DMFFTType *real, DMFFTType *imag);
+int   dmConvertFFTtoFreqAndPower(DMFFTContext *ctx, DMFFTType *buffer,
+      DMFFTType *real, DMFFTType *imag, DMFFTType *power, const DMFFTType scale);
+int   dmConvertFFTtoPowerAndSum(DMFFTContext *ctx, DMFFTType *buffer,
+      DMFFTType *power, const DMFFTType pscale, DMFFTType *sum, const DMFFTType sscale);
 int   dmConvertFFTtoTimeDomain(DMFFTContext *ctx, DMFFTType *buffer, DMFFTType *tdom);