IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 11703 for trunk/psLib/src/fft


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.

Location:
trunk/psLib/src/fft
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r11668 r11703  
    1 /** @file  psImageFFT.c
    2  *
    3  *  @brief Contains FFT transform related functions for psImage.
    4  *
    5  *  @author Robert DeSonia, MHPCC
    6  *
    7  *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2007-02-06 21:36:09 $
    9  *
    10  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    11  */
     1/// @file  psImageFFT.c
     2///
     3/// @brief Contains FFT transform related functions for psImage.
     4///
     5/// @author Paul Price, IfA
     6/// @author Robert DeSonia, MHPCC
     7///
     8/// @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
     9/// @date $Date: 2007-02-08 04:17:58 $
     10///
     11/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     12///
    1213
    1314#ifdef HAVE_CONFIG_H
    … …  
    2021#include <fftw3.h>
    2122
    22 #include "psImageFFT.h"
     23#include "psAssert.h"
    2324#include "psError.h"
    2425#include "psMemory.h"
    2526#include "psLogMsg.h"
     27#include "psConstants.h"
    2628#include "psImageStructManip.h"
    27 
    28 
    29 
    30 #define PS_FFTW_PLAN_RIGOR FFTW_ESTIMATE
    31 
    32 static bool p_fftwWisdomImported = false;
    33 
    34 psImage* psImageFFT(psImage* out, const psImage* image, psFFTFlags direction)
    35 {
    36     psU32 numCols;
    37     psU32 numRows;
    38     psElemType type;
    39     fftwf_plan plan;
    40 
    41     /* got good image data? */
    42     if (image == NULL) {
    43         psFree(out);
     29#include "psImageFFT.h"
     30
     31
     32#define FFTW_PLAN_RIGOR FFTW_ESTIMATE   // How rigorous the FFTW planning is
     33
     34static psBool fftwWisdomImported = false; // Has the FFTW wisdom been imported yet?
     35
     36bool psImageForwardFFT(psImage **real, psImage **imag, const psImage *in)
     37{
     38    PS_ASSERT_IMAGE_NON_NULL(in, false);
     39    PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, false);
     40    PS_ASSERT_PTR_NON_NULL(real, false);
     41    PS_ASSERT_PTR_NON_NULL(imag, false);
     42
     43    // Make sure the system-level wisdom information is imported.
     44    if (!fftwWisdomImported) {
     45        fftwf_import_system_wisdom();
     46        fftwWisdomImported = true;
     47    }
     48
     49    int numCols = in->numCols;          // Number of columns
     50    int numRows = in->numRows;          // Number of rows
     51
     52    psF32 *input;                       // Input data, for FFTW
     53    if (!in->parent) {
     54        // No parent --- can just use the data buffer
     55        input = psMemIncrRefCounter(in->p_rawDataBuffer);
     56    } else {
     57        // Need to copy the data
     58        input = psAlloc(numCols * numRows);
     59        for (int y = 0; y < numRows; y++) {
     60            memcpy(&input[y * numRows], in->data.F32[y], numCols);
     61        }
     62    }
     63
     64    // Do the FFT
     65    fftwf_complex *out = fftwf_malloc((numCols/2 + 1) * numRows * sizeof(fftwf_complex)); // Output data
     66    fftwf_plan plan = fftwf_plan_dft_r2c_2d(numRows, numCols, input, out, FFTW_PLAN_RIGOR);
     67    fftwf_execute(plan);
     68    fftwf_destroy_plan(plan);
     69    psFree(input);
     70
     71    // Pull the real and imaginary parts out
     72    numCols = numCols/2 + 1;            // x dimension is halved by FFTW
     73    *real = psImageRecycle(*real, numCols, numRows, PS_TYPE_F32);
     74    *imag = psImageRecycle(*imag, numCols, numRows, PS_TYPE_F32);
     75    for (int y = 0, index = 0; y < numRows; y++) {
     76        for (int x = 0; x < numCols; index++, x++) {
     77#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
     78            // C99 complex support
     79            (*real)->data.F32[y][x] = creal(out[index]);
     80            (*imag)->data.F32[y][x] = cimag(out[index]);
     81#else
     82            // FFTW's backup complex support
     83            (*real)->data.F32[y][x] = out[index][0];
     84            (*imag)->data.F32[y][x] = out[index][1];
     85#endif
     86        }
     87    }
     88    fftwf_free(out);
     89
     90    return true;
     91}
     92
     93bool psImageBackwardFFT(psImage **out, const psImage *real, const psImage *imag, int origCols)
     94{
     95    PS_ASSERT_IMAGE_NON_NULL(real, false);
     96    PS_ASSERT_IMAGE_TYPE(real, PS_TYPE_F32, false);
     97    PS_ASSERT_IMAGE_NON_NULL(imag, false);
     98    PS_ASSERT_IMAGE_TYPE(imag, PS_TYPE_F32, false);
     99    PS_ASSERT_IMAGES_SIZE_EQUAL(real, imag, false);
     100    PS_ASSERT_PTR_NON_NULL(out, false);
     101
     102    int numCols = real->numCols;        // Number of columns
     103    int numRows = real->numRows;        // Number of rows
     104
     105    // Because of the way FFT r2c and c2r work, need the number of columns in the target to be:
     106    // 2 * numCols = 2*(numCols/2 + 1) = origCols % 2 ? origCols + 1 : origCols + 2
     107    if (2 * numCols != ((origCols % 2) ? origCols + 1 : origCols + 2)) {
     108        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
     109                "Number of columns in FFT-ed images (%d) should be origCols (%d) / 2 + 1 = %d",
     110                numCols, origCols, origCols/2 + 1);
     111        return false;
     112    }
     113
     114    if (*out && (*out)->parent) {
     115        // It has a parent, so we can't write directly into the buffer.
     116        // It had better be the correct size and type, because we don't want to resize a child image
     117        PS_ASSERT_IMAGE_SIZE(*out, origCols, numRows, false);
     118        PS_ASSERT_IMAGE_TYPE(*out, PS_TYPE_F32, false);
     119    }
     120
     121
     122    // Make sure the system-level wisdom information is imported.
     123    if (!fftwWisdomImported) {
     124        fftwf_import_system_wisdom();
     125        fftwWisdomImported = true;
     126    }
     127
     128    // Stuff the real and imaginary parts in
     129    psF32 *target = psAlloc(2 * numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Target for FFTW
     130    fftwf_complex *in = fftwf_malloc(numCols * numRows * sizeof(fftwf_complex)); // Input data
     131    for (int y = 0, index = 0; y < numRows; y++) {
     132        for (int x = 0; x < numCols; index++, x++) {
     133#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
     134            // C99 complex support
     135            in[index] = real->data.F32[y][x] + imag->data.F32[y][x] * I;
     136#else
     137            // FFTW's backup complex support
     138            in[index][0] = real->data.F32[y][x];
     139            in[index][1] = imag->data.F32[y][x];
     140#endif
     141        }
     142    }
     143
     144    // Do the FFT
     145    fftwf_plan plan = fftwf_plan_dft_c2r_2d(numRows, origCols, in, target, FFTW_PLAN_RIGOR);
     146    fftwf_execute(plan);
     147    fftwf_destroy_plan(plan);
     148    fftwf_free(in);
     149
     150    if (!(*out) || !(*out)->parent) {
     151        *out = psImageRecycle(*out, origCols, numRows, PS_TYPE_F32);
     152    }
     153
     154    // Copy the target pixels into the output
     155    // There's a slight offset --- target has the additional padding required by FFTW,
     156    // while the output doesn't.
     157    for (int y = 0, index = 0; y < numRows; y++, index += origCols) {
     158        memcpy(&(*out)->data.F32[y][0], &target[index], numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
     159    }
     160    psFree(target);
     161
     162    return true;
     163}
     164
     165psImage* psImagePowerSpectrum(psImage *out, const psImage *in)
     166{
     167    PS_ASSERT_IMAGE_NON_NULL(in, NULL);
     168    PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, NULL);
     169
     170    psImage *real = NULL;              // Real component of FFT
     171    psImage *imag = NULL;              // Imaginary component of FFT
     172
     173    if (!psImageForwardFFT(&real, &imag, in)) {
     174        psError(PS_ERR_UNKNOWN, false, "Unable to perform forward FFT.");
    44175        return NULL;
    45176    }
    46177
    47     if ( ((direction & PS_FFT_FORWARD) != 0) ) {
    48         if ((direction & PS_FFT_REVERSE) != 0) {
    49             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    50                     _("Can not specify both PS_FFT_FORWARD and PS_FFT_REVERSE options."));
    51             psFree(out);
    52             return NULL;
    53         }
    54         if ((direction & PS_FFT_REAL_RESULT) != 0) {
    55             psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    56                     _("The PS_FFT_FORWARD and PS_FFT_REAL_RESULT combinition is not supported."));
    57             psFree(out);
    58             return NULL;
    59         }
    60     } else if ((direction & PS_FFT_REVERSE) == 0) {
     178    int numCols = real->numCols;        // Number of columns
     179    int numRows = real->numRows;        // Number of rows
     180
     181    float norm = 1.0 / PS_SQR(in->numCols) / PS_SQR(in->numRows);
     182    for (int y = 0; y < numRows; y++) {
     183        for (int x = 0; x < numCols; x++) {
     184            // Power spectrum is the square of the complex modulus
     185            real->data.F32[y][x] = norm * (PS_SQR(real->data.F32[y][x]) + PS_SQR(imag->data.F32[y][x]));
     186        }
     187    }
     188    psFree(imag);
     189
     190    return real;
     191}
     192
     193
     194
     195bool psImageComplexMultiply(psImage **outReal, psImage **outImag,
     196                            const psImage *in1Real, const psImage *in1Imag,
     197                            const psImage *in2Real, const psImage *in2Imag)
     198{
     199    PS_ASSERT_IMAGE_NON_NULL(in1Real, false);
     200    PS_ASSERT_IMAGE_NON_NULL(in1Imag, false);
     201    PS_ASSERT_IMAGE_NON_NULL(in2Real, false);
     202    PS_ASSERT_IMAGE_NON_NULL(in2Imag, false);
     203    PS_ASSERT_IMAGE_TYPE(in1Real, PS_TYPE_F32, false);
     204    PS_ASSERT_IMAGE_TYPE(in1Imag, PS_TYPE_F32, false);
     205    PS_ASSERT_IMAGE_TYPE(in2Real, PS_TYPE_F32, false);
     206    PS_ASSERT_IMAGE_TYPE(in2Imag, PS_TYPE_F32, false);
     207    PS_ASSERT_IMAGES_SIZE_EQUAL(in1Imag, in1Real, false);
     208    PS_ASSERT_IMAGES_SIZE_EQUAL(in2Real, in1Real, false);
     209    PS_ASSERT_IMAGES_SIZE_EQUAL(in2Imag, in1Real, false);
     210    PS_ASSERT_PTR_NON_NULL(outReal, false);
     211    PS_ASSERT_PTR_NON_NULL(outImag, false);
     212    if (*outReal) {
     213        PS_ASSERT_IMAGE_NON_NULL(*outReal, false);
     214        PS_ASSERT_IMAGE_NON_NULL(*outImag, false);
     215    } else if (*outImag) {
    61216        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    62                 _("Must specify either PS_FFT_FORWARD or PS_FFT_REVERSE option."));
    63         psFree(out);
    64         return NULL;
    65     }
    66 
    67     type = image->type.type;
    68 
    69     /* make sure the system-level wisdom information is imported. */
    70     if (!p_fftwWisdomImported) {
    71         fftwf_import_system_wisdom();
    72         p_fftwWisdomImported = true;
    73     }
    74 
    75     numRows = image->numRows;
    76     numCols = image->numCols;
    77 
    78     // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place.
    79     psS32 sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD;
    80 
    81     fftwf_complex* outBuffer = fftwf_malloc(numRows*numCols*sizeof(fftwf_complex));
    82     p_psImageCopyToRawBuffer(outBuffer, image, PS_TYPE_C32);
    83 
    84     plan = fftwf_plan_dft_2d(numRows, numCols,
    85                              outBuffer,
    86                              outBuffer,
    87                              sign,
    88                              PS_FFTW_PLAN_RIGOR);
    89 
    90     /* check if a plan exists now -- if not, it is a real problem at this point */
    91     if (plan == NULL) {
    92         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    93                 _("Could not create a valid FFT plan to perform the transform."));
    94         psFree(out);
    95         return NULL;
    96     }
    97 
    98     /* finally, call FFTW with the plan made above */
    99     fftwf_execute(plan);
    100 
    101     fftwf_destroy_plan(plan);
    102 
    103     if (direction & PS_FFT_REAL_RESULT) {
    104         // n.b., we do this instead of using fftwf_plan_dft_c2r because that
    105         // plan requires a half-image, which would require a image reordering
    106         // that is not as simple as performing a normal complex transform and
    107         // then taking the real part of the result.  If performance here
    108         // becomes an issue, the use of fftwf_plan_dft_r2c should be considered
    109         // as well as fftwf_plan_dft_c2r.
    110         out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
    111         int index = 0;
    112         for (int row=0; row < numRows; row++) {
    113             psF32* outRow = out->data.F32[row];
    114             for (int col=0; col < numCols; col++) {
    115                 outRow[col] = crealf(outBuffer[index++]); // take just the real part
    116             }
    117         }
     217                "If the output real part is provided, the output imaginary part must also be provided.");
     218        return false;
     219    }
     220
     221    int numRows = in1Real->numRows;     // Number of rows
     222    int numCols = in1Real->numCols;     // Number of columns
     223
     224    // Need to worry if the outputs are children
     225    psImage *targetReal, *targetImag;   // Target real and imaginary parts
     226
     227    if ((*outReal)->parent) {
     228        // It had better be the correct size and type, because we don't want to resize a child image
     229        PS_ASSERT_IMAGE_TYPE(*outReal, PS_TYPE_F32, false);
     230        PS_ASSERT_IMAGE_SIZE(*outReal, numCols, numRows, false);
     231        targetReal = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    118232    } else {
    119         out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
    120         int index = 0;
    121         for (int row=0; row < numRows; row++) {
    122             psC32* outRow = out->data.C32[row];
    123             for (int col=0; col < numCols; col++) {
    124                 outRow[col] = outBuffer[index++]; // take just the real part
    125             }
    126         }
    127         //        memcpy(out->rawDataBuffer, outBuffer, numRows*numCols*sizeof(fftwf_complex));
    128     }
    129 
    130     fftwf_free(outBuffer);
    131 
    132     return out;
    133 
    134 }
    135 
    136 psImage* psImageReal(psImage* out, const psImage* in)
    137 {
    138     psElemType type;
    139     psU32 numCols;
    140     psU32 numRows;
    141 
    142     if (in == NULL) {
    143         psFree(out);
    144         return NULL;
    145     }
    146 
    147     type = in->type.type;
    148     numCols = in->numCols;
    149     numRows = in->numRows;
    150 
    151     /* if not a complex number, this is logically just a copy then */
    152     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    153         return psImageCopy(out, in, type);
    154     }
    155 
    156     if (type == PS_TYPE_C32) {
    157         psF32* outRow;
    158         psC32* inRow;
    159 
    160         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
    161         for (psU32 row = 0; row < numRows; row++) {
    162             outRow = out->data.F32[row];
    163             inRow = in->data.C32[row];
    164 
    165             for (psU32 col = 0; col < numCols; col++) {
    166                 outRow[col] = crealf(inRow[col]);
    167             }
    168         }
    169     } else if (type == PS_TYPE_C64) {
    170         psF64* outRow;
    171         psC64* inRow;
    172 
    173         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
    174         for (psU32 row = 0; row < numRows; row++) {
    175             outRow = out->data.F64[row];
    176             inRow = in->data.C64[row];
    177 
    178             for (psU32 col = 0; col < numCols; col++) {
    179                 outRow[col] = creal(inRow[col]);
    180             }
    181         }
     233        *outReal = psImageRecycle(*outReal, numCols, numRows, PS_TYPE_F32);
     234        targetReal = psMemIncrRefCounter(*outReal);
     235    }
     236    if ((*outImag)->parent) {
     237        // It had better be the correct size and type, because we don't want to resize a child image
     238        if ((*outReal)->numCols != numCols || (*outReal)->numRows != numRows ||
     239            (*outImag)->type.type != PS_TYPE_F32) {
     240            // Plug potential memory leak --- need to free targetReal if there's a problem
     241            psFree(targetReal);
     242        }
     243        PS_ASSERT_IMAGE_TYPE(*outReal, PS_TYPE_F32, false);
     244        PS_ASSERT_IMAGE_SIZE(*outReal, numCols, numRows, false);
     245        targetImag = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    182246    } else {
    183         char* typeStr;
    184         PS_TYPE_NAME(typeStr,type);
    185         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    186                 _("Specified psImage type, %s, is not supported."),
    187                 typeStr);
    188 
    189         psFree(out);
    190         return NULL;
    191     }
    192 
    193     return out;
    194 }
    195 
    196 psImage* psImageImaginary(psImage* out, const psImage* in)
    197 {
    198     psElemType type;
    199     psU32 numCols;
    200     psU32 numRows;
    201 
    202     if (in == NULL) {
    203         psFree(out);
    204         return NULL;
    205     }
    206 
    207     type = in->type.type;
    208     numCols = in->numCols;
    209     numRows = in->numRows;
    210 
    211     /* if not a complex image type, this is logically just zeroed image of same size */
    212     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    213         out = psImageRecycle(out, numCols, numRows, type);
    214         memset(out->data.V[0], 0, PSELEMTYPE_SIZEOF(type) * numCols * numRows);
    215         return out;
    216     }
    217 
    218     if (type == PS_TYPE_C32) {
    219         psF32* outRow;
    220         psC32* inRow;
    221 
    222         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
    223         for (psU32 row = 0; row < numRows; row++) {
    224             outRow = out->data.F32[row];
    225             inRow = in->data.C32[row];
    226 
    227             for (psU32 col = 0; col < numCols; col++) {
    228                 outRow[col] = cimagf(inRow[col]);
    229             }
    230         }
    231     } else if (type == PS_TYPE_C64) {
    232         psF64* outRow;
    233         psC64* inRow;
    234 
    235         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
    236         for (psU32 row = 0; row < numRows; row++) {
    237             outRow = out->data.F64[row];
    238             inRow = in->data.C64[row];
    239 
    240             for (psU32 col = 0; col < numCols; col++) {
    241                 outRow[col] = cimag(inRow[col]);
    242             }
    243         }
    244     } else {
    245         char* typeStr;
    246         PS_TYPE_NAME(typeStr,type);
    247         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    248                 _("Specified psImage type, %s, is not supported."),
    249                 typeStr);
    250         psFree(out);
    251         return NULL;
    252     }
    253 
    254     return out;
    255 }
    256 
    257 psImage* psImageComplex(psImage* out, const psImage* real, const psImage* imag)
    258 {
    259     psElemType type;
    260     psU32 numCols;
    261     psU32 numRows;
    262 
    263     if (real == NULL || imag == NULL) {
    264         psFree(out);
    265         return NULL;
    266     }
    267 
    268     type = real->type.type;
    269     numCols = real->numCols;
    270     numRows = real->numRows;
    271 
    272     if (imag->type.type != type) {
    273         char* typeStrReal;
    274         char* typeStrImag;
    275         PS_TYPE_NAME(typeStrReal,type);
    276         PS_TYPE_NAME(typeStrImag,imag->type.type);
    277         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    278                 _("Real psImage type (%s) and imaginary psImage type (%s) must be the same."),
    279                 typeStrReal,typeStrImag);
    280         psFree(out);
    281         return NULL;
    282     }
    283 
    284     if (imag->numCols != numCols || imag->numRows != numRows) {
    285         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    286                 _("Real psImage size (%dx%d) and imaginary psImage size (%dx%d) must be the same."),
    287                 numCols, numRows, imag->numCols, imag->numRows);
    288         psFree(out);
    289         return NULL;
    290     }
    291 
    292     if (type == PS_TYPE_F32) {
    293         psC32* outRow;
    294         psF32* realRow;
    295         psF32* imagRow;
    296 
    297         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32);
    298 
    299         for (psU32 row = 0; row < numRows; row++) {
    300             outRow = out->data.C32[row];
    301             realRow = real->data.F32[row];
    302             imagRow = imag->data.F32[row];
    303 
    304             for (psU32 col = 0; col < numCols; col++) {
    305                 outRow[col] = realRow[col] + I * imagRow[col];
    306             }
    307         }
    308     } else if (type == PS_TYPE_F64) {
    309         psC64* outRow;
    310         psF64* realRow;
    311         psF64* imagRow;
    312 
    313         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64);
    314         for (psU32 row = 0; row < numRows; row++) {
    315             outRow = out->data.C64[row];
    316             realRow = real->data.F64[row];
    317             imagRow = imag->data.F64[row];
    318 
    319             for (psU32 col = 0; col < numCols; col++) {
    320                 outRow[col] = realRow[col] + I * imagRow[col];
    321             }
    322         }
    323     } else {
    324         char* typeStr;
    325         PS_TYPE_NAME(typeStr,type);
    326         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    327                 _("Input psImage type, %s, is required to be either psF32 or psF64."),
    328                 typeStr);
    329         psFree(out);
    330         return NULL;
    331     }
    332 
    333 
    334     return out;
    335 }
    336 
    337 psImage* psImageConjugate(psImage* out, const psImage* in)
    338 {
    339     psElemType type;
    340     psU32 numCols;
    341     psU32 numRows;
    342 
    343     if (in == NULL) {
    344         psFree(out);
    345         return NULL;
    346     }
    347 
    348     type = in->type.type;
    349     numCols = in->numCols;
    350     numRows = in->numRows;
    351 
    352     /* if not a complex image, this is logically just a image copy */
    353     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    354         return psImageCopy(out, in, type);
    355     }
    356 
    357     if (type == PS_TYPE_C32) {
    358         psC32* outRow;
    359         psC32* inRow;
    360 
    361         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32);
    362         for (psU32 row = 0; row < numRows; row++) {
    363             outRow = out->data.C32[row];
    364             inRow = in->data.C32[row];
    365 
    366             for (psU32 col = 0; col < numCols; col++) {
    367                 outRow[col] = crealf(inRow[col]) - I * cimagf(inRow[col]);
    368             }
    369         }
    370     } else if (type == PS_TYPE_C64) {
    371         psC64* outRow;
    372         psC64* inRow;
    373 
    374         out = psImageRecycle(out, numCols, numRows, PS_TYPE_C64);
    375         for (psU32 row = 0; row < numRows; row++) {
    376             outRow = out->data.C64[row];
    377             inRow = in->data.C64[row];
    378 
    379             for (psU32 col = 0; col < numCols; col++) {
    380                 outRow[col] = creal(inRow[col]) - I * cimag(inRow[col]);
    381             }
    382         }
    383     } else {
    384         char* typeStr;
    385         PS_TYPE_NAME(typeStr,type);
    386         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    387                 _("Input psImage type, %s, is required to be either psC32 or psC64."),
    388                 typeStr);
    389         psFree(out);
    390         return NULL;
    391     }
    392 
    393     return out;
    394 }
    395 
    396 psImage* psImagePowerSpectrum(psImage* out, const psImage* in)
    397 {
    398     psElemType type;
    399     psU32 numCols;
    400     psU32 numRows;
    401 
    402     if (in == NULL) {
    403         psFree(out);
    404         return NULL;
    405     }
    406 
    407     type = in->type.type;
    408     numCols = in->numCols;
    409     numRows = in->numRows;
    410 
    411     if (type == PS_TYPE_C32) {
    412         psF32* outRow;
    413         psC32* inRow;
    414         psF32 real;
    415         psF32 imag;
    416         psF32 fNumCols = numCols;
    417         psF32 fNumRows = numRows;
    418         psF32 numElementsSquared = fNumCols * fNumCols * fNumRows * fNumRows;
    419 
    420         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
    421         for (psU32 row = 0; row < numRows; row++) {
    422             outRow = out->data.F32[row];
    423             inRow = in->data.C32[row];
    424 
    425             for (psU32 col = 0; col < numCols; col++) {
    426                 real = crealf(inRow[col]);
    427                 imag = cimagf(inRow[col]);
    428                 outRow[col] = (real * real + imag * imag) / numElementsSquared;
    429             }
    430         }
    431     } else if (type == PS_TYPE_C64) {
    432         psF64* outRow;
    433         psC64* inRow;
    434         psF64 real;
    435         psF64 imag;
    436         psF64 numElementsSquared = numCols * numCols * numRows * numRows;
    437 
    438         out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64);
    439         for (psU32 row = 0; row < numRows; row++) {
    440             outRow = out->data.F64[row];
    441             inRow = in->data.C64[row];
    442 
    443             for (psU32 col = 0; col < numCols; col++) {
    444                 real = creal(inRow[col]);
    445                 imag = cimag(inRow[col]);
    446                 outRow[col] = (real * real + imag * imag) / numElementsSquared;
    447             }
    448         }
    449     } else {
    450         char* typeStr;
    451         PS_TYPE_NAME(typeStr,type);
    452         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    453                 _("Input psImage type, %s, is required to be either psC32 or psC64."),
    454                 typeStr);
    455         psFree(out);
    456         return NULL;
    457     }
    458 
    459     return out;
    460 
    461 }
     247        *outImag = psImageRecycle(*outImag, numCols, numRows, PS_TYPE_F32);
     248        targetImag = psMemIncrRefCounter(*outImag);
     249    }
     250
     251    for (int y = 0; y < numRows; y++) {
     252        for (int x = 0; x < numCols; x++) {
     253            // (a + bi) * (c + di) = (ac - bd) + (bc + ad)i
     254            float real = in1Real->data.F32[y][x] * in2Real->data.F32[y][x] -
     255                in1Imag->data.F32[y][x] * in2Imag->data.F32[y][x];
     256            float imag = in1Imag->data.F32[y][x] * in2Real->data.F32[y][x] +
     257                in1Real->data.F32[y][x] * in2Imag->data.F32[y][x];
     258            targetReal->data.F32[y][x] = real;
     259            targetImag->data.F32[y][x] = imag;
     260        }
     261    }
     262
     263    if ((*outReal)->parent) {
     264        *outReal = psImageCopy(*outReal, targetReal, PS_TYPE_F32);
     265    }
     266    if ((*outImag)->parent) {
     267        *outImag = psImageCopy(*outImag, targetImag, PS_TYPE_F32);
     268    }
     269
     270    psFree(targetReal);
     271    psFree(targetImag);
     272
     273    return true;
     274}

  • trunk/psLib/src/fft/psImageFFT.h

    r11248 r11703  
    1 /* @file  psImageFFT.h
    2  * @brief Contains FFT transform related functions for psImage
    3  *
    4  * @author Robert DeSonia, MHPCC
    5  *
    6  * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2007-01-23 22:47:22 $
    8  * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    9  */
     1/// @file  psImageFFT.h
     2/// @brief Contains FFT transform related functions for psImage
     3///
     4/// @author Paul Price, IfA
     5/// @author Robert DeSonia, MHPCC
     6///
     7/// @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     8/// @date $Date: 2007-02-08 04:17:58 $
     9/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     10///
    1011
    1112#ifndef PS_IMAGE_FFT_H
    … …  
    1314
    1415#include "psImage.h"
    15 #include "psVectorFFT.h"               // for psFFTFlags
    1616
    1717/// @addtogroup MathOps Mathematical Operations
    1818/// @{
    1919
    20 /** Forward and reverse FFT calculations.
    21  *
    22  *  This takes as input the image of interest (in) and the direction
    23  *  (direction), which is specified by an enumerated type psFftDirection.
    24  *  The input image may be of type psF32 or psC32, the result is always
    25  *  psC32. If the input vector is psF32, the direction must be forward.
    26  *
    27  *  @return psImage* the FFT transformation result
    28  */
    29 psImage* psImageFFT(
    30     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    31     const psImage* image,              ///< the psImage to apply transform to
    32     psFFTFlags direction               ///< the direction of the transform
    33 );
     20/// Forward FFT of an image
     21///
     22/// Applies a forward FFT (exponent -1), with the result returned in both real and imaginary parts.  The FFT
     23/// is not normalised (a forward followed by a reverse is the original scaled by the total number of pixels).
     24/// The FFT takes advantage of the fact that the input is purely real; hence the output number of columns is
     25/// numCols/2 + 1 (with division rounding down).  Only implemented for F32 input.
     26bool psImageForwardFFT(psImage **real,  ///< Real part of FFT
     27                       psImage **imag,  ///< Imaginary part of FFT
     28                       const psImage *in///< Input image (F32)
     29    );
    3430
    35 /** extract the real portion of a complex image
    36  *
    37  *  @return psImage*   real portion of the input image.
    38  */
    39 psImage* psImageReal(
    40     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    41     const psImage* in                  ///< the psImage to extract real portion from
    42 );
     31/// Backward FFT of an image
     32///
     33/// Applies a backward FFT (exponent +1) from the real and imaginary parts with the (purely) real result
     34/// returned.  The FFT is not normalised (a forward followed by a reverse is the original scaled by the total
     35/// number of pixels).  The FFT takes advantage of the fact that the output will be purely real; hence the
     36/// input number of columns is numCols/2 + 1 (with division rounding down); to manage the redundancy (is the
     37/// original number of columns even or odd?), we need the original number of columns (the number of columns of
     38/// the image that was input to psImageForwardFFT) to be provided.  Only implemented for F32 input.
     39bool psImageBackwardFFT(psImage **out,///< Output image
     40                        const psImage *real, ///< Real input (F32)
     41                        const psImage *imag, ///< Imaginary input (F32)
     42                        int origCols    ///< Original number of columns
     43    );
    4344
    44 /** extract the imaginary portion of a complex image
    45  *
    46  *  @return psImage*   imaginary portion of the input image.
    47  */
    48 psImage* psImageImaginary(
    49     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    50     const psImage* in                  ///< the psImage to extract imaginary portion from
    51 );
    5245
    53 /** creates a complex image from separate real and imaginary plane images
    54  *
    55  *  @return psImage*   resulting complex image
    56  */
    57 psImage* psImageComplex(
    58     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    59     const psImage* real,               ///< the real plane image
    60     const psImage* imag                ///< the imaginary plane image
    61 );
     46/// Power spectrum of an image
     47///
     48/// Generates the power spectrum of an image.  Only implemented for F32 input.
     49psImage* psImagePowerSpectrum(psImage *out, const psImage *in);
    6250
    63 /** computes the complex conjugate of an image
    64  *
    65  *  @return psImage*   the complex conjugate of the 'in' image
    66  */
    67 psImage* psImageConjugate(
    68     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    69     const psImage* in                  ///< the psImage to compute conjugate of
    70 );
    71 
    72 /** computes the power spectrum of an image
    73  *
    74  *  @return psImage*   the power spectrum of the 'in' image
    75  */
    76 psImage* psImagePowerSpectrum(
    77     psImage* out,                      ///< a psImage to recycle.  If NULL, a new psImage is made.
    78     const psImage* in                  ///< the psImage to power spectrum of
    79 );
     51/// Multiply complex images
     52///
     53/// The input images are the real and imaginary parts of each of two images.  The real and imaginary parts
     54/// of the output image are returned.  Only implemented for F32 input.
     55bool psImageComplexMultiply(psImage **outReal, ///< Real part of output
     56                            psimage **outImag, ///< Imaginary part of output
     57                            const psImage *in1Real, ///< Real part of input 1
     58                            const psImage *in1Imag, ///< Imaginary part of input 1
     59                            const psImage *in2Real, ///< Real part of input 2
     60                            const psImage *in2Imag ///< Imaginary part of input 2
     61    );
    8062
    8163/// @}

  • 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 }

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

    r11248 r11703  
    1 /* @file  psVectorFFT.h
    2  * @brief Contains FFT transform related functions for psVector
    3  *
    4  * @author Robert DeSonia, MHPCC
    5  *
    6  * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
    7  * @date $Date: 2007-01-23 22:47:22 $
    8  * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    9  */
     1/// @file  psVectorFFT.h
     2/// @brief Contains FFT transform related functions for psVector
     3///
     4/// @author Paul Price, IfA
     5/// @author Robert DeSonia, MHPCC
     6///
     7/// @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
     8/// @date $Date: 2007-02-08 04:17:58 $
     9/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
     10///
    1011
    1112#ifndef PS_VECTOR_FFT_H
    … …  
    1718/// @{
    1819
    19 /** Specify direction of FFT */
    20 typedef enum {
    21     /// psImageFFT/psVectorFFT should perform a forward FFT.
    22     PS_FFT_FORWARD = 1,
     20/// Forward FFT of a vector
     21///
     22/// Applies a forward FFT (exponent -1), with the result returned in both real and imaginary parts.  The FFT
     23/// is not normalised (a forward followed by a reverse is the original scaled by the number of elements).  The
     24/// FFT takes advantage of the fact that the input is purely real; hence the output size is N/2 + 1 (with
     25/// division rounding down).  Only implemented for F32 input.
     26bool psVectorForwardFFT(psVector **real,///< Real part of FFT
     27                        psVector **imag,///< Imaginary part of FFT
     28                        const psVector *in ///< Input vector (F32)
     29    );
    2330
    24     /// psImageFFT/psVectorFFT should perform a reverse FFT.
    25     PS_FFT_REVERSE = 2,
     31/// Backward FFT of a vector
     32///
     33/// Applies a backward FFT (exponent +1) from the real and imaginary parts with the (purely) real result
     34/// returned.  The FFT is not normalised (a forward followed by a reverse is the original scaled by the number
     35/// of elements).  The FFT takes advantage of the fact that the output will be purely real; hence the input
     36/// size is N/2 + 1 (with division rounding down); to manage the redundancy (is the original size even or
     37/// odd?), we need the original size (the size of the array that was input to psVectorForwardFFT) to be
     38/// provided.  Only implemented for F32 input.
     39bool psVectorBackwardFFT(psVector **out,///< Output vector
     40                         const psVector *real, ///< Real input (F32)
     41                         const psVector *imag, ///< Imaginary input (F32)
     42                         int origN      ///< Original number of elements
     43    );
    2644
    27     /// psImageFFT/psVectorFFT should perform a reverse FFT with a real result.
    28     PS_FFT_REAL_RESULT = 4
    29 } psFFTFlags;
     45/// Power spectrum of a vector
     46///
     47/// Generates the power spectrum of a vector.  Only implemented for F32 input.
     48psVector *psVectorPowerSpectrum(psVector *out, ///< Output power spectrum, or NULL
     49                                const psVector* in ///< Input vector (F32)
     50    );
    3051
    31 
    32 /** Forward and reverse FFT calculations.
    33  *
    34  *  This takes as input the vector of interest (in) and the direction
    35  *  (direction), which is specified by an enumerated type psFftDirection.
    36  *  The input vector may be of type psF32 or psC32, the result is always
    37  *  psC32. If the input vector is psF32, the direction must be forward.
    38  *
    39  *  @return psVector* the FFT transformation result
    40  */
    41 psVector* psVectorFFT(
    42     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    43     const psVector* in,                ///< the vector to apply transform to
    44     psFFTFlags direction               ///< the direction of the transform
    45 );
    46 
    47 /** extract the real portion of a complex vector
    48  *
    49  *  @return psVector*   real portion of the input vector.
    50  */
    51 psVector* psVectorReal(
    52     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    53     const psVector* in                ///< the psVector to extract real portion from
    54 );
    55 
    56 /** extract the imaginary portion of a complex vector
    57  *
    58  *  @return psVector*   imaginary portion of the input vector.
    59  */
    60 psVector* psVectorImaginary(
    61     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    62     const psVector* in                 ///< the psVector to extract imaginary portion from
    63 );
    64 
    65 /** creates a complex vector from separate real and imaginary vectors
    66  *
    67  *  @return psVector*   resulting complex vector
    68  */
    69 psVector* psVectorComplex(
    70     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    71     const psVector* real,              ///< the real vector
    72     const psVector* imag               ///< the imaginary vector
    73 );
    74 
    75 /** computes the complex conjugate of a vector
    76  *
    77  *  @return psVector*   the complex conjugate of the 'in' vector
    78  */
    79 psVector* psVectorConjugate(
    80     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    81     const psVector* in                 ///< the psVector to compute conjugate of
    82 );
    83 
    84 /** computes the power spectrum of a vector
    85  *
    86  *  @return psVector*   the power spectrum of the 'in' vector
    87  */
    88 psVector* psVectorPowerSpectrum(
    89     psVector* out,                     ///< a psVector to recycle.  If NULL, a new psVector is made.
    90     const psVector* in                 ///< the psVector to power spectrum of
    91 );
     52/// Multiply complex vectors
     53///
     54/// The input vectors are the real and imaginary parts of each of two vectors.  The real and imaginary parts
     55/// of the output vector are returned.  Only implemented for F32 input.
     56bool psVectorComplexMultiply(psVector **outReal, ///< Real part of output
     57                             psVector **outImag, ///< Imaginary part of output
     58                             const psVector *in1Real, ///< Real part of input 1
     59                             const psVector *in1Imag, ///< Imaginary part of input 1
     60                             const psVector *in2Real, ///< Real part of input 2
     61                             const psVector *in2Imag ///< Imaginary part of input 2
     62    );
    9263
    9364/// @}
Note: See TracChangeset for help on using the changeset viewer.