[C语⾔]矩阵运算
最近要做⼀个MFC的上位机，⽤到CSP滤波算法，这玩意⼉在MATLAB ⾥相当简单就能实现但C⾥⾯实现起来太蛋疼，写了⼀个晚上才把这个算法⽤到的矩阵运算部分的函数写的差不多，为了避免以后再重复造轮⼦，现在这⾥写⼀下备份⼀下吧。。
1.矩阵乘法
//矩阵乘法
/********参数表*******
@Parameter x: m⾏k列矩阵（⽤⼀维数组表⽰）
@Parameter y: k⾏n列矩阵（⽤⼀维数组表⽰）
@Parameter m,k,n: 矩阵⾏列参数
@Parameter z: m⾏n列输出矩阵（⽤⼀维数组表⽰）
***********************/
void MulMatrixDD(double *x,double *y, int m,int k,int n, double *z)
{
for(int nm=0; nm<m; nm++)
for(int nn=0; nn<n; nn++)
for(int nk=0; nk<k; nk++)
z[nm*n+nn] += x[nm*k+nk]*y[nk*n+nn];
}
因为输⼊的x,y可能是不同的类型（short, double）所以我只简单的复制了⼏个函数来表⽰，⽐如输⼊如果是x如果是short型，y是double型这个函数就不好⽤了。如果有什么更好的通⽤⽅法可以跟我交流下
2.⽅阵转置
//⽅阵转置
/********参数表*******
@Parameter x: m⾏m列矩阵（⽤⼀维数组表⽰）
@Parameter m: 矩阵⾏列数
***********************/
void TransSquareD(double *x, int m)
{
double temp;
for(int nm=0; nm<m; nm++){ //对原矩阵第nm⾏
for(int nn=0; nn<nm; nn++){ //对原矩阵第nn列
temp = x[nm*m+nn]; //z矩阵第nn⾏第nm列
x[nm*m+nn] = x[nn*m+nm];
x[nn*m+nm] = temp;}}
}
3.⾮⽅阵转置
//⾮⽅阵转置
/********参数表*******
@Parameter x: m⾏n列矩阵（⽤⼀维数组表⽰）
@Parameter m,n: 矩阵⾏列数
@Parameter z: n⾏m列矩阵（⽤⼀维数组表⽰）
***********************/
void TransMatrixD(double *x, int m, int n, double *z)
{
for(int nm=0; nm<m; nm++) //对原矩阵第nm⾏
for(int nn=0; nn<n; nn++) //对原矩阵第nn列
z[nn*m+nm] = x[nm*n+nn]; //z矩阵第nn⾏第nm列
}
void TransMatrixS(short *x, int m, int n, double *z)
{
for(int nm=0; nm<m; nm++) //对原矩阵第nm⾏
for(int nn=0; nn<n; nn++) //对原矩阵第nn列
z[nn*m+nm] = (double)x[nm*n+nn]; //z矩阵第nn⾏第nm列
c语言二维数组表示方法}
4.协⽅差矩阵
/
/协⽅差矩阵函数
/********参数表*******
@Parameter X[m_cov][n_cov]: m_cov⾏n_cov列矩阵（⽤⼆维数组表⽰）
***********************/
void CovMat(short X[m_cov][n_cov])
{
for(int i=0; i<m_cov; i++)
for(int j=0; j<m_cov; j++)
CovMatrix[i][j] = cov(*(X+i),*(X+j));
}
//协⽅差函数
double cov(short *x, short *y)
{
double sumx = 0;
double sumy = 0;
double sum = 0;
int kk = 0;
for(kk = 0; kk<n_cov; kk++)
{
sumx += x[kk];
sumy += y[kk];
}
sumx /= n_cov;
sumy /= n_cov;
for(kk = 0; kk<n_cov; kk++)
{
sum += (x[kk]-sumx)*(y[kk]-sumy);
}
sum /= (n_cov-1);
return sum;
}
5.求矩阵特征值和特征向量
这个短时间⾃⼰造轮⼦太⿇烦了，，参考csdn上⼀篇博客：blog.csdn/zhouxuguang236/article/details/40212143 改了⼀下/********参数表*******
* @brief 求实对称矩阵的特征值及特征向量的雅克⽐法
* 利⽤雅格⽐(Jacobi)⽅法求实对称矩阵的全部特征值及特征向量
@Parameter pMatrix 长度为n*n的数组，存放实对称矩阵
@Parameter nDim 矩阵的阶数
@Parameter pdblVects 长度为n*n的数组，返回特征向量(按列存储)
@Parameter dbEps 精度要求
@Parameter nJt 整型变量，控制最⼤迭代次数
@Parameter pdbEigenValues 特征值数组
***********************/
void JacbiCor(double * pMatrix,int nDim, double *pdblVects, double *pdbEigenValues, double dbEps,int nJt)
{
int i,j;
for(i = 0; i < nDim; i ++)
{
pdblVects[i*nDim+i] = 1.0f;
for(int j = 0; j < nDim; j ++)
{
if(i != j)
pdblVects[i*nDim+j]=0.0f;
}
}
int nCount = 0; //迭代次数
while(1)
{
//在pMatrix的⾮对⾓线上到最⼤元素
double dbMax = pMatrix[1];
int nRow = 0;
int nCol = 1;
for (i = 0; i < nDim; i ++) //⾏
{
for (j = 0; j < nDim; j ++) //列
{
double d = fabs(pMatrix[i*nDim+j]);
if((i!=j) && (d> dbMax))
{
dbMax = d;
nRow = i;
nCol = j;
}
}
}
if(dbMax < dbEps) //精度符合要求
break;
if(nCount > nJt) //迭代次数超过限制
break;
nCount++;
double dbApp = pMatrix[nRow*nDim+nRow];
double dbApq = pMatrix[nRow*nDim+nCol];
double dbAqq = pMatrix[nCol*nDim+nCol];
//计算旋转⾓度
double dbAngle = 0.5*atan2(-2*dbApq,dbAqq-dbApp);
double dbSinTheta = sin(dbAngle);
double dbCosTheta = cos(dbAngle);
double dbSin2Theta = sin(2*dbAngle);
double dbCos2Theta = cos(2*dbAngle);
pMatrix[nRow*nDim+nRow] = dbApp*dbCosTheta*dbCosTheta +
dbAqq*dbSinTheta*dbSinTheta + 2*dbApq*dbCosTheta*dbSinTheta;
pMatrix[nCol*nDim+nCol] = dbApp*dbSinTheta*dbSinTheta +
dbAqq*dbCosTheta*dbCosTheta - 2*dbApq*dbCosTheta*dbSinTheta;
pMatrix[nRow*nDim+nCol] = 0.5*(dbAqq-dbApp)*dbSin2Theta + dbApq*dbCos2Theta; pMatrix[nCol*nDim+nRow] = pMatrix[nRow*nDim+nCol];
for(i = 0; i < nDim; i ++)
{
if((i!=nCol) && (i!=nRow))
{
int u = i*nDim + nRow; //p
int w = i*nDim + nCol; //q
dbMax = pMatrix[u];
pMatrix[u]= pMatrix[w]*dbSinTheta + dbMax*dbCosTheta;
pMatrix[w]= pMatrix[w]*dbCosTheta - dbMax*dbSinTheta;
}
}
for (j = 0; j < nDim; j ++)
{
if((j!=nCol) && (j!=nRow))
{
int u = nRow*nDim + j; //p
int w = nCol*nDim + j; //q
dbMax = pMatrix[u];
pMatrix[u]= pMatrix[w]*dbSinTheta + dbMax*dbCosTheta;
pMatrix[w]= pMatrix[w]*dbCosTheta - dbMax*dbSinTheta;
}
}
//计算特征向量
for(i = 0; i < nDim; i ++)
{
int u = i*nDim + nRow; //p
int w = i*nDim + nCol; //q
dbMax = pdblVects[u];
pdblVects[u] = pdblVects[w]*dbSinTheta + dbMax*dbCosTheta;
pdblVects[w] = pdblVects[w]*dbCosTheta - dbMax*dbSinTheta;
}
}
for(i = 0; i < nDim; i ++)
{
pdbEigenValues[i] = pMatrix[i*nDim+i];
}
//设定正负号
for(i = 0; i < nDim; i ++)
{
double dSumVec = 0;
for(j = 0; j < nDim; j ++)
dSumVec += pdblVects[j * nDim + i];
if(dSumVec<0)
{
for(j = 0;j < nDim; j ++)
pdblVects[j * nDim + i] *= -1;
}
}
}
⼆维数组转⼀维数组：
const short m_cov = 3;
const short n_cov = 3;
short X[m_cov][n_cov] = {0};
short *XFlat;
/************X测试************/
X[0][0] = 7;X[0][1] = 1;X[0][2] = 8;
X[1][0] = 4;X[1][1] = 5;X[1][2] = 8;
X[2][0] = 10;X[2][1] = 4;X[2][2] = 2; /*****************************/
XFlat = (short *)X;
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论