CIMPL: Analysis.cpp Source File

00001 //Copyright (c) 2004-2005, Baris Sumengen
00002 //All rights reserved.
00003 //
00004 // CIMPL Matrix Performance Library
00005 //
00006 //Redistribution and use in source and binary
00007 //forms, with or without modification, are
00008 //permitted provided that the following
00009 //conditions are met:
00010 //
00011 //    * No commercial use is allowed. 
00012 //    This software can only be used
00013 //    for non-commercial purposes. This 
00014 //    distribution is mainly intended for
00015 //    academic research and teaching.
00016 //    * Redistributions of source code must
00017 //    retain the above copyright notice, this
00018 //    list of conditions and the following
00019 //    disclaimer.
00020 //    * Redistributions of binary form must
00021 //    mention the above copyright notice, this
00022 //    list of conditions and the following
00023 //    disclaimer in a clearly visible part 
00024 //    in associated product manual, 
00025 //    readme, and web site of the redistributed 
00026 //    software.
00027 //    * Redistributions in binary form must
00028 //    reproduce the above copyright notice,
00029 //    this list of conditions and the
00030 //    following disclaimer in the
00031 //    documentation and/or other materials
00032 //    provided with the distribution.
00033 //    * The name of Baris Sumengen may not be
00034 //    used to endorse or promote products
00035 //    derived from this software without
00036 //    specific prior written permission.
00037 //
00038 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
00039 //HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
00040 //EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
00041 //NOT LIMITED TO, THE IMPLIED WARRANTIES OF
00042 //MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00043 //PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
00044 //CONTRIBUTORS BE LIABLE FOR ANY
00045 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00046 //EXEMPLARY, OR CONSEQUENTIAL DAMAGES
00047 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
00048 //OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
00049 //DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00050 //HOWEVER CAUSED AND ON ANY THEORY OF
00051 //LIABILITY, WHETHER IN CONTRACT, STRICT
00052 //LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
00053 //OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
00054 //OF THIS SOFTWARE, EVEN IF ADVISED OF THE
00055 //POSSIBILITY OF SUCH DAMAGE.
00056 
00057 
00058 
00059 #include "./Analysis.h"
00060 #include "mkl.h"
00061 
00062 
00063 namespace Analysis
00064 {
00065 
00066         void StatusDisplay(long status)
00067         {
00068                 long classError;                                        
00069                 
00070                 classError = DftiErrorClass(status, DFTI_ERROR_CLASS);  
00071                 if(! classError)
00072                 {
00073                         cerr << "Error Status is not a member of Predefined Error Class\n"; 
00074                 }
00075                 else
00076                 {
00077                         char* errorMessage = DftiErrorMessage(status);
00078                         cerr << "Error_message = " << errorMessage << endl;
00079                 }
00080         }
00081 
00082 
00083         Vector<ComplexFloat> FFT(Vector<ComplexFloat>& x)
00084         {
00085                 int n = x.Length();
00086                 Vector<ComplexFloat> y(n);
00087                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00088                 long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
00089                 if(status != DFTI_NO_ERROR)
00090                 {
00091                         StatusDisplay(status);
00092                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00093                         Utility::RunTimeError("Problem with FFT!");
00094                 }
00095                 // Not Inplace
00096                 status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00097                 if(status != DFTI_NO_ERROR)
00098                 {
00099                         StatusDisplay(status);
00100                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00101                         Utility::RunTimeError("Problem with FFT!");
00102                 }
00103 
00104                 // Commit Dfti descriptor
00105                 status = DftiCommitDescriptor(DescHandle);
00106                 if(status != DFTI_NO_ERROR)
00107                 {
00108                         StatusDisplay(status);
00109                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00110                         Utility::RunTimeError("Problem with FFT!");
00111                 }
00112 
00113                 // Compute Forward transform
00114                 status = DftiComputeForward( DescHandle, x.Data(), y.Data());
00115                 if(status != DFTI_NO_ERROR)
00116                 {
00117                         StatusDisplay(status);
00118                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00119                         Utility::RunTimeError("Problem with FFT!");
00120                 }
00121                 status = DftiFreeDescriptor(&DescHandle);
00122                 if(status != DFTI_NO_ERROR)
00123                 {
00124                         StatusDisplay(status);
00125                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00126                 }
00127 
00128 
00129                 return y;
00130         }
00131 
00132         Vector<ComplexFloat>& FFTI(Vector<ComplexFloat>& x)
00133         {
00134                 int n = x.Length();
00135                 //Vector<ComplexFloat> y(n);
00136                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00137                 long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
00138                 if(status != DFTI_NO_ERROR)
00139                 {
00140                         StatusDisplay(status);
00141                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00142                         Utility::RunTimeError("Problem with FFT!");
00143                 }
00144                 // Not Inplace
00145                 //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00146                 //if(status != DFTI_NO_ERROR)
00147                 //{
00148                 //      StatusDisplay(status);
00149                 //      cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00150                 //      Utility::RunTimeError("Problem with FFT!");
00151                 //}
00152 
00153                 // Commit Dfti descriptor
00154                 status = DftiCommitDescriptor(DescHandle);
00155                 if(status != DFTI_NO_ERROR)
00156                 {
00157                         StatusDisplay(status);
00158                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00159                         Utility::RunTimeError("Problem with FFT!");
00160                 }
00161 
00162                 // Compute Forward transform
00163                 status = DftiComputeForward( DescHandle, x.Data());
00164                 if(status != DFTI_NO_ERROR)
00165                 {
00166                         StatusDisplay(status);
00167                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00168                         Utility::RunTimeError("Problem with FFT!");
00169                 }
00170                 status = DftiFreeDescriptor(&DescHandle);
00171                 if(status != DFTI_NO_ERROR)
00172                 {
00173                         StatusDisplay(status);
00174                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00175                 }
00176 
00177 
00178                 return x;
00179         }
00180 
00181 
00182 
00183 
00184         Vector<ComplexDouble> FFT(Vector<ComplexDouble>& x)
00185         {
00186                 int n = x.Length();
00187                 Vector<ComplexDouble> y(n);
00188                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00189                 long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
00190                 if(status != DFTI_NO_ERROR)
00191                 {
00192                         StatusDisplay(status);
00193                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00194                         Utility::RunTimeError("Problem with FFT!");
00195                 }
00196                 // Not Inplace
00197                 status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00198                 if(status != DFTI_NO_ERROR)
00199                 {
00200                         StatusDisplay(status);
00201                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00202                         Utility::RunTimeError("Problem with FFT!");
00203                 }
00204 
00205                 // Commit Dfti descriptor
00206                 status = DftiCommitDescriptor(DescHandle);
00207                 if(status != DFTI_NO_ERROR)
00208                 {
00209                         StatusDisplay(status);
00210                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00211                         Utility::RunTimeError("Problem with FFT!");
00212                 }
00213 
00214                 // Compute Forward transform
00215                 status = DftiComputeForward( DescHandle, x.Data(), y.Data());
00216                 if(status != DFTI_NO_ERROR)
00217                 {
00218                         StatusDisplay(status);
00219                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00220                         Utility::RunTimeError("Problem with FFT!");
00221                 }
00222                 status = DftiFreeDescriptor(&DescHandle);
00223                 if(status != DFTI_NO_ERROR)
00224                 {
00225                         StatusDisplay(status);
00226                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00227                 }
00228 
00229                 return y;
00230         }
00231 
00232         Vector<ComplexDouble>& FFTI(Vector<ComplexDouble>& x)
00233         {
00234                 int n = x.Length();
00235                 //Vector<ComplexDouble> y(n);
00236                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00237                 long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
00238                 if(status != DFTI_NO_ERROR)
00239                 {
00240                         StatusDisplay(status);
00241                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00242                         Utility::RunTimeError("Problem with FFT!");
00243                 }
00244                 // Not Inplace
00245                 //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00246                 //if(status != DFTI_NO_ERROR)
00247                 //{
00248                 //      StatusDisplay(status);
00249                 //      cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00250                 //      Utility::RunTimeError("Problem with FFT!");
00251                 //}
00252 
00253                 // Commit Dfti descriptor
00254                 status = DftiCommitDescriptor(DescHandle);
00255                 if(status != DFTI_NO_ERROR)
00256                 {
00257                         StatusDisplay(status);
00258                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00259                         Utility::RunTimeError("Problem with FFT!");
00260                 }
00261 
00262                 // Compute Forward transform
00263                 status = DftiComputeForward( DescHandle, x.Data());
00264                 if(status != DFTI_NO_ERROR)
00265                 {
00266                         StatusDisplay(status);
00267                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00268                         Utility::RunTimeError("Problem with FFT!");
00269                 }
00270                 status = DftiFreeDescriptor(&DescHandle);
00271                 if(status != DFTI_NO_ERROR)
00272                 {
00273                         StatusDisplay(status);
00274                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00275                 }
00276 
00277                 return x;
00278         }
00279 
00280 
00281         Vector<ComplexFloat> FFT(Vector<float>& x)
00282         {
00283                 int n = x.Length();
00284                 Vector<ComplexFloat> y(n);
00285                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00286                 long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_REAL, 1, n);
00287                 if(status != DFTI_NO_ERROR)
00288                 {
00289                         StatusDisplay(status);
00290                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00291                         Utility::RunTimeError("Problem with FFT!");
00292                 }
00293                 // Not Inplace
00294                 status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00295                 if(status != DFTI_NO_ERROR)
00296                 {
00297                         StatusDisplay(status);
00298                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00299                         Utility::RunTimeError("Problem with FFT!");
00300                 }
00301 
00302                 // Commit Dfti descriptor
00303                 status = DftiCommitDescriptor(DescHandle);
00304                 if(status != DFTI_NO_ERROR)
00305                 {
00306                         StatusDisplay(status);
00307                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00308                         Utility::RunTimeError("Problem with FFT!");
00309                 }
00310 
00311                 // Compute Forward transform
00312                 status = DftiComputeForward( DescHandle, x.Data(), y.Data());
00313                 if(status != DFTI_NO_ERROR)
00314                 {
00315                         StatusDisplay(status);
00316                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00317                         Utility::RunTimeError("Problem with FFT!");
00318                 }
00319                 status = DftiFreeDescriptor(&DescHandle);
00320                 if(status != DFTI_NO_ERROR)
00321                 {
00322                         StatusDisplay(status);
00323                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00324                 }
00325                 
00326                 for(int i=1;i<(n+1)/2;i++)
00327                 {
00328                         y[n-i] = conj(y[i]);
00329                 }
00330 
00331                 return y;
00332         }
00333 
00334 
00335 
00336         Vector<ComplexDouble> FFT(Vector<double>& x)
00337         {
00338                 int n = x.Length();
00339                 Vector<ComplexDouble> y(n);
00340                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00341                 long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_REAL, 1, n);
00342                 if(status != DFTI_NO_ERROR)
00343                 {
00344                         StatusDisplay(status);
00345                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00346                         Utility::RunTimeError("Problem with FFT!");
00347                 }
00348                 // Not Inplace
00349                 status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00350                 if(status != DFTI_NO_ERROR)
00351                 {
00352                         StatusDisplay(status);
00353                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00354                         Utility::RunTimeError("Problem with FFT!");
00355                 }
00356 
00357                 // Commit Dfti descriptor
00358                 status = DftiCommitDescriptor(DescHandle);
00359                 if(status != DFTI_NO_ERROR)
00360                 {
00361                         StatusDisplay(status);
00362                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00363                         Utility::RunTimeError("Problem with FFT!");
00364                 }
00365 
00366                 // Compute Forward transform
00367                 status = DftiComputeForward( DescHandle, x.Data(), y.Data());
00368                 if(status != DFTI_NO_ERROR)
00369                 {
00370                         StatusDisplay(status);
00371                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00372                         Utility::RunTimeError("Problem with FFT!");
00373                 }
00374                 status = DftiFreeDescriptor(&DescHandle);
00375                 if(status != DFTI_NO_ERROR)
00376                 {
00377                         StatusDisplay(status);
00378                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00379                 }
00380                 
00381                 for(int i=1;i<(n+1)/2;i++)
00382                 {
00383                         y[n-i] = conj(y[i]);
00384                 }
00385 
00386                 return y;
00387         }
00388 
00389 
00390 
00391 
00392 
00393 
00394 
00395 
00396         // ////////////////////
00397         // Backward transform
00398         // ////////////////////
00399 
00400 
00401         Vector<ComplexFloat> IFFT(Vector<ComplexFloat>& x)
00402         {
00403                 int n = x.Length();
00404                 Vector<ComplexFloat> y(n);
00405                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00406                 long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
00407                 if(status != DFTI_NO_ERROR)
00408                 {
00409                         StatusDisplay(status);
00410                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00411                         Utility::RunTimeError("Problem with FFT!");
00412                 }
00413                 // Not Inplace
00414                 status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00415                 if(status != DFTI_NO_ERROR)
00416                 {
00417                         StatusDisplay(status);
00418                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00419                         Utility::RunTimeError("Problem with FFT!");
00420                 }
00421                 // Backward scale
00422                 status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
00423                 if(status != DFTI_NO_ERROR)
00424                 {
00425                         StatusDisplay(status);
00426                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00427                         Utility::RunTimeError("Problem with FFT!");
00428                 }
00429 
00430                 // Commit Dfti descriptor
00431                 status = DftiCommitDescriptor(DescHandle);
00432                 if(status != DFTI_NO_ERROR)
00433                 {
00434                         StatusDisplay(status);
00435                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00436                         Utility::RunTimeError("Problem with FFT!");
00437                 }
00438 
00439                 // Compute Backward transform
00440                 status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
00441                 if(status != DFTI_NO_ERROR)
00442                 {
00443                         StatusDisplay(status);
00444                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00445                         Utility::RunTimeError("Problem with FFT!");
00446                 }
00447                 status = DftiFreeDescriptor(&DescHandle);
00448                 if(status != DFTI_NO_ERROR)
00449                 {
00450                         StatusDisplay(status);
00451                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00452                 }
00453 
00454 
00455                 return y;
00456         }
00457 
00458         Vector<ComplexFloat>& IFFTI(Vector<ComplexFloat>& x)
00459         {
00460                 int n = x.Length();
00461                 //Vector<ComplexFloat> y(n);
00462                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00463                 long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 1, n);
00464                 if(status != DFTI_NO_ERROR)
00465                 {
00466                         StatusDisplay(status);
00467                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00468                         Utility::RunTimeError("Problem with FFT!");
00469                 }
00470                 // Not Inplace
00471                 //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00472                 //if(status != DFTI_NO_ERROR)
00473                 //{
00474                 //      StatusDisplay(status);
00475                 //      cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00476                 //      Utility::RunTimeError("Problem with FFT!");
00477                 //}
00478                 // Backward scale
00479                 status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
00480                 if(status != DFTI_NO_ERROR)
00481                 {
00482                         StatusDisplay(status);
00483                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00484                         Utility::RunTimeError("Problem with FFT!");
00485                 }
00486 
00487                 // Commit Dfti descriptor
00488                 status = DftiCommitDescriptor(DescHandle);
00489                 if(status != DFTI_NO_ERROR)
00490                 {
00491                         StatusDisplay(status);
00492                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00493                         Utility::RunTimeError("Problem with FFT!");
00494                 }
00495 
00496                 // Compute Backward transform
00497                 status = DftiComputeBackward( DescHandle, x.Data());
00498                 if(status != DFTI_NO_ERROR)
00499                 {
00500                         StatusDisplay(status);
00501                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00502                         Utility::RunTimeError("Problem with FFT!");
00503                 }
00504                 status = DftiFreeDescriptor(&DescHandle);
00505                 if(status != DFTI_NO_ERROR)
00506                 {
00507                         StatusDisplay(status);
00508                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00509                 }
00510 
00511 
00512                 return x;
00513         }
00514 
00515 
00516 
00517         Vector<ComplexDouble> IFFT(Vector<ComplexDouble>& x)
00518         {
00519                 int n = x.Length();
00520                 Vector<ComplexDouble> y(n);
00521                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00522                 long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
00523                 if(status != DFTI_NO_ERROR)
00524                 {
00525                         StatusDisplay(status);
00526                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00527                         Utility::RunTimeError("Problem with FFT!");
00528                 }
00529                 // Not Inplace
00530                 status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00531                 if(status != DFTI_NO_ERROR)
00532                 {
00533                         StatusDisplay(status);
00534                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00535                         Utility::RunTimeError("Problem with FFT!");
00536                 }
00537                 // Backward scale
00538                 status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
00539                 if(status != DFTI_NO_ERROR)
00540                 {
00541                         StatusDisplay(status);
00542                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00543                         Utility::RunTimeError("Problem with FFT!");
00544                 }
00545 
00546                 // Commit Dfti descriptor
00547                 status = DftiCommitDescriptor(DescHandle);
00548                 if(status != DFTI_NO_ERROR)
00549                 {
00550                         StatusDisplay(status);
00551                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00552                         Utility::RunTimeError("Problem with FFT!");
00553                 }
00554 
00555                 // Compute Backward transform
00556                 status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
00557                 if(status != DFTI_NO_ERROR)
00558                 {
00559                         StatusDisplay(status);
00560                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00561                         Utility::RunTimeError("Problem with FFT!");
00562                 }
00563                 status = DftiFreeDescriptor(&DescHandle);
00564                 if(status != DFTI_NO_ERROR)
00565                 {
00566                         StatusDisplay(status);
00567                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00568                 }
00569 
00570                 return y;
00571         }
00572 
00573 
00574 
00575 
00576         Vector<ComplexDouble>& IFFTI(Vector<ComplexDouble>& x)
00577         {
00578                 int n = x.Length();
00579                 //Vector<ComplexDouble> y(n);
00580                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00581                 long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n);
00582                 if(status != DFTI_NO_ERROR)
00583                 {
00584                         StatusDisplay(status);
00585                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00586                         Utility::RunTimeError("Problem with FFT!");
00587                 }
00588                 // Not Inplace
00589                 //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00590                 //if(status != DFTI_NO_ERROR)
00591                 //{
00592                 //      StatusDisplay(status);
00593                 //      cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00594                 //      Utility::RunTimeError("Problem with FFT!");
00595                 //}
00596                 // Backward scale
00597                 status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)n);
00598                 if(status != DFTI_NO_ERROR)
00599                 {
00600                         StatusDisplay(status);
00601                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00602                         Utility::RunTimeError("Problem with FFT!");
00603                 }
00604 
00605                 // Commit Dfti descriptor
00606                 status = DftiCommitDescriptor(DescHandle);
00607                 if(status != DFTI_NO_ERROR)
00608                 {
00609                         StatusDisplay(status);
00610                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00611                         Utility::RunTimeError("Problem with FFT!");
00612                 }
00613 
00614                 // Compute Backward transform
00615                 status = DftiComputeBackward( DescHandle, x.Data());
00616                 if(status != DFTI_NO_ERROR)
00617                 {
00618                         StatusDisplay(status);
00619                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00620                         Utility::RunTimeError("Problem with FFT!");
00621                 }
00622                 status = DftiFreeDescriptor(&DescHandle);
00623                 if(status != DFTI_NO_ERROR)
00624                 {
00625                         StatusDisplay(status);
00626                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00627                 }
00628 
00629                 return x;
00630         }
00631 
00632 
00633 
00634 
00635 
00636 
00637 
00638 
00639 
00640 
00641 
00642 
00643 
00644 
00645 
00646 
00647         // /////////////////
00648         // 2D Transform
00649         // /////////////////
00650 
00651 
00652         Matrix<ComplexFloat> FFT2(Matrix<ComplexFloat>& x)
00653         {
00654                 int rows = x.Rows();
00655                 int cols = x.Columns();
00656                 Matrix<ComplexFloat> y(rows,cols);
00657                 int dims[2] = {cols,rows};
00658                 int strides[3] = {0,rows,1};
00659                 
00660                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00661                 long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
00662                 if(status != DFTI_NO_ERROR)
00663                 {
00664                         StatusDisplay(status);
00665                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00666                         Utility::RunTimeError("Problem with FFT!");
00667                 }
00668                 // Strides
00669                 status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
00670                 if(status != DFTI_NO_ERROR)
00671                 {
00672                         StatusDisplay(status);
00673                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00674                         Utility::RunTimeError("Problem with FFT!");
00675                 }
00676                 // Not Inplace
00677                 status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00678                 if(status != DFTI_NO_ERROR)
00679                 {
00680                         StatusDisplay(status);
00681                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00682                         Utility::RunTimeError("Problem with FFT!");
00683                 }
00684 
00685                 // Commit Dfti descriptor
00686                 status = DftiCommitDescriptor(DescHandle);
00687                 if(status != DFTI_NO_ERROR)
00688                 {
00689                         StatusDisplay(status);
00690                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00691                         Utility::RunTimeError("Problem with FFT!");
00692                 }
00693 
00694                 // Compute Forward transform
00695                 status = DftiComputeForward( DescHandle, x.Data(), y.Data());
00696                 if(status != DFTI_NO_ERROR)
00697                 {
00698                         StatusDisplay(status);
00699                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00700                         Utility::RunTimeError("Problem with FFT!");
00701                 }
00702                 status = DftiFreeDescriptor(&DescHandle);
00703                 if(status != DFTI_NO_ERROR)
00704                 {
00705                         StatusDisplay(status);
00706                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00707                 }
00708 
00709 
00710                 return y;
00711         }
00712 
00713 
00714         Matrix<ComplexFloat>& FFT2I(Matrix<ComplexFloat>& x)
00715         {
00716                 int rows = x.Rows();
00717                 int cols = x.Columns();
00718                 //Matrix<ComplexFloat> y(rows,cols);
00719                 int dims[2] = {cols,rows};
00720                 int strides[3] = {0,rows,1};
00721                 
00722                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00723                 long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
00724                 if(status != DFTI_NO_ERROR)
00725                 {
00726                         StatusDisplay(status);
00727                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00728                         Utility::RunTimeError("Problem with FFT!");
00729                 }
00730                 // Strides
00731                 status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
00732                 if(status != DFTI_NO_ERROR)
00733                 {
00734                         StatusDisplay(status);
00735                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00736                         Utility::RunTimeError("Problem with FFT!");
00737                 }
00738                 // Not Inplace
00739                 //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00740                 //if(status != DFTI_NO_ERROR)
00741                 //{
00742                 //      StatusDisplay(status);
00743                 //      cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00744                 //      Utility::RunTimeError("Problem with FFT!");
00745                 //}
00746 
00747                 // Commit Dfti descriptor
00748                 status = DftiCommitDescriptor(DescHandle);
00749                 if(status != DFTI_NO_ERROR)
00750                 {
00751                         StatusDisplay(status);
00752                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00753                         Utility::RunTimeError("Problem with FFT!");
00754                 }
00755 
00756                 // Compute Forward transform
00757                 status = DftiComputeForward( DescHandle, x.Data());
00758                 if(status != DFTI_NO_ERROR)
00759                 {
00760                         StatusDisplay(status);
00761                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00762                         Utility::RunTimeError("Problem with FFT!");
00763                 }
00764                 status = DftiFreeDescriptor(&DescHandle);
00765                 if(status != DFTI_NO_ERROR)
00766                 {
00767                         StatusDisplay(status);
00768                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00769                 }
00770 
00771 
00772                 return x;
00773         }
00774 
00775 
00776 
00777 
00778         Matrix<ComplexDouble> FFT2(Matrix<ComplexDouble>& x)
00779         {
00780                 int rows = x.Rows();
00781                 int cols = x.Columns();
00782                 Matrix<ComplexDouble> y(rows,cols);
00783                 int dims[2] = {cols,rows};
00784                 int strides[3] = {0,rows,1};
00785                 
00786                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00787                 long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
00788                 if(status != DFTI_NO_ERROR)
00789                 {
00790                         StatusDisplay(status);
00791                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00792                         Utility::RunTimeError("Problem with FFT!");
00793                 }
00794                 // Strides
00795                 status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
00796                 if(status != DFTI_NO_ERROR)
00797                 {
00798                         StatusDisplay(status);
00799                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00800                         Utility::RunTimeError("Problem with FFT!");
00801                 }
00802                 // Not Inplace
00803                 status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00804                 if(status != DFTI_NO_ERROR)
00805                 {
00806                         StatusDisplay(status);
00807                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00808                         Utility::RunTimeError("Problem with FFT!");
00809                 }
00810 
00811                 // Commit Dfti descriptor
00812                 status = DftiCommitDescriptor(DescHandle);
00813                 if(status != DFTI_NO_ERROR)
00814                 {
00815                         StatusDisplay(status);
00816                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00817                         Utility::RunTimeError("Problem with FFT!");
00818                 }
00819 
00820                 // Compute Forward transform
00821                 status = DftiComputeForward( DescHandle, x.Data(), y.Data());
00822                 if(status != DFTI_NO_ERROR)
00823                 {
00824                         StatusDisplay(status);
00825                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00826                         Utility::RunTimeError("Problem with FFT!");
00827                 }
00828                 status = DftiFreeDescriptor(&DescHandle);
00829                 if(status != DFTI_NO_ERROR)
00830                 {
00831                         StatusDisplay(status);
00832                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00833                 }
00834 
00835 
00836                 return y;
00837         }
00838 
00839 
00840         Matrix<ComplexDouble>& FFT2I(Matrix<ComplexDouble>& x)
00841         {
00842                 int rows = x.Rows();
00843                 int cols = x.Columns();
00844                 //Matrix<ComplexDouble> y(rows,cols);
00845                 int dims[2] = {cols,rows};
00846                 int strides[3] = {0,rows,1};
00847                 
00848                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00849                 long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
00850                 if(status != DFTI_NO_ERROR)
00851                 {
00852                         StatusDisplay(status);
00853                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00854                         Utility::RunTimeError("Problem with FFT!");
00855                 }
00856                 // Strides
00857                 status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
00858                 if(status != DFTI_NO_ERROR)
00859                 {
00860                         StatusDisplay(status);
00861                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00862                         Utility::RunTimeError("Problem with FFT!");
00863                 }
00864                 // Not Inplace
00865                 //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00866                 //if(status != DFTI_NO_ERROR)
00867                 //{
00868                 //      StatusDisplay(status);
00869                 //      cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00870                 //      Utility::RunTimeError("Problem with FFT!");
00871                 //}
00872 
00873                 // Commit Dfti descriptor
00874                 status = DftiCommitDescriptor(DescHandle);
00875                 if(status != DFTI_NO_ERROR)
00876                 {
00877                         StatusDisplay(status);
00878                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00879                         Utility::RunTimeError("Problem with FFT!");
00880                 }
00881 
00882                 // Compute Forward transform
00883                 status = DftiComputeForward( DescHandle, x.Data() );
00884                 if(status != DFTI_NO_ERROR)
00885                 {
00886                         StatusDisplay(status);
00887                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00888                         Utility::RunTimeError("Problem with FFT!");
00889                 }
00890                 status = DftiFreeDescriptor(&DescHandle);
00891                 if(status != DFTI_NO_ERROR)
00892                 {
00893                         StatusDisplay(status);
00894                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00895                 }
00896 
00897 
00898                 return x;
00899         }
00900 
00901 
00902 
00903 
00904 
00905 
00906 
00907 
00908 
00909 
00910 
00911 
00912         Matrix<ComplexFloat> IFFT2(Matrix<ComplexFloat>& x)
00913         {
00914                 int rows = x.Rows();
00915                 int cols = x.Columns();
00916                 Matrix<ComplexFloat> y(rows,cols);
00917                 int dims[2] = {cols,rows};
00918                 int strides[3] = {0,rows,1};
00919                 
00920                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00921                 long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
00922                 if(status != DFTI_NO_ERROR)
00923                 {
00924                         StatusDisplay(status);
00925                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00926                         Utility::RunTimeError("Problem with FFT!");
00927                 }
00928                 // Strides
00929                 status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
00930                 if(status != DFTI_NO_ERROR)
00931                 {
00932                         StatusDisplay(status);
00933                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00934                         Utility::RunTimeError("Problem with FFT!");
00935                 }
00936                 // Not Inplace
00937                 status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00938                 if(status != DFTI_NO_ERROR)
00939                 {
00940                         StatusDisplay(status);
00941                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00942                         Utility::RunTimeError("Problem with FFT!");
00943                 }
00944                 // Backward scale
00945                 status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
00946                 if(status != DFTI_NO_ERROR)
00947                 {
00948                         StatusDisplay(status);
00949                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00950                         Utility::RunTimeError("Problem with FFT!");
00951                 }
00952 
00953                 // Commit Dfti descriptor
00954                 status = DftiCommitDescriptor(DescHandle);
00955                 if(status != DFTI_NO_ERROR)
00956                 {
00957                         StatusDisplay(status);
00958                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00959                         Utility::RunTimeError("Problem with FFT!");
00960                 }
00961 
00962                 // Compute Forward transform
00963                 status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
00964                 if(status != DFTI_NO_ERROR)
00965                 {
00966                         StatusDisplay(status);
00967                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00968                         Utility::RunTimeError("Problem with FFT!");
00969                 }
00970                 status = DftiFreeDescriptor(&DescHandle);
00971                 if(status != DFTI_NO_ERROR)
00972                 {
00973                         StatusDisplay(status);
00974                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
00975                 }
00976 
00977 
00978                 return y;
00979         }
00980 
00981 
00982         Matrix<ComplexFloat>& IFFT2I(Matrix<ComplexFloat>& x)
00983         {
00984                 int rows = x.Rows();
00985                 int cols = x.Columns();
00986                 //Matrix<ComplexFloat> y(rows,cols);
00987                 int dims[2] = {cols,rows};
00988                 int strides[3] = {0,rows,1};
00989                 
00990                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
00991                 long status = DftiCreateDescriptor(&DescHandle, DFTI_SINGLE, DFTI_COMPLEX, 2, dims);
00992                 if(status != DFTI_NO_ERROR)
00993                 {
00994                         StatusDisplay(status);
00995                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00996                         Utility::RunTimeError("Problem with FFT!");
00997                 }
00998                 // Strides
00999                 status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
01000                 if(status != DFTI_NO_ERROR)
01001                 {
01002                         StatusDisplay(status);
01003                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01004                         Utility::RunTimeError("Problem with FFT!");
01005                 }
01006                 // Not Inplace
01007                 //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
01008                 //if(status != DFTI_NO_ERROR)
01009                 //{
01010                 //      StatusDisplay(status);
01011                 //      cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01012                 //      Utility::RunTimeError("Problem with FFT!");
01013                 //}
01014                 // Backward scale
01015                 status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
01016                 if(status != DFTI_NO_ERROR)
01017                 {
01018                         StatusDisplay(status);
01019                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01020                         Utility::RunTimeError("Problem with FFT!");
01021                 }
01022 
01023                 // Commit Dfti descriptor
01024                 status = DftiCommitDescriptor(DescHandle);
01025                 if(status != DFTI_NO_ERROR)
01026                 {
01027                         StatusDisplay(status);
01028                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01029                         Utility::RunTimeError("Problem with FFT!");
01030                 }
01031 
01032                 // Compute Forward transform
01033                 status = DftiComputeBackward( DescHandle, x.Data());
01034                 if(status != DFTI_NO_ERROR)
01035                 {
01036                         StatusDisplay(status);
01037                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01038                         Utility::RunTimeError("Problem with FFT!");
01039                 }
01040                 status = DftiFreeDescriptor(&DescHandle);
01041                 if(status != DFTI_NO_ERROR)
01042                 {
01043                         StatusDisplay(status);
01044                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
01045                 }
01046 
01047 
01048                 return x;
01049         }
01050 
01051 
01052 
01053 
01054         Matrix<ComplexDouble> IFFT2(Matrix<ComplexDouble>& x)
01055         {
01056                 int rows = x.Rows();
01057                 int cols = x.Columns();
01058                 Matrix<ComplexDouble> y(rows,cols);
01059                 int dims[2] = {cols,rows};
01060                 int strides[3] = {0,rows,1};
01061                 
01062                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
01063                 long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
01064                 if(status != DFTI_NO_ERROR)
01065                 {
01066                         StatusDisplay(status);
01067                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01068                         Utility::RunTimeError("Problem with FFT!");
01069                 }
01070                 // Strides
01071                 status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
01072                 if(status != DFTI_NO_ERROR)
01073                 {
01074                         StatusDisplay(status);
01075                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01076                         Utility::RunTimeError("Problem with FFT!");
01077                 }
01078                 // Not Inplace
01079                 status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
01080                 if(status != DFTI_NO_ERROR)
01081                 {
01082                         StatusDisplay(status);
01083                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01084                         Utility::RunTimeError("Problem with FFT!");
01085                 }
01086                 // Backward scale
01087                 status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
01088                 if(status != DFTI_NO_ERROR)
01089                 {
01090                         StatusDisplay(status);
01091                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01092                         Utility::RunTimeError("Problem with FFT!");
01093                 }
01094 
01095                 // Commit Dfti descriptor
01096                 status = DftiCommitDescriptor(DescHandle);
01097                 if(status != DFTI_NO_ERROR)
01098                 {
01099                         StatusDisplay(status);
01100                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01101                         Utility::RunTimeError("Problem with FFT!");
01102                 }
01103 
01104                 // Compute Forward transform
01105                 status = DftiComputeBackward( DescHandle, x.Data(), y.Data());
01106                 if(status != DFTI_NO_ERROR)
01107                 {
01108                         StatusDisplay(status);
01109                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01110                         Utility::RunTimeError("Problem with FFT!");
01111                 }
01112                 status = DftiFreeDescriptor(&DescHandle);
01113                 if(status != DFTI_NO_ERROR)
01114                 {
01115                         StatusDisplay(status);
01116                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
01117                 }
01118 
01119 
01120                 return y;
01121         }
01122 
01123 
01124 
01125 
01126         Matrix<ComplexDouble>& IFFT2I(Matrix<ComplexDouble>& x)
01127         {
01128                 int rows = x.Rows();
01129                 int cols = x.Columns();
01130                 //Matrix<ComplexDouble> y(rows,cols);
01131                 int dims[2] = {cols,rows};
01132                 int strides[3] = {0,rows,1};
01133                 
01134                 DFTI_DESCRIPTOR_HANDLE DescHandle = 0;
01135                 long status = DftiCreateDescriptor(&DescHandle, DFTI_DOUBLE, DFTI_COMPLEX, 2, dims);
01136                 if(status != DFTI_NO_ERROR)
01137                 {
01138                         StatusDisplay(status);
01139                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01140                         Utility::RunTimeError("Problem with FFT!");
01141                 }
01142                 // Strides
01143                 status = DftiSetValue(DescHandle, DFTI_INPUT_STRIDES, strides);
01144                 if(status != DFTI_NO_ERROR)
01145                 {
01146                         StatusDisplay(status);
01147                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01148                         Utility::RunTimeError("Problem with FFT!");
01149                 }
01150                 // Not Inplace
01151                 //status = DftiSetValue(DescHandle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
01152                 //if(status != DFTI_NO_ERROR)
01153                 //{
01154                 //      StatusDisplay(status);
01155                 //      cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01156                 //      Utility::RunTimeError("Problem with FFT!");
01157                 //}
01158                 // Backward scale
01159                 status = DftiSetValue(DescHandle, DFTI_BACKWARD_SCALE, 1.0/(double)(rows*cols));
01160                 if(status != DFTI_NO_ERROR)
01161                 {
01162                         StatusDisplay(status);
01163                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01164                         Utility::RunTimeError("Problem with FFT!");
01165                 }
01166 
01167                 // Commit Dfti descriptor
01168                 status = DftiCommitDescriptor(DescHandle);
01169                 if(status != DFTI_NO_ERROR)
01170                 {
01171                         StatusDisplay(status);
01172                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01173                         Utility::RunTimeError("Problem with FFT!");
01174                 }
01175 
01176                 // Compute Forward transform
01177                 status = DftiComputeBackward( DescHandle, x.Data());
01178                 if(status != DFTI_NO_ERROR)
01179                 {
01180                         StatusDisplay(status);
01181                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01182                         Utility::RunTimeError("Problem with FFT!");
01183                 }
01184                 status = DftiFreeDescriptor(&DescHandle);
01185                 if(status != DFTI_NO_ERROR)
01186                 {
01187                         StatusDisplay(status);
01188                         Utility::Warning("Problem when trying to free the FFT descriptor handle!");
01189                 }
01190 
01191 
01192                 return x;
01193         }
01194 
01195 
01196 
01197 
01198 
01199 
01200 
01201 
01202 
01203 
01204 
01205 
01206 //=========================================================================
01207 // Convolution
01208 //=========================================================================
01209 
01210 
01211         Vector<float> Conv(Vector<float>& input, Vector<float>& filter, FilterArea fa, Border b, bool PowerOfTwo)
01212         {
01213                 return Real( Conv(ToComplexFloat(input), ToComplexFloat(filter), fa, b, PowerOfTwo) );
01214         }
01215 
01216 
01217         Vector<double> Conv(Vector<double>& input, Vector<double>& filter, FilterArea fa, Border b, bool PowerOfTwo)
01218         {
01219                 return Real( Conv(ToComplexDouble(input), ToComplexDouble(filter), fa, b, PowerOfTwo) );
01220         }
01221 
01222 
01223 
01224         Vector<ComplexFloat> Conv(Vector<ComplexFloat>& input, Vector<ComplexFloat>& filter, FilterArea fa, Border b, bool PowerOfTwo)
01225         {
01226                 // Check if the filter size is smaller than the input size
01227                 if(filter.Length() > input.Length())
01228                 {
01229                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01230                         Utility::RunTimeError("Filter length should not be larger than the input signal!");
01231                 }
01232 
01233                 int fftlength = input.Length()+filter.Length()-1;
01234                 if(PowerOfTwo)
01235                 {
01236                         fftlength = NextPow2(fftlength);
01237                 }
01238                 
01239 
01240                 Vector<ComplexFloat> temp1;
01241                 if(b == ZeroPad || b == Symmetric || b == Replicate) 
01242                 {
01243                         temp1 = Vector<ComplexFloat>(fftlength);
01244                         for(int i=0; i<input.Length(); i++)
01245                         {
01246                                 temp1[i] = input[i];
01247                         }
01248                 }
01249                 else if(b == Circular)
01250                 {
01251                         temp1 = input;
01252                 }
01253                 else
01254                 {
01255                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01256                         Utility::RunTimeError("Invalid Border type!");
01257                 }
01258                 
01259                 
01260                 if(b == Symmetric)
01261                 {
01262                         int wid1 = filter.Length()/2;
01263                         int wid2 = filter.Length() - wid1 - 1;
01264                         if(wid1 != 0)
01265                         {
01266                                 Vector<ComplexFloat> c1 = Reverse( input.Slice(0, wid1-1) );
01267                                 temp1.ReadFromVector(c1, temp1.Length()-wid1);
01268                         }
01269                         if(wid2 != 0)
01270                         {
01271                                 Vector<ComplexFloat> c2 = Reverse( input.Slice(input.Length()-wid2, input.Length()-1) );
01272                                 temp1.ReadFromVector(c2, input.Length());
01273                         }
01274                 }
01275                 else if(b == Replicate)
01276                 {
01277                         int wid1 = filter.Length()/2;
01278                         int wid2 = filter.Length() - wid1 - 1;
01279                         if(wid1 != 0)
01280                         {
01281                                 Vector<ComplexFloat> c1(wid1, input.ElemNC(0));
01282                                 temp1.ReadFromVector(c1, temp1.Length()-wid1);
01283                         }
01284                         if(wid2 != 0)
01285                         {
01286                                 Vector<ComplexFloat> c2(wid2, input.ElemNC(input.Length()-1));
01287                                 temp1.ReadFromVector(c2, input.Length());
01288                         }
01289                 }
01290 
01291                 Vector<ComplexFloat> temp2;
01292                 if(b == ZeroPad || b == Symmetric || b == Replicate)
01293                 {
01294                         temp2 = Vector<ComplexFloat>(fftlength);
01295                         for(int i=0; i<filter.Length(); i++)
01296                         {
01297                                 temp2[i] = filter[i];
01298                         }
01299                 }
01300                 else if(b == Circular)
01301                 {
01302                         temp2 = Vector<ComplexFloat>(input.Length());
01303                         for(int i=0; i<filter.Length(); i++)
01304                         {
01305                                 temp2[i] = filter[i];
01306                         }
01307                 }
01308                 else
01309                 {
01310                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01311                         Utility::RunTimeError("Invalid Border type!");
01312                 }
01313 
01314                 
01315                 Vector<ComplexFloat> o = IFFT( FFT(temp1) * FFT(temp2) );
01316                 
01317                 if(fa == Same)
01318                 {
01319                         o = o.Slice(filter.Length()/2, filter.Length()/2+input.Length()-1);
01320                 }
01321                 else if(fa == Valid)
01322                 {
01323                         o = o.Slice(filter.Length()-1, filter.Length()-1+input.Length()-filter.Length());
01324                 }
01325                 else if(fa != Full)
01326                 {
01327                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01328                         Utility::RunTimeError("Invalid FilterArea type!");
01329                 }
01330                 
01331                 return o;
01332                 
01333         }
01334 
01335         Vector<ComplexDouble> Conv(Vector<ComplexDouble>& input, Vector<ComplexDouble>& filter, FilterArea fa, Border b, bool PowerOfTwo)
01336         {
01337                 // Check if the filter size is smaller than the input size
01338                 if(filter.Length() > input.Length())
01339                 {
01340                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01341                         Utility::RunTimeError("Filter length should not be larger than the input signal!");
01342                 }
01343 
01344                 int fftlength = input.Length()+filter.Length()-1;
01345                 if(PowerOfTwo)
01346                 {
01347                         fftlength = NextPow2(fftlength);
01348                 }
01349                 
01350 
01351                 Vector<ComplexDouble> temp1;
01352                 if(b == ZeroPad || b == Symmetric || b == Replicate) 
01353                 {
01354                         temp1 = Vector<ComplexDouble>(fftlength);
01355                         for(int i=0; i<input.Length(); i++)
01356                         {
01357                                 temp1[i] = input[i];
01358                         }
01359                 }
01360                 else if(b == Circular)
01361                 {
01362                         temp1 = input;
01363                 }
01364                 else
01365                 {
01366                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01367                         Utility::RunTimeError("Invalid Border type!");
01368                 }
01369                 
01370                 
01371                 if(b == Symmetric)
01372                 {
01373                         int wid1 = filter.Length()/2;
01374                         int wid2 = filter.Length() - wid1 - 1;
01375                         if(wid1 != 0)
01376                         {
01377                                 Vector<ComplexDouble> c1 = Reverse( input.Slice(0, wid1-1) );
01378                                 temp1.ReadFromVector(c1, temp1.Length()-wid1);
01379                         }
01380                         if(wid2 != 0)
01381                         {
01382                                 Vector<ComplexDouble> c2 = Reverse( input.Slice(input.Length()-wid2, input.Length()-1) );
01383                                 temp1.ReadFromVector(c2, input.Length());
01384                         }
01385                 }
01386                 else if(b == Replicate)
01387                 {
01388                         int wid1 = filter.Length()/2;
01389                         int wid2 = filter.Length() - wid1 - 1;
01390                         if(wid1 != 0)
01391                         {
01392                                 Vector<ComplexDouble> c1(wid1, input.ElemNC(0));
01393                                 temp1.ReadFromVector(c1, temp1.Length()-wid1);
01394                         }
01395                         if(wid2 != 0)
01396                         {
01397                                 Vector<ComplexDouble> c2(wid2, input.ElemNC(input.Length()-1));
01398                                 temp1.ReadFromVector(c2, input.Length());
01399                         }
01400                 }
01401 
01402                 Vector<ComplexDouble> temp2;
01403                 if(b == ZeroPad || b == Symmetric || b == Replicate)
01404                 {
01405                         temp2 = Vector<ComplexDouble>(fftlength);
01406                         for(int i=0; i<filter.Length(); i++)
01407                         {
01408                                 temp2[i] = filter[i];
01409                         }
01410                 }
01411                 else if(b == Circular)
01412                 {
01413                         temp2 = Vector<ComplexDouble>(input.Length());
01414                         for(int i=0; i<filter.Length(); i++)
01415                         {
01416                                 temp2[i] = filter[i];
01417                         }
01418                 }
01419                 else
01420                 {
01421                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01422                         Utility::RunTimeError("Invalid Border type!");
01423                 }
01424 
01425                 
01426                 Vector<ComplexDouble> o = IFFT( FFT(temp1) * FFT(temp2) );
01427                 
01428                 if(fa == Same)
01429                 {
01430                         o = o.Slice(filter.Length()/2, filter.Length()/2+input.Length()-1);
01431                 }
01432                 else if(fa == Valid)
01433                 {
01434                         o = o.Slice(filter.Length()-1, filter.Length()-1+input.Length()-filter.Length());
01435                 }
01436                 else if(fa != Full)
01437                 {
01438                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01439                         Utility::RunTimeError("Invalid FilterArea type!");
01440                 }
01441                 
01442                 return o;
01443                 
01444         }
01445         
01446         
01447         // =============
01448         // Conv2
01449         // =============
01450 
01451 
01452         Matrix<float> Conv2(Matrix<float>& input, Matrix<float>& filter, FilterArea fa, Border b, bool PowerOfTwo)
01453         {
01454                 return Real( Conv2(ToComplexFloat(input), ToComplexFloat(filter), fa, b, PowerOfTwo) );
01455         }
01456 
01457         Matrix<float> Conv2(Vector<float>& filter1, Vector<float>& filter2, Matrix<float>& input, FilterArea fa, Border b, bool PowerOfTwo)
01458         {
01459                 return Real( Conv2(ToComplexFloat(filter1), ToComplexFloat(filter2), ToComplexFloat(input), fa, b, PowerOfTwo) );
01460         }
01461 
01462         Matrix<double> Conv2(Matrix<double>& input, Matrix<double>& filter, FilterArea fa, Border b, bool PowerOfTwo)
01463         {
01464                 return Real( Conv2(ToComplexDouble(input), ToComplexDouble(filter), fa, b, PowerOfTwo) );
01465         }
01466 
01467         Matrix<double> Conv2(Vector<double>& filter1, Vector<double>& filter2, Matrix<double>& input, FilterArea fa, Border b, bool PowerOfTwo)
01468         {
01469                 return Real( Conv2(ToComplexDouble(filter1), ToComplexDouble(filter2), ToComplexDouble(input), fa, b, PowerOfTwo) );
01470         }
01471 
01472 
01473         
01474         Matrix<ComplexFloat> Conv2(Matrix<ComplexFloat>& input, Matrix<ComplexFloat>& filter, FilterArea fa, Border b, bool PowerOfTwo)
01475         {
01476 
01477                 // Check if the filter size is smaller than the input size
01478                 if(filter.Rows() > input.Rows() || filter.Columns() > input.Columns())
01479                 {
01480                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01481                         Utility::RunTimeError("Filter dimensions should not be larger than the input signal!");
01482                 }
01483                 
01484                 
01485                 // Allocate area for input using the border
01486                 
01487                 int fftrows = input.Rows()+filter.Rows()-1;
01488                 int fftcols = input.Columns()+filter.Columns()-1;
01489                 if(PowerOfTwo)
01490                 {
01491                         fftrows = NextPow2(fftrows);
01492                         fftcols = NextPow2(fftcols);
01493                 }
01494                 
01495                 Matrix<ComplexFloat> temp1;
01496                 if(b == ZeroPad || b == Symmetric || b == Replicate) 
01497                 {
01498                         temp1 = Matrix<ComplexFloat>(fftrows, fftcols,0);
01499                         for(int j=0; j<input.Columns(); j++)
01500                         {
01501                                 for(int i=0; i<input.Rows(); i++)
01502                                 {
01503                                         temp1.ElemNC(i,j) = input.ElemNC(i,j);
01504                                 }
01505                         }
01506                 }
01507                 else if(b == Circular)
01508                 {
01509                         temp1 = input;
01510                 }
01511                 else
01512                 {
01513                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01514                         Utility::RunTimeError("Invalid Border type!");
01515                 }
01516 
01517 
01518                 if(b == Symmetric)
01519                 {
01520                         int wid1 = filter.Columns()/2;
01521                         int wid2 = filter.Columns() - wid1 - 1;
01522                         if(wid1 != 0)
01523                         {
01524                                 Matrix<ComplexFloat> c1 = FlipLRI( input.Slice(0, input.Rows()-1, 0, wid1-1) );
01525                                 temp1.ReadFromMatrix(c1, 0, temp1.Columns()-wid1);
01526                         }
01527                         if(wid2 != 0)
01528                         {
01529                                 Matrix<ComplexFloat> c2 = FlipLRI( input.Slice(0, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) );
01530                                 temp1.ReadFromMatrix(c2, 0, input.Columns());
01531                         }
01532                         
01533                         int hei1 = filter.Rows()/2;
01534                         int hei2 = filter.Rows() - hei1 - 1;
01535                         if(hei1 != 0)
01536                         {
01537                                 Matrix<ComplexFloat> r1 = FlipUDI( input.Slice(0, hei1-1, 0, input.Columns()-1) );
01538                                 temp1.ReadFromMatrix(r1, temp1.Rows()-hei1, 0);
01539                         }
01540                         if(hei2 != 0)
01541                         {
01542                                 Matrix<ComplexFloat> r2 = FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, input.Columns()-1) );
01543                                 temp1.ReadFromMatrix(r2, input.Rows(), 0);
01544                         }
01545                         
01546                         
01547                         if(hei1 != 0 && wid1 != 0)
01548                         {
01549                                 Matrix<ComplexFloat> q11 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, 0, wid1-1) ));
01550                                 temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
01551                         }
01552 
01553                         if(hei1 != 0 && wid2 != 0)
01554                         {
01555                                 Matrix<ComplexFloat> q12 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, input.Columns()-wid2, input.Columns()-1) ));
01556                                 temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
01557                         }
01558 
01559                         if(hei2 != 0 && wid1 != 0)
01560                         {
01561                                 Matrix<ComplexFloat> q21 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, wid1-1) ));
01562                                 temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
01563                         }
01564 
01565                         if(hei2 != 0 && wid2 != 0)
01566                         {
01567                                 Matrix<ComplexFloat> q22 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) ));
01568                                 temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
01569                         }
01570                 }
01571                 else if(b == Replicate)
01572                 {
01573                         int wid1 = filter.Columns()/2;
01574                         int wid2 = filter.Columns() - wid1 - 1;
01575                         if(wid1 != 0)
01576                         {
01577                                 Matrix<ComplexFloat> c1 = input.Slice(0, input.Rows()-1, 0, 0);
01578                                 temp1.ReadFromMatrix(RepMat(c1,1,wid1), 0, temp1.Columns()-wid1);
01579                         }
01580                         if(wid2 != 0)
01581                         {
01582                                 Matrix<ComplexFloat> c2 = input.Slice(0, input.Rows()-1, input.Columns()-1, input.Columns()-1);
01583                                 temp1.ReadFromMatrix(RepMat(c2,1,wid2), 0, input.Columns());
01584                         }
01585                         
01586                         int hei1 = filter.Rows()/2;
01587                         int hei2 = filter.Rows() - hei1 - 1;
01588                         if(hei1 != 0)
01589                         {
01590                                 Matrix<ComplexFloat> r1 = input.Slice(0, 0, 0, input.Columns()-1);
01591                                 temp1.ReadFromMatrix(RepMat(r1,hei1,1), temp1.Rows()-hei1, 0);
01592                         }
01593                         if(hei2 != 0)
01594                         {
01595                                 Matrix<ComplexFloat> r2 = input.Slice(input.Rows()-1, input.Rows()-1, 0, input.Columns()-1);
01596                                 temp1.ReadFromMatrix(RepMat(r2,hei2,1), input.Rows(), 0);
01597                         }
01598                         
01599                         
01600                         if(hei1 != 0 && wid1 != 0)
01601                         {
01602                                 Matrix<ComplexFloat> q11( hei1, wid1, input.ElemNC(0,0) );
01603                                 temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
01604                         }
01605 
01606                         if(hei1 != 0 && wid2 != 0)
01607                         {
01608                                 Matrix<ComplexFloat> q12( hei1, wid2, input.ElemNC(0,input.Columns()-1) );
01609                                 temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
01610                         }
01611 
01612                         if(hei2 != 0 && wid1 != 0)
01613                         {
01614                                 Matrix<ComplexFloat> q21( hei2, wid1, input.ElemNC(input.Rows()-1,0) );
01615                                 temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
01616                         }
01617 
01618                         if(hei2 != 0 && wid2 != 0)
01619                         {
01620                                 Matrix<ComplexFloat> q22(hei2,wid2,input.ElemNC(input.Rows()-1,input.Columns()-1));
01621                                 temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
01622                         }
01623                 }
01624 
01625                 
01626                 // Allocate area for the filter
01627                 Matrix<ComplexFloat> temp2;
01628                 if(b == ZeroPad || b == Symmetric || b == Replicate)
01629                 {
01630                         temp2 = Matrix<ComplexFloat>(fftrows, fftcols,0);
01631                         for(int j=0; j<filter.Columns(); j++)
01632                         {
01633                                 for(int i=0; i<filter.Rows(); i++)
01634                                 {
01635                                         temp2.ElemNC(i,j) = filter.ElemNC(i,j);
01636                                 }
01637                         }
01638                 }
01639                 else if(b == Circular)
01640                 {
01641                         temp2 = Matrix<ComplexFloat>(input.Rows(), input.Columns(),0);
01642                         for(int j=0; j<filter.Columns(); j++)
01643                         {
01644                                 for(int i=0; i<filter.Rows(); i++)
01645                                 {
01646                                         temp2.ElemNC(i,j) = filter.ElemNC(i,j);
01647                                 }
01648                         }
01649                 }
01650                 else
01651                 {
01652                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01653                         Utility::RunTimeError("Invalid Border type!");
01654                 }
01655 
01656 
01657                 // FFT and inverse FFT
01658                 Matrix<ComplexFloat> o = IFFT2( FFT2(temp1)*FFT2(temp2) );
01659 
01660                 // Based on the filter area, crop the border.
01661                 if(fa == Same)
01662                 {
01663                         o = o.Slice(filter.Rows()/2, filter.Rows()/2+input.Rows()-1, filter.Columns()/2, filter.Columns()/2+input.Columns()-1);
01664                 }
01665                 else if(fa == Valid)
01666                 {
01667                         o = o.Slice(filter.Rows()-1, filter.Rows()-1+input.Rows()-filter.Rows(), filter.Columns()-1, filter.Columns()-1+input.Columns()-filter.Columns());
01668                 }
01669                 else if(fa != Full)
01670                 {
01671                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01672                         Utility::RunTimeError("Invalid FilterArea type!");
01673                 }
01674                 
01675                 return o;
01676                 
01677         }
01678 
01679 
01680         
01681         Matrix<ComplexFloat> Conv2(Vector<ComplexFloat>& filter1, Vector<ComplexFloat>& filter2, Matrix<ComplexFloat>& input, FilterArea fa, Border b, bool PowerOfTwo)
01682         {
01683                 // Check if the filter sizes is smaller than the input size
01684                 if(filter1.Length() > input.Rows() || filter2.Length() > input.Columns())
01685                 {
01686                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01687                         Utility::RunTimeError("Filter dimensions should not be larger than the input signal! filter1 <= # of rows, filter2 <= # of columns");
01688                 }
01689                 
01690                 Vector<ComplexFloat> *cols = new (std::nothrow) Vector<ComplexFloat>[input.Columns()];
01691                 Utility::CheckPointer(cols);
01692                 
01693                 for(int i=0; i<input.Columns(); i++)
01694                 {
01695                         cols[i] = Conv(input[i], filter1, fa, b, PowerOfTwo);
01696                 }
01697                 
01698                 int rlength = cols[0].Length();
01699                 
01700                 Vector<ComplexFloat> *rows = new (std::nothrow) Vector<ComplexFloat>[rlength];
01701                 Utility::CheckPointer(rows);
01702                 
01703                 for(int i=0; i<rlength; i++)
01704                 {
01705                         rows[i] = Vector<ComplexFloat>(input.Columns());
01706                         for(int j=0; j<input.Columns(); j++)
01707                         {
01708                                 rows[i][j] = cols[j][i];
01709                         }
01710                 }
01711                 
01712                 delete [] cols;
01713                 
01714                 for(int i=0; i<rlength; i++)
01715                 {
01716                         rows[i] = Conv(rows[i], filter2, fa, b, PowerOfTwo);
01717                 }
01718                 
01719                 int clength = rows[0].Length();
01720                 Matrix<ComplexFloat> temp(rlength, clength);
01721                 for(int i=0; i<rlength; i++)
01722                 {
01723                         for(int j=0; j<clength; j++)
01724                         {
01725                                 temp.ElemNC(i,j) = rows[i][j];
01726                         }
01727                 }
01728                 
01729                 delete [] rows;
01730                 
01731                 return temp;
01732         }
01733 
01734 
01735 
01736 
01737         Matrix<ComplexDouble> Conv2(Matrix<ComplexDouble>& input, Matrix<ComplexDouble>& filter, FilterArea fa, Border b, bool PowerOfTwo)
01738         {
01739                 // Check if the filter size is smaller than the input size
01740                 if(filter.Rows() > input.Rows() || filter.Columns() > input.Columns())
01741                 {
01742                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01743                         Utility::RunTimeError("Filter dimensions should not be larger than the input signal!");
01744                 }
01745                 
01746                 
01747                 // Allocate area for input using the border
01748                 
01749                 int fftrows = input.Rows()+filter.Rows()-1;
01750                 int fftcols = input.Columns()+filter.Columns()-1;
01751                 if(PowerOfTwo)
01752                 {
01753                         fftrows = NextPow2(fftrows);
01754                         fftcols = NextPow2(fftcols);
01755                 }
01756                 
01757                 Matrix<ComplexDouble> temp1;
01758                 if(b == ZeroPad || b == Symmetric || b == Replicate) 
01759                 {
01760                         temp1 = Matrix<ComplexDouble>(fftrows, fftcols,0);
01761                         for(int j=0; j<input.Columns(); j++)
01762                         {
01763                                 for(int i=0; i<input.Rows(); i++)
01764                                 {
01765                                         temp1.ElemNC(i,j) = input.ElemNC(i,j);
01766                                 }
01767                         }
01768                 }
01769                 else if(b == Circular)
01770                 {
01771                         temp1 = input;
01772                 }
01773                 else
01774                 {
01775                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01776                         Utility::RunTimeError("Invalid Border type!");
01777                 }
01778                 
01779                 if(b == Symmetric)
01780                 {
01781                         int wid1 = filter.Columns()/2;
01782                         int wid2 = filter.Columns() - wid1 - 1;
01783                         if(wid1 != 0)
01784                         {
01785                                 Matrix<ComplexDouble> c1 = FlipLRI( input.Slice(0, input.Rows()-1, 0, wid1-1) );
01786                                 temp1.ReadFromMatrix(c1, 0, temp1.Columns()-wid1);
01787                         }
01788                         if(wid2 != 0)
01789                         {
01790                                 Matrix<ComplexDouble> c2 = FlipLRI( input.Slice(0, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) );
01791                                 temp1.ReadFromMatrix(c2, 0, input.Columns());
01792                         }
01793                         
01794                         int hei1 = filter.Rows()/2;
01795                         int hei2 = filter.Rows() - hei1 - 1;
01796                         if(hei1 != 0)
01797                         {
01798                                 Matrix<ComplexDouble> r1 = FlipUDI( input.Slice(0, hei1-1, 0, input.Columns()-1) );
01799                                 temp1.ReadFromMatrix(r1, temp1.Rows()-hei1, 0);
01800                         }
01801                         if(hei2 != 0)
01802                         {
01803                                 Matrix<ComplexDouble> r2 = FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, input.Columns()-1) );
01804                                 temp1.ReadFromMatrix(r2, input.Rows(), 0);
01805                         }
01806                         
01807                         
01808                         if(hei1 != 0 && wid1 != 0)
01809                         {
01810                                 Matrix<ComplexDouble> q11 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, 0, wid1-1) ));
01811                                 temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
01812                         }
01813 
01814                         if(hei1 != 0 && wid2 != 0)
01815                         {
01816                                 Matrix<ComplexDouble> q12 = FlipLRI(FlipUDI( input.Slice(0, hei1-1, input.Columns()-wid2, input.Columns()-1) ));
01817                                 temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
01818                         }
01819 
01820                         if(hei2 != 0 && wid1 != 0)
01821                         {
01822                                 Matrix<ComplexDouble> q21 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, 0, wid1-1) ));
01823                                 temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
01824                         }
01825 
01826                         if(hei2 != 0 && wid2 != 0)
01827                         {
01828                                 Matrix<ComplexDouble> q22 = FlipLRI(FlipUDI( input.Slice(input.Rows()-hei2, input.Rows()-1, input.Columns()-wid2, input.Columns()-1) ));
01829                                 temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
01830                         }
01831                 }
01832                 else if(b == Replicate)
01833                 {
01834                         int wid1 = filter.Columns()/2;
01835                         int wid2 = filter.Columns() - wid1 - 1;
01836                         if(wid1 != 0)
01837                         {
01838                                 Matrix<ComplexDouble> c1 = input.Slice(0, input.Rows()-1, 0, 0);
01839                                 temp1.ReadFromMatrix(RepMat(c1,1,wid1), 0, temp1.Columns()-wid1);
01840                         }
01841                         if(wid2 != 0)
01842                         {
01843                                 Matrix<ComplexDouble> c2 = input.Slice(0, input.Rows()-1, input.Columns()-1, input.Columns()-1);
01844                                 temp1.ReadFromMatrix(RepMat(c2,1,wid2), 0, input.Columns());
01845                         }
01846                         
01847                         int hei1 = filter.Rows()/2;
01848                         int hei2 = filter.Rows() - hei1 - 1;
01849                         if(hei1 != 0)
01850                         {
01851                                 Matrix<ComplexDouble> r1 = input.Slice(0, 0, 0, input.Columns()-1);
01852                                 temp1.ReadFromMatrix(RepMat(r1,hei1,1), temp1.Rows()-hei1, 0);
01853                         }
01854                         if(hei2 != 0)
01855                         {
01856                                 Matrix<ComplexDouble> r2 = input.Slice(input.Rows()-1, input.Rows()-1, 0, input.Columns()-1);
01857                                 temp1.ReadFromMatrix(RepMat(r2,hei2,1), input.Rows(), 0);
01858                         }
01859                         
01860                         
01861                         if(hei1 != 0 && wid1 != 0)
01862                         {
01863                                 Matrix<ComplexDouble> q11( hei1, wid1, input.ElemNC(0,0) );
01864                                 temp1.ReadFromMatrix(q11, temp1.Rows()-hei1, temp1.Columns()-wid1);
01865                         }
01866 
01867                         if(hei1 != 0 && wid2 != 0)
01868                         {
01869                                 Matrix<ComplexDouble> q12( hei1, wid2, input.ElemNC(0,input.Columns()-1) );
01870                                 temp1.ReadFromMatrix(q12, temp1.Rows()-hei1, input.Columns());
01871                         }
01872 
01873                         if(hei2 != 0 && wid1 != 0)
01874                         {
01875                                 Matrix<ComplexDouble> q21( hei2, wid1, input.ElemNC(input.Rows()-1,0) );
01876                                 temp1.ReadFromMatrix(q21, input.Rows(), temp1.Columns()-wid1);
01877                         }
01878 
01879                         if(hei2 != 0 && wid2 != 0)
01880                         {
01881                                 Matrix<ComplexDouble> q22(hei2,wid2,input.ElemNC(input.Rows()-1,input.Columns()-1));
01882                                 temp1.ReadFromMatrix(q22, input.Rows(), input.Columns());
01883                         }
01884                 }
01885 
01886                 
01887                 // Allocate area for the filter
01888                 Matrix<ComplexDouble> temp2;
01889                 if(b == ZeroPad || b == Symmetric || b == Replicate)
01890                 {
01891                         temp2 = Matrix<ComplexDouble>(fftrows, fftcols,0);
01892                         for(int j=0; j<filter.Columns(); j++)
01893                         {
01894                                 for(int i=0; i<filter.Rows(); i++)
01895                                 {
01896                                         temp2.ElemNC(i,j) = filter.ElemNC(i,j);
01897                                 }
01898                         }
01899                 }
01900                 else if(b == Circular)
01901                 {
01902                         temp2 = Matrix<ComplexDouble>(input.Rows(), input.Columns(),0);
01903                         for(int j=0; j<filter.Columns(); j++)
01904                         {
01905                                 for(int i=0; i<filter.Rows(); i++)
01906                                 {
01907                                         temp2.ElemNC(i,j) = filter.ElemNC(i,j);
01908                                 }
01909                         }
01910                 }
01911                 else
01912                 {
01913                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01914                         Utility::RunTimeError("Invalid Border type!");
01915                 }
01916 
01917 
01918                 // FFT and inverse FFT
01919                 Matrix<ComplexDouble> o = IFFT2( FFT2(temp1) * FFT2(temp2) );
01920                 
01921                 // Based on the filter area, crop the border.
01922                 if(fa == Same)
01923                 {
01924                         o = o.Slice(filter.Rows()/2, filter.Rows()/2+input.Rows()-1, filter.Columns()/2, filter.Columns()/2+input.Columns()-1);
01925                 }
01926                 else if(fa == Valid)
01927                 {
01928                         o = o.Slice(filter.Rows()-1, filter.Rows()-1+input.Rows()-filter.Rows(), filter.Columns()-1, filter.Columns()-1+input.Columns()-filter.Columns());
01929                 }
01930                 else if(fa != Full)
01931                 {
01932                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01933                         Utility::RunTimeError("Invalid FilterArea type!");
01934                 }
01935                 
01936                 return o;               
01937         }
01938         
01939         Matrix<ComplexDouble> Conv2(Vector<ComplexDouble>& filter1, Vector<ComplexDouble>& filter2, Matrix<ComplexDouble>& input, FilterArea fa, Border b, bool PowerOfTwo)
01940         {
01941                 // Check if the filter sizes is smaller than the input size
01942                 if(filter1.Length() > input.Rows() || filter2.Length() > input.Columns())
01943                 {
01944                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
01945                         Utility::RunTimeError("Filter dimensions should not be larger than the input signal! filter1 <= # of rows, filter2 <= # of columns");
01946                 }
01947                 
01948                 Vector<ComplexDouble> *cols = new (std::nothrow) Vector<ComplexDouble>[input.Columns()];
01949                 Utility::CheckPointer(cols);
01950                 
01951                 for(int i=0; i<input.Columns(); i++)
01952                 {
01953                         cols[i] = Conv(input[i], filter1, fa, b, PowerOfTwo);
01954                 }
01955                 
01956                 int rlength = cols[0].Length();
01957                 
01958                 Vector<ComplexDouble> *rows = new (std::nothrow) Vector<ComplexDouble>[rlength];
01959                 Utility::CheckPointer(rows);
01960                 
01961                 for(int i=0; i<rlength; i++)
01962                 {
01963                         rows[i] = Vector<ComplexDouble>(input.Columns());
01964                         for(int j=0; j<input.Columns(); j++)
01965                         {
01966                                 rows[i][j] = cols[j][i];
01967                         }
01968                 }
01969                 
01970                 delete [] cols;
01971                 
01972                 for(int i=0; i<rlength; i++)
01973                 {
01974                         rows[i] = Conv(rows[i], filter2, fa, b, PowerOfTwo);
01975                 }
01976                 
01977                 int clength = rows[0].Length();
01978                 Matrix<ComplexDouble> temp(rlength, clength);
01979                 for(int i=0; i<rlength; i++)
01980                 {
01981                         for(int j=0; j<clength; j++)
01982                         {
01983                                 temp.ElemNC(i,j) = rows[i][j];
01984                         }
01985                 }
01986                 
01987                 delete [] rows;
01988                 
01989                 return temp;
01990         }
01991 
01992 
01993 
01994 
01995 
01996 };
01997 
01998 
01999 
02000 
02001