阵列信号处理中若干算法的C++实现之基础函数篇
来源:互联网 发布:淘宝司法拍卖房产 编辑:程序博客网 时间:2024/06/06 04:37
引言
在阵列信号处理中,多通道信号处理大多基于矩阵方法,本文主要介绍一些常见算法的
矩阵存储
二维矩阵在文中按照线性存储,对于
希尔伯特变换
在信号处理中,计算和信息传递方面的优势,复数的使用不可避免。在输入信号时,原始信号为实数,希尔伯特变换可以实现将实数转换为复数的功能。希尔伯特变换的频域计算公式如下
希尔伯特变换的物理意义十分明确,其将信号的正频率部分相移
1、做
2、由于
3、
利用希尔伯特变换实现信号从实数域到复数域的转换,算法实现如下:
1.实数信号作为实部传递给复数,复数的虚部置
2.将复信号做
3.
4.实部虚部交换位置得到新复信号。
5.对新的复信号做
//*****************************************************/**功能:* 利用希尔伯特变换将单通道信号从实数域变换到复数域* *参数:* 参数1---输入,信号长度* 参数2---输入,单通道实信号* 参数3---输出,恢复后的单通道复数信号* *返回值:* 无* *附加说明:* 程序运行需FFTW库,请先行配置好。*///******************************************************void matreat::heribtToComplex(int n,double* indata,complex<double>* outdata){ ... for(int i = 0;i<length;i++) { in[i][0] = indata[i];//输入实信号作为实部 in[i][1] = 0;//虚部置0 } p = fftw_plan_dft_1d(length,in,out,FFTW_FORWARD,FFTW_ESTIMATE);//傅里叶正变换 ... //正频率部分 for(int i = 0;i<midle;i++) { //实虚部互换 in[i][0] = out[i][1]; in[i][1] = -out[i][0];//实部取反 } //负频率部分 for(int i = midle;i<length;i++) { //实虚部互换 in[i][0] = -out[i][1];//虚部取反 in[i][1] = out[i][0]; } p = fftw_plan_dft_1d(length,in,out,FFTW_BACKWARD,FFTW_ESTIMATE);//傅里叶反变换 ...}
这里仅展示了单通道信号的处理方法,对于多通道信号方法类似,依次对每个通道信号进行希尔伯特变换即可。
Matlab 数据读/写
在save 'dataA.txt' -ascii -double dataA;
语句将变量
类似的load dataA.txt;
语句将
对于多维数据和复数数据,本文并未给出存取方案,在实际运用中可分为各部分依次存储即可。 下文的方案中,在
//*****************************************************/**功能:* 将数据保存成matlab可读取的.txt格式文件* *参数:* 参数1---输入,待保存数据* 参数2---输入,数据的总长度* 参数3---输出,文档保存位置* *返回值:* 无* *附加说明:* *///******************************************************void saveForMat(double indata[],int N,const string &dir){ ofstream datafile; datafile.open(dir);//若文件不存在会自动创建文件 for(int i=0;i<N;i++) { datafile<<indata[i]; datafile<<" "; } datafile.close();}
//*****************************************************/**功能:* 将读取matlab导出.txt格式文件* *参数:* 参数1---输入,待读取文档的位置* 参数2---输出,数据变量* *返回值:* 无* *附加说明:* *///******************************************************void readMatalbtext(const string &str,vector<double> &wdata){ ifstream file; file.open(str); if(!file) { cout<<"错误!!!!,请检查文本输入输出路径是否正确"<<endl; } else { double d; while (file >> d) { wdata.push_back(d);//将数据压入堆栈 } file.close();//关闭文件// }}
协方差矩阵计算
在阵列信号处理中,计算协方差如下式
原始信号输入数据数据大小为
举例,协方差矩阵数据如下:
利用该性质,计算时将矩阵分为上三角区域和下三角区域,先计算上三角区域数据,下三角区域对应位置的数据直接共轭赋值。
//*****************************************************/**功能:* 计算输入数据的协方差矩阵* *参数:* 参数1---输入,信号的通道数* 参数2---输入,每个通道数据长度* 参数3---输入,输入多通道信号* 参数3---输出,处理后矩阵* *返回值:* 无* *附加说明:***///******************************************************void getCovComplexMat(int m,int len,complex<double>*indata,complex<double>*outdata){ ... getConjComplexMat(m,len,indata,conjdata);//计算矩阵indata的共轭矩阵,详见完整代码 ... //计算上三角阵 for(int i=0;i<m;i++) { for(int j=i;j<m;j++) { indexA = i*len; indexB = j*len; ... while(num < len) { temp += indata[indexA]*conjdata[indexB]; indexA++; indexB++; num++; } ... } ... } //由对称性,下三角直接赋值 for(int i=1;i<m;i++) { for(int j=0;j<i;j++) { outdata[i*m+j].real(outdata[j*m+i].real()); outdata[i*m+j].imag(-outdata[j*m+i].imag()); } } //释放内存 ...}
方阵Toeplitz 化处理
均匀线阵的协方差矩阵为
对于
1、不同对角线,上三角区域,对角线起始元素的索引号相差
2、不同对角线,下三角区域,对角线起始元素的索引号相差
3、相同对角线上,相邻元素的索引号相差
在计算次序上,依次计算主对角线,上三角区域的第二斜对角斜对角线,第三对角线…,下三角区域的第二对角线,第三对角线…
//*****************************************************/**功能:* 方阵的Toeplitz化处理* *参数:* 参数1---输入,方阵的行/列数(方阵大小为m*m)* 参数2---输入,待处理协方差矩阵* 参数3---输出,处理后矩阵* *返回值:* 无* *附加说明:* **///******************************************************void toeplitzComplexMat(int m,complex<double>* indata,complex<double>* outdata){ ... int upIndex = m+1;//对角线相邻元素间索引间隔 ... //上斜对角线 for(int i = 0;i<m;i++) { //计算斜对角线平均值 indexA = i;//对角线起始元素的索引号 ... contrlIndex = m-i;//当前对角线上元素总数 while(useIndex<contrlIndex) { tempData += indata[indexA]; useIndex++; indexA += upIndex;//下一元素位置 } avgData.real(tempData.real()/useIndex);//幅度平均 avgData.imag(tempData.imag()/useIndex); //斜对角线赋值 ... while(useIndex<contrlIndex) { outdata[indexA] = avgData; indexA += upIndex; ... } } //下斜对角线 for(int i = 1;i<m;i++) { //计算斜对角线平均值 ... //斜对角线赋值 ... }}
共轭转置相乘
波束形成的主要任务是计算波束加权值,使波束输出信号具有指向性从而抑制非期望方向的干扰。常规波束形成的波束输出计算式如下:
注意这里的
//*****************************************************/**功能:* 计算A的共轭转置矩阵与B矩阵的乘积* *参数:* 参数1---输入,A和B矩阵的行数(A^H若能与B相乘,则A、B矩阵的行数必须相等)* 参数2---输入,A矩阵的列数* 参数3---输入,B矩阵的列数* 参数4---输出,A^H*B处理后的结果矩阵* *返回值:* 无* *附加说明:* *///******************************************************void conjugateTranMulComplexMat(int m,int aLen,int bLen,complex<double> *dataInA,complex<double>* dataInB,complex<double>* dataOut)//计算A矩阵共轭转置乘上B矩阵{ ... //计算A^H*B int indexA = 0;//A矩阵的共轭矩阵的索引值 int indexB = 0;//B矩阵的索引值 int dataAlen = m*aLen; for(int i=0;i<dataAlen;i++) { tempData[i].imag(-dataInA[i].imag());//拷贝并进行共轭处理 tempData[i].real(dataInA[i].real()); } for(int i=0;i<aLen;i++) { for(int j=0;j<bLen;j++) { indexA = i; indexB = j; ... while(num < m) { temp += tempData[indexA]*dataInB[indexB]; indexA += aLen; indexB += bLen; num++; } dataOut[index] = temp; index++; } } delete[] tempData; tempData = NULL;}
复厄米矩阵特征值求解
这里注意,协方差矩阵为复数矩阵,需要进行复数矩阵的特征分解。这里使用
//*****************************************************/**功能:* 计算m*m复数方阵的特征值和特征向量,并进行排序* *参数:* 参数1---输入,方阵的行数/列数* 参数2---输入,协方差数据* 参数3---输出,输入多通道信号* 参数4---输出,处理后矩阵* *返回值:* 无* *附加说明:* 使用前需提前配置GSL库,并包含对应头文件*///******************************************************void eigenComplexMat(int m,complex<double>*indata,double*outEigenValue,complex<double>*outEigenVec,int mode){ ... gsl_eigen_hermv_workspace *w = gsl_eigen_hermv_alloc(m);//分配计算特征值和特征向量的工作区 ... //复制一次,防止后续计算修改数据 tempA->data = (double*)indata; gsl_matrix_complex_memcpy(dataM,tempA); ... gsl_eigen_hermv(dataM,eval,evec,w);//复厄米特矩阵求解 switch(mode) { case 0: gsl_eigen_hermv_sort (eval, evec,GSL_EIGEN_SORT_ABS_ASC);//将特征列向量,按特征值升排序组成新矩阵 break; case 1: gsl_eigen_hermv_sort (eval, evec,GSL_EIGEN_SORT_ABS_DESC);//将特征列向量,按特征值降排序组成新矩阵 break; ... } //释放内存 ...}
复矩阵求逆
1、对矩阵
2、分别对
3、将2中计算得到的逆矩阵相乘。
//*****************************************************/**功能:* 计算m*m复数方阵的逆矩阵* *参数:* 参数1---输入,方阵的行数/列数* 参数2---输入,协方差数据* 参数3---输出,处理后逆矩阵* *返回值:* 无* *附加说明:* 使用前需提前配置GSL库,并包含对应头文件*///******************************************************void inverseComplexMat(int m,complex<double>* indata,complex<double>* outdata)//计算m*m的复数矩阵的逆矩阵{ ... tempA->data = (double*)indata; gsl_matrix_complex_memcpy(dataM,tempA);//复制一次,防止后续计算修改数据 int sign = 0; gsl_linalg_complex_LU_decomp(dataM, p, &sign); inverse->data = (double*)outdata; gsl_linalg_complex_LU_invert(dataM, p, inverse); //释放内存 ...}
详细代码
链接:https://pan.baidu.com/s/1c1R7DRi 密码:2wb8
参考文献
[1]http://blog.sciencenet.cn/blog-999739-779994.html
[2]https://wenku.baidu.com/view/8d513ff0941ea76e58fa04b5.html
[3]http://blog.csdn.net/piaoxuezhong/article/details/78357905
[4]https://wenku.baidu.com/view/a550fe0403d8ce2f006623da.html
[5]http://blog.csdn.net/xx_123_1_rj/article/details/39553809
[6]https://wenku.baidu.com/view/d8cf99300b1c59eef9c7b433.html
欢迎转载
转载请保留声明,并注明出处:http://blog.csdn.net/qq_31436943/article/details/78813594
- 阵列信号处理中若干算法的C++实现之基础函数篇
- C中信号处理函数
- 阵列信号处理
- 阵列信号处理
- 基于chirp信号的阵列信号处理研究随记
- Linux C 函数参考之信号处理篇
- 基于阵列信号处理的矩阵基础知识心得(持续更新中)
- 2信号处理之:信号产生原因,进程处理信号行为,信号集处理函数,PCB的信号集,sigprocmask()和sigpending(),信号捕捉设定,sigaction,C标准库信号处理函数,可重入函数,
- 阵列信号处理工具箱DBT
- 基于麦克风阵列的语音信号处理技术
- linux常用c函数 信号处理篇
- 信号和槽函数的基础实现
- 用Ruby和C实现的算法若干
- Opencv图像处理实现的虚拟16车摆阵列算法
- 【C语言】【unix c】改变信号的处理函数
- c函数之【信号函数】
- 协方差矩阵 阵列处理基础
- 浅谈语音信号处理系列之二 语音信号处理的基础
- PHP读取文件内容的五种方式
- Tomca详解
- cool -- 安全软件公司Avast开源化机器码反编译器RetDec
- hello
- iOS仿keep5.3.0版本健身等级界面实现UITableView头部图片下拉纵向拉伸效果
- 阵列信号处理中若干算法的C++实现之基础函数篇
- 使用线程和线程池
- 【Python】线程
- SpringBoot用JdbcTemplates访问Mysql
- 使用POI进行Excel导入时解决的一些问题
- POJ1459
- ’17读经感言集
- 关于Apache ad 下载总结
- 培训第八天,接口与多态