Mercurial > hg > dmlib
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);