/*------------------------------------------------------------------------------* * File Name:matrix.c * * Creation:July 22th, 2002 * * Purpose: OriginC Source C file * * Copyright (c) OriginLab Corp.2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 * * All Rights Reserved * * * * Modification Log: * * TCZ 07/29/02 QA70-2497 ADD_GLOBAL_STATS_SUMMARY_FUNCTION * * ER 01/29/03 QA70_3802 ADD_VECTOR_FFT_IFFT * CPY 7/7/03 changed to include only origin.h to take advantage of PCH *------------------------------------------------------------------------------*/ #include // main Origin C header that is precompiled and already include most headers // this file include most of the other header files except the NAG header, which takes longer to compile // NAG routines #include // this contains the FFT NAG headers, #include //this contain the linear algebra function #include //this file contain the simple statistics function //////////////////////////////////////////////////////////////////////////////////// #include // start your functions here //local functions for the internal use // This function will do the 2d fft with the data form required by nag funtion. static int fft2_for_complex_matrix(int nNumRows, int nNumCols, matrix& matReal, matrix& matImage, matrix& matFFT2Result); // This function will take a matrix and set the result matrix accoring to the parameter values. // in other words, it will do the neccessory padding with zero or truncting. static int set_matrix_with_padding_truncting(matrix & matSource, matrix & matResult, int nRowProcess , int nColProcess); // This function will take care of the complex data matrix and sperate it into two double matrix. static int seperate_complex_mat_to_real_and_imag_mat(matrix& matSource, matrix& matReal, matrix& matImage, int nRowProcess, int nColProcess); /** int FFT2(matrix& matSource, matrix& matFFT2Result, int nRowProcess = -1, int nColProcess = -1); int FFT2(matrix& matSource, matrix& matFFT2Result, int nRowProcess = -1, int nColProcess = -1); These two Functions will do the 2 dimension Fast Fourier Transform. Remarks: Right now, we only support the data matrix in double and complex form. If user have the data matrix in integer form, user should first transform it to double. It is a member function of matrix. Also, if user has a vector data, it should firstly be transformed into a matrix and then do the FFT. Example: Parameters: matSource: The data source matrix matFFT2Result: The matrix containing the 2 dimension FFT results nRowProcess: nColProcess: The number of rows or columns user provided. If they are -1, there is no padding or truncating. If they are more than the sizes of the matrix, the function will pad zero to meet the size If less than, it will truncate the data matrix. Return: If succeed, it will return 0; Otherwise it returns error indications. Error: -1: Provided column number or row number. Must be -1 or positive integer -2: The data source matrix is empty -3: When calling member function, it fails. -4: When calling GetReal or GetImage member function of matrix class, it failed. Nag Error: 11: The column or the row size of matrix is less then 1 73: System failed to allocate memory */ int FFT2(matrix& matSource, matrix& matFFT2Result, int nRowProcess, int nColProcess) { matrix matReal; matrix matTemp(matSource); int nRet; if(nRet = set_matrix_with_padding_truncting(matTemp, matReal, nRowProcess, nColProcess)) return nRet; matrix matImage; int nNumCols = matReal.GetNumCols(); int nNumRows = matReal.GetNumRows(); matImage.SetSize(nNumRows, nNumCols); matImage = 0; if(nRet = fft2_for_complex_matrix(nNumRows, nNumCols, matReal, matImage, matFFT2Result)) return nRet; return nRet; } static int fft2_for_complex_matrix(int nNumRows, int nNumCols, matrix& matReal, matrix& matImage, matrix& matFFT2Result) { vector vecTrigRows; vecTrigRows.SetSize(2 * nNumRows); //cal nage function to prepare the trig value int nRet; if( nRet = nag_fft_init_trig(nNumRows, vecTrigRows)) return nRet; vector vecTrigCols; vecTrigCols.SetSize(2 * nNumCols); if( nRet = nag_fft_init_trig(nNumCols, vecTrigCols)) return nRet; //call the nag function to do the fft2 if(nRet = nag_fft_2d_complex(nNumRows, nNumCols, matReal, matImage, vecTrigRows, vecTrigCols)) return nRet; //set the data back if( nRet = matFFT2Result.MakeComplex(matReal, matImage)) return nRet; } int FFT2(matrix& matSource, matrix& matFFT2Result, int nRowProcess, int nColProcess ) { int nRet; matrix matReal, matImage; if(nRet = seperate_complex_mat_to_real_and_imag_mat(matSource, matReal, matImage, nRowProcess, nColProcess)) return nRet; int nNumCols = matImage.GetNumCols(); int nNumRows = matImage.GetNumRows(); if(nRet = fft2_for_complex_matrix(nNumRows, nNumCols, matReal, matImage, matFFT2Result)) return nRet; } static int set_matrix_with_padding_truncting(matrix & matSource, matrix & matResult, int nRowProcess , int nColProcess) { if(nRowProcess < -1 || nRowProcess == 0 ||nColProcess < -1 ||nColProcess == 0) return FFT_ERROR_COULUMN_ROW_NUMBER_NEGATIVE; int nNumCols = matSource.GetNumCols(); int nNumRows = matSource.GetNumRows(); if(nNumCols == 0 || nNumRows == 0) return FFT_ERROR_SOURCE_MATRIX_EMPTY; int nIntendNumCols = nNumCols; int nIntendNumRows = nNumRows; if(nRowProcess != -1) nIntendNumRows = nRowProcess; if(nColProcess != -1) nIntendNumCols = nColProcess; int nRet; matResult = matSource; if(nIntendNumRows > nNumRows)//so we should increase the dimension and do the padding { matrix matTemp; matTemp.SetSize(nIntendNumRows - nNumRows, nNumCols ); matTemp = 0; if(nRet = matSource.Concatenate(1,matTemp, matResult)) return nRet; } BOOL bOK = TRUE; if (nIntendNumRows < nNumRows) bOK = matSource.GetSubMatrix(matResult, 0, -1, 0, nIntendNumRows - 1); if(!bOK) return FFT_ERROR_GETSUBMATRIX_FUNCTION_ERROR; if(nIntendNumCols > nNumCols) { matrix matTemp; matTemp.SetSize(nIntendNumRows, nIntendNumCols - nNumCols ); matTemp = 0; if(nRet = matResult.Concatenate(2, matTemp, matResult)) return nRet; } if(nIntendNumCols < nNumCols) { matrix matResultTemp(matResult); bOK = matResultTemp.GetSubMatrix(matResult, 0, nIntendNumCols -1, 0, -1); if(!bOK) return FFT_ERROR_GETSUBMATRIX_FUNCTION_ERROR; } return 0; } /* int IFFT2(matrix& matSource, matrix& matIFFT2Result, int nRowProcess = -1, int nColProcess = -1); Remarks: Right now, we only support the data matrix in complex form. If user have the data matrix in double form, user should add the imaginary part as zero to make a complex matrix. Matrix class has a member function to make the complex matrix. Also, if user has a vector data, it should firstly be transformed into a matrix and then one can do the inverse FFT. Example: Parameters: matSource: The data source matrix matIFFT2Result: The matrix containing the 2 dimension inverse FFT results nRowProcess: nColProcess:The number of rows or columns user provided. if they are -1, there are no padding or truncating. if they are larger than the size of the matrix, the function will pad zero to meet the size. If less than, it will truncate the data matrix. Return: If succeed, it will return 0; Otherwise it returns error indications. Error: -1: Provided column number or row number. Must be -1 or positive integer -2: The data source matrix is empty -3: When calling member function, it fails. -4: When calling GetReal or GetImage member function of matrix class, it failed. Nag Error: 11: The column or the row size of matrix is less then 1 73: System failed to allocate memory */ int IFFT2(matrix& matSource, matrix& matFFT2Result, int nRowProcess, int nColProcess) { int nRet; matrix matSourceCopy(matSource); if(nRet = matSourceCopy.Conjugate()) return nRet; if(nRet = FFT2(matSourceCopy, matFFT2Result,nRowProcess, nColProcess )) return nRet; if(nRet = matFFT2Result.Conjugate()) return nRet; } /* int FFT(matrix& matSource, matrix& matFFTResult, int nColProcess = -1); int FFT(matrix& matSource, matrix& matFFTResult, int nColProcess = -1); These two Functions will do the 1 dimension Fast Fourier Transform. Remarks: Right now, we only support the data matrix in double and complex form. If user has the data matrix in integer form, user should firstly transform it to double. It is a member function of matrix. Also, if user has a vector data, he should firstly transform data into a matrix and then one can do the FFT. Unlike the matlab, this function will do the FFT on each row (Matlab do FFT on each column). Example: Parameters: matSource: The data source matrix matFFTResult: The matrix containing the 1 dimension FFT results nColProcess: The number of columns user provided. If it is -1, function will not do any padding or truncating. But if it is larger than the size of the matrix, it will pad zero to meet the size if less than, it will truncate the data matrix. Return: If succeed, it will return 0; Otherwise it returns error indications. Error: -1: Provided column number or row number. Must be -1 or positive integer -2: The data source matrix is empty -3: When calling member function, it fails. -4: When calling GetReal or GetImage member function of matrix class, it failed. Nag Error: 11: The column or the row size of matrix is less then 1 73: System failed to allocate memory */ int FFT(matrix& matSource, matrix& matFFTResult, int nColProcess) { matrix matTemp; int nRet; if(nRet = set_matrix_with_padding_truncting(matSource, matTemp, -1, nColProcess)) return nRet; int nNumCols = matTemp.GetNumCols(); int nNumRows = matTemp.GetNumRows(); vector vecTrig; vecTrig.SetSize(2 * nNumCols); if(nRet = nag_fft_init_trig(nNumCols, vecTrig)) return nRet; if(nRet = nag_fft_multiple_real(nNumRows, nNumCols, matTemp, vecTrig)) return nRet; matrix matReal, matImage; matReal.SetSize(nNumRows, nNumCols); matImage.SetSize(nNumRows, nNumCols); if(nRet = nag_multiple_hermitian_to_complex(nNumRows, nNumCols, matTemp, matReal, matImage)) return nRet; if(nRet = matFFTResult.MakeComplex(matReal, matImage)) return nRet; return nRet; } int FFT(matrix& matSource, matrix& matFFTResult, int nColProcess) { int nRet; matrix matReal, matImage; if(nRet = seperate_complex_mat_to_real_and_imag_mat(matSource, matReal, matImage, -1, nColProcess)) return nRet; int nNumCols = matImage.GetNumCols(); int nNumRows = matImage.GetNumRows(); vector vecTrig; vecTrig.SetSize(2 * nNumCols); if(nRet = nag_fft_init_trig(nNumCols, vecTrig)) return nRet; if(nRet = nag_fft_multiple_complex(nNumRows, nNumCols, matReal, matImage, vecTrig )) return nRet; if(nRet = matFFTResult.MakeComplex(matReal, matImage)) return nRet; return nRet; } /* int IFFT(matrix& matSource, matrix& matFFTResult, int nColProcess = -1); This Function will do the 1 dimension inverse Fast Fourier Transform. Remarks: Right now, we only support the data matrix in complex form. If user has the data matrix in double form, user should add the imaginary part as zero to make a complex matrix. Matrix class has a member function to make the complex matrix. also, if user has a vector data, it should be firstly transformed into a matrix and then one can do the inverse FFT. Example: Parameters: matSource:The data source matrix matIFFT2Result: The matrix containing the 1 dimension inverse FFT results nColProcess:The number of columns user provided. if it is -1, there is no padding or truncting. if it is more than the size of the matrix, it will pad zero to meet the size if less than, it will truncate the data matrix. Return: If succeed, it will return 0; Otherwise it returns error indications. Error: -1: Provided column number or row number must be -1 or positve integer -2: The data source matrix is empty -3: When calling member function, it fails. -4: When calling GetReal or GetImage member function of matrix class, it failed. Nag Error: 11: The coloumn or the row size of matrix is less then 1; 73: System failed to allocate memory */ int IFFT(matrix& matSource, matrix& matFFTResult, int nColProcess) { int nRet; matrix matSourceCopy(matSource); if(nRet = matSourceCopy.Conjugate()) return nRet; if(nRet = FFT(matSourceCopy, matFFTResult, nColProcess)) return nRet; if(nRet = matFFTResult.Conjugate()) return nRet; } static int seperate_complex_mat_to_real_and_imag_mat(matrix& matSource, matrix& matReal, matrix& matImage, int nRowProcess, int nColProcess) { matrix matTemp; BOOL bOK; bOK = matSource.GetReal(matTemp); if(!bOK) return FFT_ERROR_GETREAL_FUNCTION_FAIL; int nRet; if(nRet = set_matrix_with_padding_truncting(matTemp, matReal, nRowProcess, nColProcess)) return nRet; bOK = matSource.GetImaginary(matTemp); if(!bOK) return FFT_ERROR_GETREAL_FUNCTION_FAIL; if(nRet = set_matrix_with_padding_truncting(matTemp, matImage, nRowProcess, nColProcess)) return nRet; } // ER 01/29/03 QA70_3802 ADD_VECTOR_FFT_IFFT // The following function performs 1-D FFT on a complex vector. // It uses the matrix-based 1-D FFT function. int FFT(vector& vecSource, vector& vecFFTResult) { // Get size of source vector int iSize = vecSource.GetSize(); if(iSize < 1) return -5; // Create a complex matrix of size 1xn and copy source vector into matrix matrix matSource; matSource.SetSize(1,iSize); matSource.SetByVector(vecSource); // Create complex matrix to hold result of FFT operation matrix matFFTResult; // Call the matrix-based FFT function to perform FFT on this one row int iRet = FFT(matSource, matFFTResult); if(iRet) return iRet; // Extract FFT result from result matrix, into a complex matrix, and return matFFTResult.GetAsVector(vecFFTResult); // return 0 on success return 0; } // ER 01/29/03 QA70_3802 ADD_VECTOR_FFT_IFFT // The following function performs 1-D inverse FFT on a complex vector. // It uses the matrix-based 1-D inverse FFT function. int IFFT(vector& vecSource, vector& vecFFTResult) { // Get size of source vector int iSize = vecSource.GetSize(); if(iSize < 1) return -5; // Create a complex matrix of size 1xn and copy source vector into matrix matrix matSource; matSource.SetSize(1,iSize); matSource.SetByVector(vecSource); // Create complex matrix to hold result of IFFT operation matrix matFFTResult; // Call the matrix-based IFFT function to perform IFFT on this one row int iRet = IFFT(matSource, matFFTResult); if(iRet) return iRet; // Extract IFFT result from result matrix, into a complex matrix, and return matFFTResult.GetAsVector(vecFFTResult); // return 0 on success return 0; } /** int SVD(matrix& matSource, matrix& matS, matrix& matU, matrix& matV); int SVD(matrix& matSource, matrix& matS); int SVD(matrix& matSource, matrix& matS); int SVD(matrix& matSource, matrix& matS, matrix& matU, matrix& matV); These four functions are used to do the singular value decomposition. matSource = matU* matS* matV' Example: Parameters: matSource: The m * n source data matrix matS: The min(m, n) diaganol sigular matrix matU: The left side m * m unitary matrix matV: The right side n* n unitary matrix return: if succeed, it will return 0 */ int SVD(matrix& matSource, matrix & matS, matrix& matU, matrix& matV) { //fist check the size of the matrix int nNumCols = matSource.GetNumCols(); int nNumRows = matSource.GetNumRows(); //no need to check size, because nag will take care for it //set the result matrix size matU.SetSize(nNumRows, nNumRows); matV.SetSize(nNumCols, nNumCols); vector vecE; int nMinSize = nNumCols > nNumRows ? nNumRows : nNumCols; vecE.SetSize(nMinSize -1); vector vecS; vecS.SetSize(nMinSize); //call the nag function int nRet; int nNumIter, nfailIfor; //we make a copy of source data here becasue we do not want to overwirte the data matrix matSourceCopy(matSource); if(nRet = nag_real_svd(nNumRows, nNumCols, matSourceCopy, nNumCols, 0, NULL, 0, TRUE, matU, nNumRows, vecS, TRUE, matV, nNumCols, &nNumIter, vecE, &nfailIfor)) return nRet; //set back the U V matrix if(nNumRows >= nNumCols) matU = matSourceCopy; else matV = matSourceCopy; matV.Transpose(); // Nag return MatV' //set back matS nRet = matS.SetDiagonal(vecS, 0); //here we do not consider nfailinfo, becasue it is rarely occurs. return nRet; } int SVD(matrix& matSource, matrix& matS, matrix& matU, matrix& matV) { //fist check the size of the matrix int nNumCols = matSource.GetNumCols(); int nNumRows = matSource.GetNumRows(); //no need to check size, because nag will take care for it //set the result matrix size matU.SetSize(nNumRows, nNumRows); matV.SetSize(nNumCols, nNumCols); vector vecE; int nMinSize = nNumCols > nNumRows ? nNumRows : nNumCols; vecE.SetSize(nMinSize -1 ); vector vecS; vecS.SetSize(nMinSize); //call the nag function int nRet; int nNumIter, nfailIfor; //we make a copy of source data here becasue we do not want to overwirte the data matrix matSourceCopy(matSource); if(nRet = nag_complex_svd(nNumRows, nNumCols, matSourceCopy, nNumCols, 0, NULL, 0, TRUE, matU, nNumRows, vecS, TRUE, matV, nNumCols, &nNumIter, vecE, &nfailIfor)) return nRet; //set back the U V matrix if(nNumRows >= nNumCols) matU = matSourceCopy; else matV = matSourceCopy; matV.Transpose(); // Nag return MatV' //set back matS nRet = matS.SetDiagonal(vecS, 0); //here we do not consider nfailinfo, becasue it is rarely occurs. return nRet; } int SVD(matrix& matSource, matrix& matS) { matrix matU, matV; int nRet = 0; nRet = SVD(matSource, matS, matU, matV); return nRet; } int SVD(matrix& matSource, matrix& matS) { matrix matU, matV; int nRet = 0; nRet = SVD(matSource, matS, matU, matV); return nRet; } /** int Trace(matrix & matSource, int & nSum); int Trace(matrix & matSource, int & nSum); int Trace(matrix & matSource, double & dSum); int Trace(matrix & matSource, complex & cSum); These four functions will give the trace of a matrix Remarks: we do not support every matrix type. Example: Parameters: matSource: The source matrix nSum: dSum: cSum: These are the variables to store the trace. Return: if succeed, it will return 0 */ #define CHECK_MATRIX_FOR_TRACE {\ int nNumCols = matSource.GetNumCols(); \ int nNumRows = matSource.GetNumRows(); \ if(!nNumCols || !nNumRows) \ return FFT_ERROR_SOURCE_MATRIX_EMPTY; \ } int Trace(matrix &matSource, int& nSum) { CHECK_MATRIX_FOR_TRACE; vector vecDiag; int nRet; if(nRet = matSource.GetDiagonal(vecDiag)) return nRet; if(nRet = vecDiag.Sum(nSum)) return nRet; } int Trace(matrix &matSource, int& nSum) { CHECK_MATRIX_FOR_TRACE; vector vecDiag; int nRet; if(nRet = matSource.GetDiagonal(vecDiag)) return nRet; if(nRet = vecDiag.Sum(nSum)) return nRet; } int Trace(matrix &matSource, double& dSum) { CHECK_MATRIX_FOR_TRACE; vector vecDiag; int nRet; if(nRet = matSource.GetDiagonal(vecDiag)) return nRet; if(nRet = vecDiag.Sum(dSum)) return nRet; } int Trace(matrix & matSource, complex& cSum) { CHECK_MATRIX_FOR_TRACE; vector vecDiag; int nRet; if(nRet = matSource.GetDiagonal(vecDiag)) return nRet; if(nRet = vecDiag.Sum(cSum)) return nRet; } ///TCZ 07/29/02 QA70-2497 ADD_GLOBAL_STATS_SUMMARY_FUNCTION //#define MAT_ERR_WEIGHTS_DATA_SIZE_NOT_MATCH (-5) //defined in the matrix.h /** This function will give the summary of the statistics of a source data contained in a matrix Remarks: If user does not have the weight, then just supply a matix without any initialization for weight Example: Parameters: matData: The Data Source MatrixStatsSummary: The Results matWeights: The weight if any Return: Return 0 if succeed. MAT_ERR_WEIGHTS_DATA_SIZE_NOT_MATCH: The weight matrix size does not match that of the source int MatrixBasicStats(matrix & mData, MatrixStats & sMatStatsSum, matrix * mpWeights = NULL); */ int MatrixBasicStats(matrix & mData, MatrixStats & sMatStatsSum, matrix * mpWeights ) { //first check the size of the weight and data size int nNumCols = mData.GetNumCols(); int nNumRows = mData.GetNumRows(); if(mpWeights != NULL) { int nNumColsWeights = mpWeights->GetNumCols(); int nNumRowsWeights = mpWeights->GetNumRows(); if(nNumCols != nNumColsWeights || nNumRows != nNumRowsWeights) return MAT_ERR_WEIGHTS_DATA_SIZE_NOT_MATCH; } int nNumPoints = nNumCols * nNumRows; int nValid; double dMean, dSD, dMax, dMin, dWSum; int nRet; if(mpWeights == NULL) { if( nRet = nag_summary_stats_1var(nNumPoints, mData, NULL, &nValid, &dMean, &dSD, NULL, NULL, &dMin, &dMax, &dWSum)) return nRet; } else { if( nRet = nag_summary_stats_1var(nNumPoints, mData, *mpWeights, &nValid, &dMean, &dSD, NULL, NULL, &dMin, &dMax, &dWSum)) return nRet; } //call the second nag function to get the median and hinge double dRes[5]; if( nRet = nag_5pt_summary_stats(nNumPoints, mData, dRes)) return nRet; //set the structure value sMatStatsSum.nMissingValue = nNumPoints - nValid; sMatStatsSum.dMean = dMean; sMatStatsSum.dSD = dSD; sMatStatsSum.dMax = dMax; sMatStatsSum.dMin = dMin; sMatStatsSum.dMedian = dRes[2]; sMatStatsSum.dLowHinge = dRes[1]; sMatStatsSum.dUpHinge = dRes[3]; sMatStatsSum.dWSum = dWSum; return nRet; } //END-------------------------------ADD_GLOBAL_STATS_SUMMARY_FUNCTION