IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 4, 2004, 1:37:39 PM (22 years ago)
Author:
desonia
Message:

found the server astyle upgrade was faulty, so the format was reset.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psFFT.c

    r1366 r1385  
    55*  @author Robert DeSonia, MHPCC
    66*
    7 *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
    8 *  @date $Date: 2004-08-02 19:43:59 $
     7*  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
     8*  @date $Date: 2004-08-04 23:37:39 $
    99*
    1010*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3333    psElemType type;
    3434    fftwf_plan plan;
    35    
     35
    3636    /* got good image data? */
    3737    if ( in == NULL ) {
    38             psFree( out );
    39             return NULL;
    40         }
    41        
    42     type = in->type.type;
    43    
     38        psFree( out );
     39        return NULL;
     40    }
     41
     42    type = in->type.type;
     43
    4444    if ( ( type != PS_TYPE_F32 ) && ( type != PS_TYPE_C32 ) ) {
    45             psError( __func__, "Input image must be a 32-bit float or complex image (type=%d)",
    46                      type );
    47             psFree( out );
    48             return NULL;
    49         }
    50        
     45        psError( __func__, "Input image must be a 32-bit float or complex image (type=%d)",
     46                 type );
     47        psFree( out );
     48        return NULL;
     49    }
     50
    5151    if ( type != PS_TYPE_C32 && direction == PS_FFT_REVERSE ) {
    52             psError( __func__, "Input image must be complex image for reverse FFT (type=%d).",
    53                      type );
    54             psFree( out );
    55             return NULL;
    56            
    57         }
    58        
     52        psError( __func__, "Input image must be complex image for reverse FFT (type=%d).",
     53                 type );
     54        psFree( out );
     55        return NULL;
     56
     57    }
     58
    5959    if ( type != PS_TYPE_F32 && direction == PS_FFT_FORWARD ) {
    60             psError( __func__, "Input image must be real image for forward FFT (type=%d).",
    61                      type );
    62             psFree( out );
    63             return NULL;
    64         }
    65        
     60        psError( __func__, "Input image must be real image for forward FFT (type=%d).",
     61                 type );
     62        psFree( out );
     63        return NULL;
     64    }
     65
    6666    /* make sure the system-level wisdom information is imported. */
    6767    if ( ! p_fftwWisdomImported ) {
    68             fftwf_import_system_wisdom();
    69             p_fftwWisdomImported = true;
    70         }
    71        
     68        fftwf_import_system_wisdom();
     69        p_fftwWisdomImported = true;
     70    }
     71
    7272    numRows = in->numRows;
    7373    numCols = in->numCols;
    74    
     74
    7575    out = psImageCopy( out, in, PS_TYPE_C32 );
    76    
     76
    7777    plan = fftwf_plan_dft_2d( numCols, numRows,
    7878                              ( fftwf_complex* ) out->data.C32[ 0 ],
    … …  
    8080                              direction,
    8181                              P_FFTW_PLAN_RIGOR );
    82                              
     82
    8383    /* check if a plan exists now*/
    8484    if ( plan == NULL ) {
    85             psError( __func__, "Failed to create FFTW plan." );
    86             psFree( out );
    87             return NULL;
    88         }
    89        
     85        psError( __func__, "Failed to create FFTW plan." );
     86        psFree( out );
     87        return NULL;
     88    }
     89
    9090    /* finally, call FFTW with the plan made above */
    9191    fftwf_execute( plan );
    92    
     92
    9393    fftwf_destroy_plan( plan );
    94    
    95     return out;
    96    
     94
     95    return out;
     96
    9797}
    9898
    … …  
    103103    unsigned int numCols;
    104104    unsigned int numRows;
    105    
    106    
    107     if ( in == NULL ) {
    108             psFree( out );
    109             return NULL;
    110         }
    111        
     105
     106
     107    if ( in == NULL ) {
     108        psFree( out );
     109        return NULL;
     110    }
     111
    112112    type = in->type.type;
    113113    numCols = in->numCols;
    114114    numRows = in->numRows;
    115    
     115
    116116    /* if not a complex number, this is logically just a copy */
    117117    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    118             // Warn user, as this is probably not expected
    119             psLogMsg( __func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
    120                       "Just an image copy was performed." );
    121             return psImageCopy( out, in, type );
    122         }
    123        
     118        // Warn user, as this is probably not expected
     119        psLogMsg( __func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
     120                  "Just an image copy was performed." );
     121        return psImageCopy( out, in, type );
     122    }
     123
    124124    if ( type == PS_TYPE_C32 ) {
    125             psF32 * outRow;
    126             psC32* inRow;
    127            
    128             out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
    129             for ( unsigned int row = 0;row < numRows;row++ ) {
    130                     outRow = out->data.F32[ row ];
    131                     inRow = in->data.C32[ row ];
    132                    
    133                     for ( unsigned int col = 0;col < numCols;col++ ) {
    134                             outRow[ col ] = crealf( inRow[ col ] );
    135                         }
    136                 }
    137         } else if ( type == PS_TYPE_C64 ) {
     125        psF32 * outRow;
     126        psC32* inRow;
     127
     128        out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
     129        for ( unsigned int row = 0;row < numRows;row++ ) {
     130            outRow = out->data.F32[ row ];
     131            inRow = in->data.C32[ row ];
     132
     133            for ( unsigned int col = 0;col < numCols;col++ ) {
     134                outRow[ col ] = crealf( inRow[ col ] );
     135            }
     136        }
     137    } else
     138        if ( type == PS_TYPE_C64 ) {
    138139            psF64 * outRow;
    139140            psC64* inRow;
    140            
     141
    141142            out = psImageRecycle( out, numCols, numRows, PS_TYPE_F64 );
    142143            for ( unsigned int row = 0;row < numRows;row++ ) {
    143                     outRow = out->data.F64[ row ];
    144                     inRow = in->data.C64[ row ];
    145                    
    146                     for ( unsigned int col = 0;col < numCols;col++ ) {
    147                             outRow[ col ] = creal( inRow[ col ] );
    148                         }
     144                outRow = out->data.F64[ row ];
     145                inRow = in->data.C64[ row ];
     146
     147                for ( unsigned int col = 0;col < numCols;col++ ) {
     148                    outRow[ col ] = creal( inRow[ col ] );
    149149                }
     150            }
    150151        } else {
    151152            psError( __func__, "Can not extract real component from given image type (%d).",
    … …  
    154155            return NULL;
    155156        }
    156        
     157
    157158    return out;
    158159}
    … …  
    163164    unsigned int numCols;
    164165    unsigned int numRows;
    165    
    166    
    167     if ( in == NULL ) {
    168             psFree( out );
    169             return NULL;
    170         }
    171        
     166
     167
     168    if ( in == NULL ) {
     169        psFree( out );
     170        return NULL;
     171    }
     172
    172173    type = in->type.type;
    173174    numCols = in->numCols;
    174175    numRows = in->numRows;
    175    
     176
    176177    /* if not a complex number, this is logically just zeroed image of same size */
    177178    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    178             // Warn user, as this is probably not expected
    179             psLogMsg( __func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
    180                       "A zero image was returned." );
    181             out = psImageRecycle( out, numCols, numRows, type );
    182             memset( out->data.V[ 0 ], 0, PSELEMTYPE_SIZEOF( type ) * numCols * numRows );
    183             return out;
    184         }
    185        
     179        // Warn user, as this is probably not expected
     180        psLogMsg( __func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
     181                  "A zero image was returned." );
     182        out = psImageRecycle( out, numCols, numRows, type );
     183        memset( out->data.V[ 0 ], 0, PSELEMTYPE_SIZEOF( type ) * numCols * numRows );
     184        return out;
     185    }
     186
    186187    if ( type == PS_TYPE_C32 ) {
    187             psF32 * outRow;
    188             psC32* inRow;
    189            
    190             out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
    191             for ( unsigned int row = 0;row < numRows;row++ ) {
    192                     outRow = out->data.F32[ row ];
    193                     inRow = in->data.C32[ row ];
    194                    
    195                     for ( unsigned int col = 0;col < numCols;col++ ) {
    196                             outRow[ col ] = cimagf( inRow[ col ] );
    197                         }
    198                 }
    199         } else if ( type == PS_TYPE_C64 ) {
     188        psF32 * outRow;
     189        psC32* inRow;
     190
     191        out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
     192        for ( unsigned int row = 0;row < numRows;row++ ) {
     193            outRow = out->data.F32[ row ];
     194            inRow = in->data.C32[ row ];
     195
     196            for ( unsigned int col = 0;col < numCols;col++ ) {
     197                outRow[ col ] = cimagf( inRow[ col ] );
     198            }
     199        }
     200    } else
     201        if ( type == PS_TYPE_C64 ) {
    200202            psF64 * outRow;
    201203            psC64* inRow;
    202            
     204
    203205            out = psImageRecycle( out, numCols, numRows, PS_TYPE_F64 );
    204206            for ( unsigned int row = 0;row < numRows;row++ ) {
    205                     outRow = out->data.F64[ row ];
    206                     inRow = in->data.C64[ row ];
    207                    
    208                     for ( unsigned int col = 0;col < numCols;col++ ) {
    209                             outRow[ col ] = cimag( inRow[ col ] );
    210                         }
     207                outRow = out->data.F64[ row ];
     208                inRow = in->data.C64[ row ];
     209
     210                for ( unsigned int col = 0;col < numCols;col++ ) {
     211                    outRow[ col ] = cimag( inRow[ col ] );
    211212                }
     213            }
    212214        } else {
    213215            psError( __func__, "Can not extract imaginary component from given image type (%d).",
    … …  
    216218            return NULL;
    217219        }
    218        
     220
    219221    return out;
    220222}
    … …  
    225227    unsigned int numCols;
    226228    unsigned int numRows;
    227    
    228    
     229
     230
    229231    if ( real == NULL || imag == NULL ) {
    230             psFree( out );
    231             return NULL;
    232         }
    233        
     232        psFree( out );
     233        return NULL;
     234    }
     235
    234236    type = real->type.type;
    235237    numCols = real->numCols;
    236238    numRows = real->numRows;
    237    
     239
    238240    if ( imag->type.type != type ) {
    239             psError( __func__, "The inputs to psImageComplex must be the same type." );
    240             psFree( out );
    241             return NULL;
    242         }
    243        
     241        psError( __func__, "The inputs to psImageComplex must be the same type." );
     242        psFree( out );
     243        return NULL;
     244    }
     245
    244246    if ( imag->numCols != numCols ||
    245247            imag->numRows != numRows ) {
    246             psError( __func__, "The inputs to psImageComplex must be the same dimensions." );
    247             psFree( out );
    248             return NULL;
    249         }
    250        
     248        psError( __func__, "The inputs to psImageComplex must be the same dimensions." );
     249        psFree( out );
     250        return NULL;
     251    }
     252
    251253    if ( PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    252             psError( __func__, "The inputs to psImageComplex can not be complex." );
    253             psFree( out );
    254             return NULL;
    255         }
    256        
     254        psError( __func__, "The inputs to psImageComplex can not be complex." );
     255        psFree( out );
     256        return NULL;
     257    }
     258
    257259    if ( type != PS_TYPE_F32 && type != PS_TYPE_F64 ) {
    258             psError( __func__, "The input type to psImageComplex must be a floating point." );
    259             psFree( out );
    260             return NULL;
    261         }
    262        
     260        psError( __func__, "The input type to psImageComplex must be a floating point." );
     261        psFree( out );
     262        return NULL;
     263    }
     264
    263265    if ( type == PS_TYPE_F32 ) {
    264             psC32 * outRow;
    265             psF32* realRow;
    266             psF32* imagRow;
    267            
    268             out = psImageRecycle( out, numCols, numRows, PS_TYPE_C32 );
    269            
    270             for ( unsigned int row = 0;row < numRows;row++ ) {
    271                     outRow = out->data.C32[ row ];
    272                     realRow = real->data.F32[ row ];
    273                     imagRow = imag->data.F32[ row ];
    274                    
    275                     for ( unsigned int col = 0;col < numCols;col++ ) {
    276                             outRow[ col ] = realRow[ col ] + I * imagRow[ col ];
    277                         }
    278                 }
    279         } else if ( type == PS_TYPE_F64 ) {
     266        psC32 * outRow;
     267        psF32* realRow;
     268        psF32* imagRow;
     269
     270        out = psImageRecycle( out, numCols, numRows, PS_TYPE_C32 );
     271
     272        for ( unsigned int row = 0;row < numRows;row++ ) {
     273            outRow = out->data.C32[ row ];
     274            realRow = real->data.F32[ row ];
     275            imagRow = imag->data.F32[ row ];
     276
     277            for ( unsigned int col = 0;col < numCols;col++ ) {
     278                outRow[ col ] = realRow[ col ] + I * imagRow[ col ];
     279            }
     280        }
     281    } else
     282        if ( type == PS_TYPE_F64 ) {
    280283            psC64 * outRow;
    281284            psF64* realRow;
    282285            psF64* imagRow;
    283            
     286
    284287            out = psImageRecycle( out, numCols, numRows, PS_TYPE_C64 );
    285288            for ( unsigned int row = 0;row < numRows;row++ ) {
    286                     outRow = out->data.C64[ row ];
    287                     realRow = real->data.F64[ row ];
    288                     imagRow = imag->data.F64[ row ];
    289                    
    290                     for ( unsigned int col = 0;col < numCols;col++ ) {
    291                             outRow[ col ] = realRow[ col ] + I * imagRow[ col ];
    292                         }
     289                outRow = out->data.C64[ row ];
     290                realRow = real->data.F64[ row ];
     291                imagRow = imag->data.F64[ row ];
     292
     293                for ( unsigned int col = 0;col < numCols;col++ ) {
     294                    outRow[ col ] = realRow[ col ] + I * imagRow[ col ];
    293295                }
     296            }
    294297        } else {
    295298            psError( __func__, "Can not merge real and imaginary portions for given image type (%d).",
    … …  
    298301            return NULL;
    299302        }
    300        
     303
    301304    return out;
    302305}
    … …  
    307310    unsigned int numCols;
    308311    unsigned int numRows;
    309    
    310    
    311     if ( in == NULL ) {
    312             psFree( out );
    313             return NULL;
    314         }
    315        
     312
     313
     314    if ( in == NULL ) {
     315        psFree( out );
     316        return NULL;
     317    }
     318
    316319    type = in->type.type;
    317320    numCols = in->numCols;
    318321    numRows = in->numRows;
    319    
     322
    320323    /* if not a complex number, this is logically just a image copy */
    321324    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    322             // Warn user, as this is probably not expected
    323             psLogMsg( __func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
    324                       "Image copy was performed instead." );
    325             return psImageCopy( out, in, type );
    326         }
    327        
     325        // Warn user, as this is probably not expected
     326        psLogMsg( __func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
     327                  "Image copy was performed instead." );
     328        return psImageCopy( out, in, type );
     329    }
     330
    328331    if ( type == PS_TYPE_C32 ) {
    329             psC32 * outRow;
    330             psC32* inRow;
    331            
    332             out = psImageRecycle( out, numCols, numRows, PS_TYPE_C32 );
    333             for ( unsigned int row = 0;row < numRows;row++ ) {
    334                     outRow = out->data.C32[ row ];
    335                     inRow = in->data.C32[ row ];
    336                    
    337                     for ( unsigned int col = 0;col < numCols;col++ ) {
    338                             outRow[ col ] = crealf( inRow[ col ] ) - I * cimagf( inRow[ col ] );
    339                         }
    340                 }
    341         } else if ( type == PS_TYPE_C64 ) {
     332        psC32 * outRow;
     333        psC32* inRow;
     334
     335        out = psImageRecycle( out, numCols, numRows, PS_TYPE_C32 );
     336        for ( unsigned int row = 0;row < numRows;row++ ) {
     337            outRow = out->data.C32[ row ];
     338            inRow = in->data.C32[ row ];
     339
     340            for ( unsigned int col = 0;col < numCols;col++ ) {
     341                outRow[ col ] = crealf( inRow[ col ] ) - I * cimagf( inRow[ col ] );
     342            }
     343        }
     344    } else
     345        if ( type == PS_TYPE_C64 ) {
    342346            psC64 * outRow;
    343347            psC64* inRow;
    344            
     348
    345349            out = psImageRecycle( out, numCols, numRows, PS_TYPE_C64 );
    346350            for ( unsigned int row = 0;row < numRows;row++ ) {
    347                     outRow = out->data.C64[ row ];
    348                     inRow = in->data.C64[ row ];
    349                    
    350                     for ( unsigned int col = 0;col < numCols;col++ ) {
    351                             outRow[ col ] = creal( inRow[ col ] ) - I * cimag( inRow[ col ] );
    352                         }
     351                outRow = out->data.C64[ row ];
     352                inRow = in->data.C64[ row ];
     353
     354                for ( unsigned int col = 0;col < numCols;col++ ) {
     355                    outRow[ col ] = creal( inRow[ col ] ) - I * cimag( inRow[ col ] );
    353356                }
     357            }
    354358        } else {
    355359            psError( __func__, "Can not compute complex conjugate for given image type (%d).",
    … …  
    358362            return NULL;
    359363        }
    360        
     364
    361365    return out;
    362366}
    … …  
    368372    unsigned int numRows;
    369373    int numElementsSquared;
    370    
    371     if ( in == NULL ) {
    372             psFree( out );
    373             return NULL;
    374         }
    375        
     374
     375    if ( in == NULL ) {
     376        psFree( out );
     377        return NULL;
     378    }
     379
    376380    type = in->type.type;
    377381    numCols = in->numCols;
    378382    numRows = in->numRows;
    379383    numElementsSquared = numCols * numCols * numRows * numRows;
    380    
     384
    381385    /* if not a complex number, this is not implemented */
    382386    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    383             psError( __func__, "Power Spectrum for non-complex inputs is not implemented." );
    384             psFree( out );
    385             return NULL;
    386         }
    387        
     387        psError( __func__, "Power Spectrum for non-complex inputs is not implemented." );
     388        psFree( out );
     389        return NULL;
     390    }
     391
    388392    if ( type == PS_TYPE_C32 ) {
    389             psF32 * outRow;
    390             psC32* inRow;
    391             psF32 real;
    392             psF32 imag;
    393            
    394            
    395             out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
    396             for ( unsigned int row = 0;row < numRows;row++ ) {
    397                     outRow = out->data.F32[ row ];
    398                     inRow = in->data.C32[ row ];
    399                    
    400                     for ( unsigned int col = 0;col < numCols;col++ ) {
    401                             real = crealf( inRow[ col ] );
    402                             imag = cimagf( inRow[ col ] );
    403                             outRow[ col ] = ( real * real + imag * imag ) / numElementsSquared;
    404                         }
    405                 }
    406         } else if ( type == PS_TYPE_C64 ) {
     393        psF32 * outRow;
     394        psC32* inRow;
     395        psF32 real;
     396        psF32 imag;
     397
     398
     399        out = psImageRecycle( out, numCols, numRows, PS_TYPE_F32 );
     400        for ( unsigned int row = 0;row < numRows;row++ ) {
     401            outRow = out->data.F32[ row ];
     402            inRow = in->data.C32[ row ];
     403
     404            for ( unsigned int col = 0;col < numCols;col++ ) {
     405                real = crealf( inRow[ col ] );
     406                imag = cimagf( inRow[ col ] );
     407                outRow[ col ] = ( real * real + imag * imag ) / numElementsSquared;
     408            }
     409        }
     410    } else
     411        if ( type == PS_TYPE_C64 ) {
    407412            psF64 * outRow;
    408413            psC64* inRow;
    409414            psF64 real;
    410415            psF64 imag;
    411            
    412            
     416
     417
    413418            out = psImageRecycle( out, numCols, numRows, PS_TYPE_F64 );
    414419            for ( unsigned int row = 0;row < numRows;row++ ) {
    415                     outRow = out->data.F64[ row ];
    416                     inRow = in->data.C64[ row ];
    417                    
    418                     for ( unsigned int col = 0;col < numCols;col++ ) {
    419                             real = crealf( inRow[ col ] );
    420                             imag = cimagf( inRow[ col ] );
    421                             outRow[ col ] = real * real + imag * imag / numElementsSquared;
    422                         }
     420                outRow = out->data.F64[ row ];
     421                inRow = in->data.C64[ row ];
     422
     423                for ( unsigned int col = 0;col < numCols;col++ ) {
     424                    real = crealf( inRow[ col ] );
     425                    imag = cimagf( inRow[ col ] );
     426                    outRow[ col ] = real * real + imag * imag / numElementsSquared;
    423427                }
     428            }
    424429        } else {
    425430            psError( __func__, "Can not power spectrum for given image type (%d).",
    … …  
    428433            return NULL;
    429434        }
    430        
    431     return out;
    432    
     435
     436    return out;
     437
    433438}
    434439
    … …  
    440445    psElemType type;
    441446    fftwf_plan plan;
    442    
     447
    443448    /* got good image data? */
    444449    if ( in == NULL ) {
    445             psFree( out );
    446             return NULL;
    447         }
    448        
    449     type = in->type.type;
    450    
     450        psFree( out );
     451        return NULL;
     452    }
     453
     454    type = in->type.type;
     455
    451456    if ( ( type != PS_TYPE_F32 ) && ( type != PS_TYPE_C32 ) ) {
    452             psError( __func__, "Input image must be a 32-bit float or complex image (type=%d)",
    453                      type );
    454             psFree( out );
    455             return NULL;
    456         }
    457        
     457        psError( __func__, "Input image must be a 32-bit float or complex image (type=%d)",
     458                 type );
     459        psFree( out );
     460        return NULL;
     461    }
     462
    458463    if ( ( type != PS_TYPE_C32 ) && ( direction == PS_FFT_REVERSE ) ) {
    459             psError( __func__, "Input image must be complex image for reverse FFT (type=%d).",
    460                      type );
    461             psFree( out );
    462             return NULL;
    463            
    464         }
    465        
     464        psError( __func__, "Input image must be complex image for reverse FFT (type=%d).",
     465                 type );
     466        psFree( out );
     467        return NULL;
     468
     469    }
     470
    466471    if ( ( type != PS_TYPE_F32 ) && ( direction == PS_FFT_FORWARD ) ) {
    467             psError( __func__, "Input image must be real image for forward FFT (type=%d).",
    468                      type );
    469             psFree( out );
    470             return NULL;
    471         }
    472        
     472        psError( __func__, "Input image must be real image for forward FFT (type=%d).",
     473                 type );
     474        psFree( out );
     475        return NULL;
     476    }
     477
    473478    /* make sure the system-level wisdom information is imported. */
    474479    if ( ! p_fftwWisdomImported ) {
    475             fftwf_import_system_wisdom();
    476             p_fftwWisdomImported = true;
    477         }
    478        
     480        fftwf_import_system_wisdom();
     481        p_fftwWisdomImported = true;
     482    }
     483
    479484    numElements = in->n;
    480    
     485
    481486    out = psVectorRecycle( out, PS_TYPE_C32, numElements );
    482487    out->n = numElements;
    483    
     488
    484489    if ( type == PS_TYPE_F32 ) {
    485             // need to convert to complex
    486             psC32 * outVec = out->data.C32;
    487             psF32* inVec = in->data.F32;
    488             for ( unsigned int i = 0;i < numElements;i++ ) {
    489                     outVec[ i ] = inVec[ i ];
    490                 }
    491         } else {
    492             psC32* outVec = out->data.C32;
    493             psC32* inVec = in->data.C32;
    494             for ( unsigned int i = 0;i < numElements;i++ ) {
    495                     outVec[ i ] = inVec[ i ];
    496                 }
    497         }
    498        
     490        // need to convert to complex
     491        psC32 * outVec = out->data.C32;
     492        psF32* inVec = in->data.F32;
     493        for ( unsigned int i = 0;i < numElements;i++ ) {
     494            outVec[ i ] = inVec[ i ];
     495        }
     496    } else {
     497        psC32* outVec = out->data.C32;
     498        psC32* inVec = in->data.C32;
     499        for ( unsigned int i = 0;i < numElements;i++ ) {
     500            outVec[ i ] = inVec[ i ];
     501        }
     502    }
     503
    499504    plan = fftwf_plan_dft_1d( numElements,
    500505                              ( fftwf_complex* ) out->data.C32,
    … …  
    502507                              direction,
    503508                              P_FFTW_PLAN_RIGOR );
    504                              
     509
    505510    /* check if a plan exists now*/
    506511    if ( plan == NULL ) {
    507             psError( __func__, "Failed to create FFTW plan." );
    508             psFree( out );
    509             return NULL;
    510         }
    511        
     512        psError( __func__, "Failed to create FFTW plan." );
     513        psFree( out );
     514        return NULL;
     515    }
     516
    512517    /* finally, call FFTW with the plan made above */
    513518    fftwf_execute( plan );
    514    
     519
    515520    fftwf_destroy_plan( plan );
    516    
     521
    517522    return out;
    518523}
    … …  
    523528    psElemType type;
    524529    unsigned int numElements;
    525    
    526     if ( in == NULL ) {
    527             psFree( out );
    528             return NULL;
    529         }
    530        
     530
     531    if ( in == NULL ) {
     532        psFree( out );
     533        return NULL;
     534    }
     535
    531536    type = in->type.type;
    532537    numElements = in->n;
    533    
     538
    534539    /* if not a complex number, this is logically just a copy */
    535540    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    536             // Warn user, as this is probably not expected
    537             psLogMsg( __func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
    538                       "Just a vector copy was performed." );
    539             out = psVectorRecycle( out, type, numElements );
    540             out->n = numElements;
    541             memcpy( out->data.V, in->data.V, numElements * PSELEMTYPE_SIZEOF( type ) );
    542             return out;
    543         }
    544        
     541        // Warn user, as this is probably not expected
     542        psLogMsg( __func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
     543                  "Just a vector copy was performed." );
     544        out = psVectorRecycle( out, type, numElements );
     545        out->n = numElements;
     546        memcpy( out->data.V, in->data.V, numElements * PSELEMTYPE_SIZEOF( type ) );
     547        return out;
     548    }
     549
    545550    if ( type == PS_TYPE_C32 ) {
    546             psF32 * outVec;
    547             psC32* inVec = in->data.C32;
    548            
    549             out = psVectorRecycle( out, PS_TYPE_F32, numElements );
    550             out->n = numElements;
    551             outVec = out->data.F32;
    552            
    553             for ( unsigned int i = 0;i < numElements;i++ ) {
    554                     outVec[ i ] = crealf( inVec[ i ] );
    555                 }
    556         } else {
    557             psError( __func__, "Can not extract real component from given vector type (%d).",
    558                      type );
    559             psFree( out );
    560             return NULL;
    561         }
    562        
     551        psF32 * outVec;
     552        psC32* inVec = in->data.C32;
     553
     554        out = psVectorRecycle( out, PS_TYPE_F32, numElements );
     555        out->n = numElements;
     556        outVec = out->data.F32;
     557
     558        for ( unsigned int i = 0;i < numElements;i++ ) {
     559            outVec[ i ] = crealf( inVec[ i ] );
     560        }
     561    } else {
     562        psError( __func__, "Can not extract real component from given vector type (%d).",
     563                 type );
     564        psFree( out );
     565        return NULL;
     566    }
     567
    563568    return out;
    564569}
    … …  
    568573    psElemType type;
    569574    unsigned int numElements;
    570    
    571    
    572     if ( in == NULL ) {
    573             psFree( out );
    574             return NULL;
    575         }
    576        
     575
     576
     577    if ( in == NULL ) {
     578        psFree( out );
     579        return NULL;
     580    }
     581
    577582    type = in->type.type;
    578583    numElements = in->n;
    579    
     584
    580585    /* if not a complex number, this is logically just zeroed image of same size */
    581586    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    582             // Warn user, as this is probably not expected
    583             psLogMsg( __func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
    584                       "A zeroed vector was returned." );
    585             out = psVectorRecycle( out, type, numElements );
    586             out->n = numElements;
    587             memset( out->data.V, 0, PSELEMTYPE_SIZEOF( type ) * numElements );
    588             return out;
    589         }
    590        
     587        // Warn user, as this is probably not expected
     588        psLogMsg( __func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
     589                  "A zeroed vector was returned." );
     590        out = psVectorRecycle( out, type, numElements );
     591        out->n = numElements;
     592        memset( out->data.V, 0, PSELEMTYPE_SIZEOF( type ) * numElements );
     593        return out;
     594    }
     595
    591596    if ( type == PS_TYPE_C32 ) {
    592             psF32 * outVec;
    593             psC32* inVec = in->data.C32;
    594            
    595             out = psVectorRecycle( out, PS_TYPE_F32, numElements );
    596             out->n = numElements;
    597             outVec = out->data.F32;
    598            
    599             for ( unsigned int i = 0;i < numElements;i++ ) {
    600                     outVec[ i ] = cimagf( inVec[ i ] );
    601                 }
    602         } else {
    603             psError( __func__, "Can not extract imaginary component from given vector type (%d).",
    604                      type );
    605             psFree( out );
    606             return NULL;
    607         }
    608        
     597        psF32 * outVec;
     598        psC32* inVec = in->data.C32;
     599
     600        out = psVectorRecycle( out, PS_TYPE_F32, numElements );
     601        out->n = numElements;
     602        outVec = out->data.F32;
     603
     604        for ( unsigned int i = 0;i < numElements;i++ ) {
     605            outVec[ i ] = cimagf( inVec[ i ] );
     606        }
     607    } else {
     608        psError( __func__, "Can not extract imaginary component from given vector type (%d).",
     609                 type );
     610        psFree( out );
     611        return NULL;
     612    }
     613
    609614    return out;
    610615}
    … …  
    614619    psElemType type;
    615620    unsigned int numElements;
    616    
    617    
     621
     622
    618623    if ( real == NULL || imag == NULL ) {
    619             psFree( out );
    620             return NULL;
    621         }
    622        
     624        psFree( out );
     625        return NULL;
     626    }
     627
    623628    type = real->type.type;
    624629    if ( real->n < imag->n ) {
    625             numElements = real->n;
    626         } else {
    627             numElements = imag->n;
    628         }
    629        
     630        numElements = real->n;
     631    } else {
     632        numElements = imag->n;
     633    }
     634
    630635    if ( imag->type.type != type ) {
    631             psError( __func__, "The inputs to psVectorComplex must be the same type." );
    632             psFree( out );
    633             return NULL;
    634         }
    635        
     636        psError( __func__, "The inputs to psVectorComplex must be the same type." );
     637        psFree( out );
     638        return NULL;
     639    }
     640
    636641    if ( PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    637             psError( __func__, "The inputs to psVectorComplex can not be complex." );
    638             psFree( out );
    639             return NULL;
    640         }
    641        
     642        psError( __func__, "The inputs to psVectorComplex can not be complex." );
     643        psFree( out );
     644        return NULL;
     645    }
     646
    642647    if ( type == PS_TYPE_F32 ) {
    643             psC32 * outVec;
    644             psF32* realVec = real->data.F32;
    645             psF32* imagVec = imag->data.F32;
    646            
    647             out = psVectorRecycle( out, PS_TYPE_C32, numElements );
    648             out->n = numElements;
    649             outVec = out->data.C32;
    650            
    651             for ( unsigned int i = 0;i < numElements;i++ ) {
    652                     outVec[ i ] = realVec[ i ] + I * imagVec[ i ];
    653                 }
    654         } else {
    655             psError( __func__, "Can not merge real and imaginary portions for given vector type (%d).",
    656                      type );
    657             psFree( out );
    658             return NULL;
    659         }
    660        
     648        psC32 * outVec;
     649        psF32* realVec = real->data.F32;
     650        psF32* imagVec = imag->data.F32;
     651
     652        out = psVectorRecycle( out, PS_TYPE_C32, numElements );
     653        out->n = numElements;
     654        outVec = out->data.C32;
     655
     656        for ( unsigned int i = 0;i < numElements;i++ ) {
     657            outVec[ i ] = realVec[ i ] + I * imagVec[ i ];
     658        }
     659    } else {
     660        psError( __func__, "Can not merge real and imaginary portions for given vector type (%d).",
     661                 type );
     662        psFree( out );
     663        return NULL;
     664    }
     665
    661666    return out;
    662667}
    … …  
    666671    psElemType type;
    667672    unsigned int numElements;
    668    
    669    
    670     if ( in == NULL ) {
    671             psFree( out );
    672             return NULL;
    673         }
    674        
     673
     674
     675    if ( in == NULL ) {
     676        psFree( out );
     677        return NULL;
     678    }
     679
    675680    type = in->type.type;
    676681    numElements = in->n;
    677    
     682
    678683    /* if not a complex number, this is logically just a image copy */
    679684    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    680             // Warn user, as this is probably not expected
    681             psLogMsg( __func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
    682                       "Vector copy was performed instead." );
    683                      
    684             out = psVectorRecycle( out, type, numElements );
    685             out->n = numElements;
    686             memcpy( out->data.V, in->data.V, PSELEMTYPE_SIZEOF( type ) * numElements );
    687             return out;
    688         }
    689        
     685        // Warn user, as this is probably not expected
     686        psLogMsg( __func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
     687                  "Vector copy was performed instead." );
     688
     689        out = psVectorRecycle( out, type, numElements );
     690        out->n = numElements;
     691        memcpy( out->data.V, in->data.V, PSELEMTYPE_SIZEOF( type ) * numElements );
     692        return out;
     693    }
     694
    690695    if ( type == PS_TYPE_C32 ) {
    691             psC32 * outVec;
    692             psC32* inVec = in->data.C32;
    693            
    694             out = psVectorRecycle( out, PS_TYPE_C32, numElements );
    695             out->n = numElements;
    696             outVec = out->data.C32;
    697            
    698             for ( unsigned int i = 0;i < numElements;i++ ) {
    699                     outVec[ i ] = crealf( inVec[ i ] ) - I * cimagf( inVec[ i ] );
    700                 }
    701         } else {
    702             psError( __func__, "Can not compute complex conjugate for given vector type (%d).",
    703                      type );
    704             psFree( out );
    705             return NULL;
    706         }
    707        
     696        psC32 * outVec;
     697        psC32* inVec = in->data.C32;
     698
     699        out = psVectorRecycle( out, PS_TYPE_C32, numElements );
     700        out->n = numElements;
     701        outVec = out->data.C32;
     702
     703        for ( unsigned int i = 0;i < numElements;i++ ) {
     704            outVec[ i ] = crealf( inVec[ i ] ) - I * cimagf( inVec[ i ] );
     705        }
     706    } else {
     707        psError( __func__, "Can not compute complex conjugate for given vector type (%d).",
     708                 type );
     709        psFree( out );
     710        return NULL;
     711    }
     712
    708713    return out;
    709714}
    … …  
    716721    unsigned int inHalfNumElements;
    717722    unsigned int inNumElementsSquared;
    718    
    719     if ( in == NULL ) {
    720             psFree( out );
    721             return NULL;
    722         }
    723        
     723
     724    if ( in == NULL ) {
     725        psFree( out );
     726        return NULL;
     727    }
     728
    724729    type = in->type.type;
    725730    inNumElements = in->n;
    … …  
    727732    inHalfNumElements = inNumElements / 2;
    728733    outNumElements = inHalfNumElements + 1;
    729    
     734
    730735    /* if not a complex number, this is not implemented */
    731736    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
    732             psError( __func__, "Power Spectrum for non-complex inputs is not implemented." );
    733             psFree( out );
    734             return NULL;
    735         }
    736        
     737        psError( __func__, "Power Spectrum for non-complex inputs is not implemented." );
     738        psFree( out );
     739        return NULL;
     740    }
     741
    737742    if ( type == PS_TYPE_C32 ) {
    738             psF32 * outVec;
    739             psC32* inVec = in->data.C32;
    740             psF32 inAbs1;
    741             psF32 inAbs2;
    742            
    743             out = psVectorRecycle( out, PS_TYPE_F32, outNumElements );
    744             out->n = outNumElements;
    745             outVec = out->data.F32;
    746            
    747             // from ADD: P_0 = |C_0|^2/N^2
    748             inAbs1 = cabsf( inVec[ 0 ] );
    749             outVec[ 0 ] = inAbs1 * inAbs1 / inNumElementsSquared;
    750            
    751             // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
    752             for ( unsigned int i = 1;i < inHalfNumElements;i++ ) {
    753                     inAbs1 = cabsf( inVec[ i ] );
    754                     inAbs2 = cabsf( inVec[ inNumElements - i ] );
    755                     outVec[ i ] = ( inAbs1 * inAbs1 + inAbs2 * inAbs2 ) / inNumElementsSquared;
    756                 }
    757                
    758             // from ADD: P_N/2 = |C_N/2|^2/N^2
    759             inAbs1 = cabsf( inVec[ inHalfNumElements ] );
    760             outVec[ inHalfNumElements ] = inAbs1 * inAbs1 / inNumElementsSquared;
    761         } else {
    762             psError( __func__, "Can not power spectrum for given vector type (%d).",
    763                      type );
    764             psFree( out );
    765             return NULL;
    766         }
    767        
    768     return out;
    769    
    770 }
     743        psF32 * outVec;
     744        psC32* inVec = in->data.C32;
     745        psF32 inAbs1;
     746        psF32 inAbs2;
     747
     748        out = psVectorRecycle( out, PS_TYPE_F32, outNumElements );
     749        out->n = outNumElements;
     750        outVec = out->data.F32;
     751
     752        // from ADD: P_0 = |C_0|^2/N^2
     753        inAbs1 = cabsf( inVec[ 0 ] );
     754        outVec[ 0 ] = inAbs1 * inAbs1 / inNumElementsSquared;
     755
     756        // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
     757        for ( unsigned int i = 1;i < inHalfNumElements;i++ ) {
     758            inAbs1 = cabsf( inVec[ i ] );
     759            inAbs2 = cabsf( inVec[ inNumElements - i ] );
     760            outVec[ i ] = ( inAbs1 * inAbs1 + inAbs2 * inAbs2 ) / inNumElementsSquared;
     761        }
     762
     763        // from ADD: P_N/2 = |C_N/2|^2/N^2
     764        inAbs1 = cabsf( inVec[ inHalfNumElements ] );
     765        outVec[ inHalfNumElements ] = inAbs1 * inAbs1 / inNumElementsSquared;
     766    } else {
     767        psError( __func__, "Can not power spectrum for given vector type (%d).",
     768                 type );
     769        psFree( out );
     770        return NULL;
     771    }
     772
     773    return out;
     774
     775}
Note: See TracChangeset for help on using the changeset viewer.