IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 7, 2007, 6:17:58 PM (19 years ago)
Author:
Paul Price
Message:

Reworked image convolution (split direct and FFT methods into separate functions) and FFT functions. psVectorFFT needs a bit more work to clean up the N/2 issues from FFTW's r2c and c2r.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/fft/psVectorFFT.c

    r11668 r11703  
    33 *  @brief Contains FFT transform related functions for psVector
    44 *
     5 *  @author Paul Price, IfA
    56 *  @author Robert DeSonia, MHPCC
    67 *
    7  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-02-06 21:36:09 $
     8 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2007-02-08 04:17:58 $
    910 *
    1011 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    1213
    1314#ifdef HAVE_CONFIG_H
    14 # include "config.h"
     15#include "config.h"
    1516#endif
    1617
    … …  
    2122#include <fftw3.h>
    2223
    23 #include "psVectorFFT.h"
     24#include "psAssert.h"
    2425#include "psError.h"
    2526#include "psMemory.h"
    2627#include "psLogMsg.h"
     28#include "psConstants.h"
     29#include "psVectorFFT.h"
     30
     31#define FFTW_PLAN_RIGOR FFTW_ESTIMATE   // How rigorous the FFTW planning is
     32
     33static psBool fftwWisdomImported = false; // Has the system wisdom been imported?
     34
     35bool psVectorForwardFFT(psVector **real, psVector **imag, const psVector *in)
     36{
     37    PS_ASSERT_VECTOR_NON_NULL(in, false);
     38    PS_ASSERT_VECTOR_TYPE(in, PS_TYPE_F32, false);
     39    PS_ASSERT_PTR_NON_NULL(real, false);
     40    PS_ASSERT_PTR_NON_NULL(imag, false);
     41
     42    // Make sure the system-level wisdom information is imported.
     43    if (!fftwWisdomImported) {
     44        fftwf_import_system_wisdom();
     45        fftwWisdomImported = true;
     46    }
     47
     48    long num = in->n;                   // Number of elements
     49
     50    // Do the FFT
     51    fftwf_complex *out = fftwf_malloc(num * sizeof(fftwf_complex)); // Output data
     52    fftwf_plan plan = fftwf_plan_dft_r2c_1d(num, in->data.F32, out, FFTW_PLAN_RIGOR);
     53    fftwf_execute(plan);
     54    fftwf_destroy_plan(plan);
     55
     56    // Pull the real and imaginary parts out
     57    *real = psVectorRecycle(*real, num/2 + 1, PS_TYPE_F32);
     58    *imag = psVectorRecycle(*imag, num/2 + 1, PS_TYPE_F32);
     59    for (int i = 0; i < num; i++) {
     60#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
     61        // C99 complex support
     62        (*real)->data.F32[i] = creal(out[i]);
     63        (*imag)->data.F32[i] = cimag(out[i]);
     64#else
     65        // FFTW's backup complex support
     66        (*real)->data.F32[i] = out[i][0];
     67        (*imag)->data.F32[i] = out[i][1];
     68#endif
     69    }
     70    fftwf_free(out);
     71
     72    return true;
     73}
     74
     75bool psVectorBackwardFFT(psVector **out, const psVector *real, const psVector *imag)
     76{
     77    PS_ASSERT_VECTOR_NON_NULL(real, false);
     78    PS_ASSERT_VECTOR_TYPE(real, PS_TYPE_F32, false);
     79    PS_ASSERT_VECTOR_NON_NULL(imag, false);
     80    PS_ASSERT_VECTOR_TYPE(imag, PS_TYPE_F32, false);
     81    PS_ASSERT_VECTORS_SIZE_EQUAL(real, imag, false);
     82    PS_ASSERT_PTR_NON_NULL(out, false);
     83
     84    // Make sure the system-level wisdom information is imported.
     85    if (!fftwWisdomImported) {
     86        fftwf_import_system_wisdom();
     87        fftwWisdomImported = true;
     88    }
     89
     90    long num = real->n;                 // Number of elements
     91
     92    // Stuff the real and imaginary parts in
     93    fftwf_complex *in = fftwf_malloc(num * sizeof(fftwf_complex)); // Input data
     94    for (int i = 0; i < num; i++) {
     95#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
     96        // C99 complex support
     97        in[i] = real->data.F32[i] + imag->data.F32[i] * I;
     98#else
     99        // FFTW's backup complex support
     100        in[i][0] = real->data.F32[i];
     101        in[i][1] = imag->data.F32[i];
     102#endif
     103    }
    27104
    28105
     106    // Do the FFT
     107    *out = psVectorRecycle(*out, 2 * num, PS_TYPE_F32);
     108    fftwf_plan plan = fftwf_plan_dft_c2r_1d(2 * num, in, (*out)->data.F32, FFTW_PLAN_RIGOR);
     109    fftwf_execute(plan);
     110    fftwf_destroy_plan(plan);
    29111
    30 #define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE
     112    fftwf_free(in);
    31113
    32 static bool p_fftwWisdomImported = false;
     114    return true;
     115}
    33116
    34 psVector* psVectorFFT(psVector* out, const psVector* in, psFFTFlags direction)
     117psVector *psVectorPowerSpectrum(psVector* out, const psVector* in)
    35118{
    36     psU32 numElements;
    37     psElemType type;
    38     fftwf_plan plan;
     119    PS_ASSERT_VECTOR_NON_NULL(in, NULL);
     120    PS_ASSERT_VECTOR_TYPE(in, PS_TYPE_F32, NULL);
    39121
    40     /* got good image data? */
    41     if (in == NULL) {
    42         psFree(out);
     122    psVector *real = NULL;              // Real component of FFT
     123    psVector *imag = NULL;              // Imaginary component of FFT
     124
     125    if (!psVectorForwardFFT(&real, &imag, in)) {
     126        psError(PS_ERR_UNKNOWN, false, "Unable to perform forward FFT.");
    43127        return NULL;
    44128    }
    45129
    46     type = in->type.type;
     130    int num = real->n;                  // Number of elements
    47131
    48     /* make sure the system-level wisdom information is imported. */
    49     if (!p_fftwWisdomImported) {
    50         fftwf_import_system_wisdom();
    51         p_fftwWisdomImported = true;
     132    float norm = 1.0 / PS_SQR(in->n);
     133    for (int i = 0; i < num; i++) {
     134        // Power spectrum is the square of the complex modulus
     135        real->data.F32[i] = norm * (PS_SQR(real->data.F32[i]) + PS_SQR(imag->data.F32[i]));
     136    }
     137    psFree(imag);
     138
     139    return real;
     140}
     141
     142bool psVectorComplexMultiply(psVector **outReal, psVector **outImag,
     143                             const psVector *in1Real, const psVector *in1Imag,
     144                             const psVector *in2Real, const psVector *in2Imag)
     145{
     146    PS_ASSERT_VECTOR_NON_NULL(in1Real, false);
     147    PS_ASSERT_VECTOR_NON_NULL(in1Imag, false);
     148    PS_ASSERT_VECTOR_NON_NULL(in2Real, false);
     149    PS_ASSERT_VECTOR_NON_NULL(in2Imag, false);
     150    PS_ASSERT_VECTOR_TYPE(in1Real, PS_TYPE_F32, false);
     151    PS_ASSERT_VECTOR_TYPE(in1Imag, PS_TYPE_F32, false);
     152    PS_ASSERT_VECTOR_TYPE(in2Real, PS_TYPE_F32, false);
     153    PS_ASSERT_VECTOR_TYPE(in2Imag, PS_TYPE_F32, false);
     154    PS_ASSERT_VECTORS_SIZE_EQUAL(in1Imag, in1Real, false);
     155    PS_ASSERT_VECTORS_SIZE_EQUAL(in2Real, in1Real, false);
     156    PS_ASSERT_VECTORS_SIZE_EQUAL(in2Imag, in1Real, false);
     157    PS_ASSERT_PTR_NON_NULL(outReal, false);
     158    PS_ASSERT_PTR_NON_NULL(outImag, false);
     159    if (*outReal) {
     160        PS_ASSERT_VECTOR_NON_NULL(*outReal, false);
     161        PS_ASSERT_VECTOR_NON_NULL(*outImag, false);
     162    } else if (*outImag) {
     163        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     164                "If the output real part is provided, the output imaginary part must also be provided.");
     165        return false;
    52166    }
    53167
    54     numElements = in->n;
     168    int num = in1Real->n;               // Number of elements
    55169
    56     out = psVectorCopy(out, in, PS_TYPE_C32);
    57     out->n = numElements;
     170    *outReal = psVectorRecycle(*outReal, num, PS_TYPE_F32);
     171    *outImag = psVectorRecycle(*outImag, num, PS_TYPE_F32);
    58172
    59     if ((direction & PS_FFT_FORWARD) != 0) {
    60         plan = fftwf_plan_dft_1d(numElements,
    61                                  (fftwf_complex *) out->data.C32,
    62                                  (fftwf_complex *) out->data.C32, FFTW_FORWARD, P_FFTW_PLAN_RIGOR);
    63     } else if ((direction & PS_FFT_REVERSE) != 0) {
    64         plan = fftwf_plan_dft_1d(numElements,
    65                                  (fftwf_complex *) out->data.C32,
    66                                  (fftwf_complex *) out->data.C32, FFTW_BACKWARD, P_FFTW_PLAN_RIGOR);
    67     } else {
    68         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    69                 _("Must specify the direction as either PS_FFT_FORWARD or PS_FFT_REVERSE."));
    70         psFree(out);
    71         return NULL;
     173    for (int i = 0; i < num; i++) {
     174        // (a + bi) * (c + di) = (ac - bd) + (bc + ad)i
     175        (*outReal)->data.F32[i] = in1Real->data.F32[i] * in2Real->data.F32[i] -
     176            in1Imag->data.F32[i] * in2Imag->data.F32[i];
     177        (*outImag)->data.F32[i] = in1Imag->data.F32[i] * in2Real->data.F32[i] +
     178            in1Real->data.F32[i] * in2Imag->data.F32[i];
    72179    }
    73180
    74     /* check if a plan exists now */
    75     if (plan == NULL) {
    76         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    77                 _("Could not create a valid FFT plan to perform the transform."));
    78         psFree(out);
    79         return NULL;
    80     }
    81 
    82     /* finally, call FFTW with the plan made above */
    83     fftwf_execute(plan);
    84 
    85     fftwf_destroy_plan(plan);
    86 
    87     if ((direction & PS_FFT_REAL_RESULT) != 0) {
    88         for (psS32 i = 0; i < numElements; i++) {
    89             out->data.F32[i] = out->data.C32[i];
    90         }
    91         out->type.type = PS_TYPE_F32;
    92         out->data.U8 = psRealloc(out->data.U8,PSELEMTYPE_SIZEOF(PS_TYPE_F32)*out->nalloc);
    93     }
    94 
    95     return out;
     181    return true;
    96182}
    97 
    98 psVector* psVectorReal(psVector* out, const psVector* in)
    99 {
    100     psElemType type;
    101     psU32 numElements;
    102 
    103     if (in == NULL) {
    104         psFree(out);
    105         return NULL;
    106     }
    107 
    108     type = in->type.type;
    109     numElements = in->n;
    110 
    111     /* if not a complex number, this is logically just a copy */
    112     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    113         // Warn user, as this is probably not expected
    114         psLogMsg(__func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
    115                  "Just a vector copy was performed.");
    116         out = psVectorRecycle(out, numElements, type);
    117         out->n = numElements;
    118         memcpy(out->data.U8, in->data.U8, numElements * PSELEMTYPE_SIZEOF(type));
    119         return out;
    120     }
    121 
    122     if (type == PS_TYPE_C32) {
    123         psF32* outVec;
    124         psC32* inVec = in->data.C32;
    125 
    126         out = psVectorRecycle(out, numElements, PS_TYPE_F32);
    127         out->n = numElements;
    128         outVec = out->data.F32;
    129 
    130         for (psU32 i = 0; i < numElements; i++) {
    131             outVec[i] = crealf(inVec[i]);
    132         }
    133     } else if (type == PS_TYPE_C64) {
    134         psF64* outVec;
    135         psC64* inVec = in->data.C64;
    136 
    137         out = psVectorRecycle(out, numElements, PS_TYPE_F64);
    138         out->n = numElements;
    139         outVec = out->data.F64;
    140 
    141         for (psU32 i = 0; i < numElements; i++) {
    142             outVec[i] = creal(inVec[i]);
    143         }
    144     } else {
    145         char* typeStr;
    146         PS_TYPE_NAME(typeStr,type);
    147         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    148                 _("Specified psVector type, %s, is not supported."),
    149                 typeStr);
    150         psFree(out);
    151         return NULL;
    152     }
    153 
    154     return out;
    155 }
    156 
    157 psVector* psVectorImaginary(psVector* out, const psVector* in)
    158 {
    159     psElemType type;
    160     psU32 numElements;
    161 
    162     if (in == NULL) {
    163         psFree(out);
    164         return NULL;
    165     }
    166 
    167     type = in->type.type;
    168     numElements = in->n;
    169 
    170     /* if not a complex number, this is logically just zeroed image of same size */
    171     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    172         // Warn user, as this is probably not expected
    173         psLogMsg(__func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
    174                  "A zeroed vector was returned.");
    175         out = psVectorRecycle(out, numElements, type);
    176         out->n = numElements;
    177         memset(out->data.U8, 0, PSELEMTYPE_SIZEOF(type) * numElements);
    178         return out;
    179     }
    180 
    181     if (type == PS_TYPE_C32) {
    182         psF32* outVec;
    183         psC32* inVec = in->data.C32;
    184 
    185         out = psVectorRecycle(out, numElements, PS_TYPE_F32);
    186         out->n = numElements;
    187         outVec = out->data.F32;
    188 
    189         for (psU32 i = 0; i < numElements; i++) {
    190             outVec[i] = cimagf(inVec[i]);
    191         }
    192     } else if (type == PS_TYPE_C64) {
    193         psF64* outVec;
    194         psC64* inVec = in->data.C64;
    195 
    196         out = psVectorRecycle(out, numElements, PS_TYPE_F64);
    197         out->n = numElements;
    198         outVec = out->data.F64;
    199 
    200         for (psU32 i = 0; i < numElements; i++) {
    201             outVec[i] = cimag(inVec[i]);
    202         }
    203     } else {
    204         char* typeStr;
    205         PS_TYPE_NAME(typeStr,type);
    206         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    207                 _("Specified psVector type, %s, is not supported."),
    208                 typeStr);
    209         psFree(out);
    210         return NULL;
    211     }
    212 
    213     return out;
    214 }
    215 
    216 psVector* psVectorComplex(psVector* out, const psVector* real, const psVector* imag)
    217 {
    218     psElemType type;
    219     psU32 numElements;
    220 
    221     if (real == NULL || imag == NULL) {
    222         psFree(out);
    223         return NULL;
    224     }
    225 
    226     type = real->type.type;
    227     if (real->n < imag->n) {
    228         numElements = real->n;
    229     } else {
    230         numElements = imag->n;
    231     }
    232 
    233     if (imag->type.type != type) {
    234         char* typeStrReal;
    235         char* typeStrImag;
    236         PS_TYPE_NAME(typeStrReal,type);
    237         PS_TYPE_NAME(typeStrImag,imag->type.type);
    238         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    239                 _("Real psVector type, %s, and imaginary psVector type, %s, must be the same."),
    240                 typeStrReal,typeStrImag);
    241         psFree(out);
    242         return NULL;
    243     }
    244 
    245     if (type == PS_TYPE_F32) {
    246         psC32* outVec;
    247         psF32* realVec = real->data.F32;
    248         psF32* imagVec = imag->data.F32;
    249 
    250         out = psVectorRecycle(out, numElements, PS_TYPE_C32);
    251         out->n = numElements;
    252         outVec = out->data.C32;
    253 
    254         for (psU32 i = 0; i < numElements; i++) {
    255             outVec[i] = realVec[i] + I * imagVec[i];
    256         }
    257     } else if (type == PS_TYPE_F64) {
    258         psC64* outVec;
    259         psF64* realVec = real->data.F64;
    260         psF64* imagVec = imag->data.F64;
    261 
    262         out = psVectorRecycle(out, numElements, PS_TYPE_C64);
    263         out->n = numElements;
    264         outVec = out->data.C64;
    265 
    266         for (psU32 i = 0; i < numElements; i++) {
    267             outVec[i] = realVec[i] + I * imagVec[i];
    268         }
    269     } else {
    270         char* typeStr;
    271         PS_TYPE_NAME(typeStr,type);
    272         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    273                 _("Input psVector type, %s, is required to be either psF32 or psF64."),
    274                 typeStr);
    275         psFree(out);
    276         return NULL;
    277     }
    278 
    279     return out;
    280 }
    281 
    282 psVector* psVectorConjugate(psVector* out, const psVector* in)
    283 {
    284     psElemType type;
    285     psU32 numElements;
    286 
    287     if (in == NULL) {
    288         psFree(out);
    289         return NULL;
    290     }
    291 
    292     type = in->type.type;
    293     numElements = in->n;
    294 
    295     /* if not a complex number, this is logically just a image copy */
    296     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    297         // Warn user, as this is probably not expected
    298         psLogMsg(__func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
    299                  "Vector copy was performed instead.");
    300 
    301         out = psVectorRecycle(out, numElements, type);
    302         out->n = numElements;
    303         memcpy(out->data.U8, in->data.U8, PSELEMTYPE_SIZEOF(type) * numElements);
    304         return out;
    305     }
    306 
    307     if (type == PS_TYPE_C32) {
    308         psC32* outVec;
    309         psC32* inVec = in->data.C32;
    310 
    311         out = psVectorRecycle(out, numElements, PS_TYPE_C32);
    312         out->n = numElements;
    313         outVec = out->data.C32;
    314 
    315         for (psU32 i = 0; i < numElements; i++) {
    316             outVec[i] = crealf(inVec[i]) - I * cimagf(inVec[i]);
    317         }
    318     } else if (type == PS_TYPE_C64) {
    319         psC64* outVec;
    320         psC64* inVec = in->data.C64;
    321 
    322         out = psVectorRecycle(out, numElements, PS_TYPE_C64);
    323         out->n = numElements;
    324         outVec = out->data.C64;
    325 
    326         for (psU32 i = 0; i < numElements; i++) {
    327             outVec[i] = creal(inVec[i]) - I * cimag(inVec[i]);
    328         }
    329     } else {
    330         char* typeStr;
    331         PS_TYPE_NAME(typeStr,type);
    332         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    333                 _("Input psVector type, %s, is required to be either psC32 or psC64."),
    334                 typeStr);
    335         psFree(out);
    336         return NULL;
    337     }
    338 
    339     return out;
    340 }
    341 
    342 psVector* psVectorPowerSpectrum(psVector* out, const psVector* in)
    343 {
    344     psElemType type;
    345     psU32 outNumElements;
    346     psU32 inNumElements;
    347     psU32 inHalfNumElements;
    348     psU32 inNumElementsSquared;
    349 
    350     if (in == NULL) {
    351         psFree(out);
    352         return NULL;
    353     }
    354 
    355     type = in->type.type;
    356     inNumElements = in->n;
    357     inNumElementsSquared = inNumElements * inNumElements;
    358     inHalfNumElements = inNumElements / 2;
    359     outNumElements = inHalfNumElements + 1;
    360 
    361     if (type == PS_TYPE_C32) {
    362         psF32* outVec;
    363         psC32* inVec = in->data.C32;
    364         psF32 inAbs1;
    365         psF32 inAbs2;
    366 
    367         out = psVectorRecycle(out, outNumElements, PS_TYPE_F32);
    368         out->n = outNumElements;
    369         outVec = out->data.F32;
    370 
    371         // from ADD: P_0 = |C_0|^2/N^2
    372         inAbs1 = cabsf(inVec[0]);
    373         outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared;
    374 
    375         // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
    376         for (psU32 i = 1; i < inHalfNumElements; i++) {
    377             inAbs1 = cabsf(inVec[i]);
    378             inAbs2 = cabsf(inVec[inNumElements - i]);
    379             outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared;
    380         }
    381 
    382         // from ADD: P_N/2 = |C_N/2|^2/N^2
    383         inAbs1 = cabsf(inVec[inHalfNumElements]);
    384         outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared;
    385     } else if (type == PS_TYPE_C64) {
    386         psF64* outVec;
    387         psC64* inVec = in->data.C64;
    388         psF64 inAbs1;
    389         psF64 inAbs2;
    390 
    391         out = psVectorRecycle(out, outNumElements, PS_TYPE_F64);
    392         out->n = outNumElements;
    393         outVec = out->data.F64;
    394 
    395         // from ADD: P_0 = |C_0|^2/N^2
    396         inAbs1 = cabs(inVec[0]);
    397         outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared;
    398 
    399         // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
    400         for (psU32 i = 1; i < inHalfNumElements; i++) {
    401             inAbs1 = cabs(inVec[i]);
    402             inAbs2 = cabs(inVec[inNumElements - i]);
    403             outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared;
    404         }
    405 
    406         // from ADD: P_N/2 = |C_N/2|^2/N^2
    407         inAbs1 = cabs(inVec[inHalfNumElements]);
    408         outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared;
    409     } else {
    410         char* typeStr;
    411         PS_TYPE_NAME(typeStr,type);
    412         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    413                 _("Input psVector type, %s, is required to be either psC32 or psC64."),
    414                 typeStr);
    415         psFree(out);
    416         return NULL;
    417     }
    418 
    419     return out;
    420 
    421 }
Note: See TracChangeset for help on using the changeset viewer.