利用Matlab工具箱進(jìn)行數(shù)字信號(hào)處理略談
相關(guān)鏈接： 中國安全網(wǎng) 中國質(zhì)量網(wǎng) 中國論文網(wǎng) 中國資訊網(wǎng)
| 摘要： | 本文詳述了在Vc環(huán)境下如何利用Matlab工具箱進(jìn)行數(shù)字信號(hào)處理，全文以Matlab工具箱中功率譜密度分析函數(shù)為例，介紹了通過Matlab自帶的引擎、Matlab自身的編譯器以及利用MathTools公司的Matcom進(jìn)行對(duì)工具箱函數(shù)的調(diào)用。 |
關(guān)鍵詞：Matlab M-文件 引擎 編譯器 Matcom Vc++
Matlab的信號(hào)處理工具箱是信號(hào)算法文件的集合，它處理的基本對(duì)象是信號(hào)與系統(tǒng)，信號(hào)處理工具箱位于目錄、ToolboxSignal下，利用工具箱中的文件可以實(shí)現(xiàn)信號(hào)的變換、濾波、譜估計(jì)、濾波器設(shè)計(jì)等。在其它的環(huán)境如Vc下如果能調(diào)用Matlab工具箱中的文件，會(huì)大大地加快一些算法的實(shí)現(xiàn)，同時(shí)其可靠性也很高。
利用Matlab引擎Matlab引擎采用客戶和服務(wù)器計(jì)算方式，在運(yùn)用中，Vc的C語言或C++語言的程序作為前端客戶機(jī)，它向Matlab引擎?zhèn)鬟f命令和數(shù)據(jù)信息，并從Matlab引擎接收數(shù)據(jù)信息，它提供了下列幾個(gè)函數(shù)： engOpen, engGetArray, engPutArray, engEvaString,
engOutputBuffer ,engClose與客戶機(jī)進(jìn)行交互。
下面例程是在Vc下建一個(gè)基于對(duì)話框的應(yīng)用程序，在對(duì)話框中設(shè)置一個(gè)Button控件OnMatlabEngine.,在對(duì)話框 .cpp文件中加入”engine.h” 和“math.h” 頭文件，下面給出部分程序清單。
Void CtestmatlabDlg::OnMatlabEngine(){
Engine *ep;
mxArray* T=NULL,*result=NULL,*mFs=NULL,*mnfft= NULL;
double datax[1024];
char buffer[1024];
for(int j=0;j<1024;j++)//注：如通過采集卡采集數(shù)據(jù)可將采集的數(shù)據(jù)放在datax[]數(shù)組中,此循環(huán)就不需要
{
double samt=(double)(1.0/1024);
datax[j]=sin(2.0*63.0*samt*3.1415926+1.15*3.1415926);
}
double *pPxx,*pFxx;
if(!(ep=engOpen(" "))){//打開Matlab引擎，建立與本地Matlab的連接
fprintf(stderr,"n Can't start MATLAB enginen");
exit(-1);
}
double Fs[1]={1024};//因?yàn)镸atlab所有參與運(yùn)算的參數(shù)都是矩陣的形式，因而下列幾行將參數(shù)轉(zhuǎn)變
double nfft[1]={1024};//成Matlab可接受的矩陣形式。
T=mxCreateDoubleMatrix(1,1024,mxREAL);
mnfft=mxCreateDoubleMatrix(1,1,mxREAL);
mFs=mxCreateDoubleMatrix(1,1,mxREAL);
mxSetName(T,"T");
mxSetName(mnfft,"mnfft");
mxSetName(mFs,"mFs");
memcpy((char*)mxGetPr(T),(char*)datax, 1024*sizeof(double));
memcpy((char*)mxGetPr(mnfft),(char*)nfft, sizeof(double));
memcpy((char*)mxGetPr(mFs),(char*)Fs,1*sizeof(double));
engPutArray(ep,T); //將轉(zhuǎn)化的參數(shù)放入引擎中，此時(shí)可在Matlab command窗口下查看此參數(shù)
engPutArray(ep,mnfft);
engPutArray(ep,mFs);
engEvalString(ep,"[pxx,fo]=psd(T,mnfft,mFs);"); //利用引擎執(zhí)行工具箱中文件
engOutputBuffer(ep,buffer,512); //如只想看顯示圖形，可將返回參數(shù)去掉，psd無返回參數(shù)缺省情況下會(huì)自動(dòng)畫圖形
result=engGetArray(ep,"pxx");//取出引擎中的數(shù)據(jù)放在所指的區(qū)域中供后續(xù)處理
pPxx=mxGetPr(result);
result=engGetArray(ep,"fo");
pFxx=mxGetPr(result);
engEvalString(ep,"plot(fo,10*log10(pxx));");//利用引擎畫圖
engEvalString(ep,"title('功率譜分析');");
engEvalString(ep,"xlabel('Hz');");
engEvalString(ep,"ylable('db');");
mxDestroyArray(T); //釋放內(nèi)存
mxDestroyArray(mFs);
mxDestroyArray(mnfft);
mxDestroyArray(result);
engEvalString(ep,"close;");
engClose(ep);
}
上述程序在Vc下編譯需要將 libeng.dll和libmx.dll兩個(gè)動(dòng)態(tài)庫利用以下的命令：
lib/def:<自己的Matlab的安裝路徑,下同>e: Matlabexterninclude*.def /machine:ix86 /out:*.lib來生成程序所需的靜態(tài)連接庫libeng.lib和libmx.lib,將libeng.lib和libmx.lib所在的目錄加入Vc++ project/link/object/library modules下即可。
利用Matlab自身的編譯器調(diào)用工具箱中的函數(shù)
Matlab的編譯器可將Matlab的M文件轉(zhuǎn)換為為C或C++的源代碼以產(chǎn)生完全脫離Matlab運(yùn)行環(huán)境的獨(dú)立的運(yùn)用程序，但Matlab本身的資料說明編譯器如用來建立獨(dú)立的運(yùn)用程序，不能調(diào)用Matlab工具箱中的函數(shù)，這非常不利于搞一些特殊的算法。本人研究了一段時(shí)間發(fā)現(xiàn)，工具箱中的函數(shù)既然是M文件就一定可以用編譯器來編譯，以提供如Vc的調(diào)用函數(shù)，但是編譯器只能編譯一個(gè)獨(dú)立的M文件，即這個(gè) M文件不依賴于其他的M文件。如果M文件中又調(diào)用了其他的M文件，可將被調(diào)用的M文件拷貝到調(diào)用M文件的相應(yīng)位置，作適當(dāng)?shù)母膭?dòng)就可以用于編譯器編譯。編譯器不支持圖形函數(shù)，所以M文件中如有圖形函數(shù)需注釋掉。
當(dāng)Matlab的編譯器mcc加入適當(dāng)?shù)膮?shù)-e(mcc –e *.*)或-p(mcc –p *.*)就可生成將輸入的M文件轉(zhuǎn)換為適用于特定運(yùn)用的C或C++源代碼。這樣如果要在Vc下編譯通過,還需連入以下幾個(gè)庫libmmfile.dll, libmatlb.dll, libmcc.dll, libmat.dll. libmx.dll. mibut.dll 以及Matlab C MATH庫，建議采用前述的方法將動(dòng)態(tài)連接改為靜態(tài)連接。對(duì)于C/C++編譯環(huán)境的設(shè)置，在Matlab command窗口下運(yùn)行mex –setup 然后依提示操作，而對(duì)于C/C++連接環(huán)境的設(shè)置，運(yùn)行mbuild –setup依提示操作即可。
下面給出利用編譯器將Matlab工具箱中psd.m文件生成可供Vc調(diào)用的函數(shù)。
將psd.m文件拷貝一份至Matlabbin目錄下，改寫相應(yīng)調(diào)用的M文件如nargchk.m, hanning.m等。為生成的代碼簡潔，對(duì)于采集數(shù)據(jù)處理輸入?yún)?shù)很明了的情況下可作大量的刪減，最終使psd.m成為一個(gè)不依賴于其他M文件的獨(dú)立的M文件，注意千萬注釋掉作圖代碼，最終改成如下形式，限于篇幅給出關(guān)鍵的幾步：function [Pxx,f]=psd(Fs,nfft,noverlap,x)
window=o,5*(1-cos(2*pi*(1:nfft)’/(nffft+1)));//hanning 窗
dflag=’none’;
window=window(;)
………………………………….
以上只要稍懂Matlab語言和信號(hào)處理知識(shí)就可完成這項(xiàng)工作。
假設(shè)上述代碼重新存為testwin.m,在Matlab command 窗口下設(shè)置好環(huán)境參數(shù)運(yùn)行mcc –e testwin,則可在Matlabbin下生成testwin.c ,如運(yùn)行mcc –p testwin 則生成testwin.cpp. Vc下建立一個(gè)基于對(duì)話框的文件，然后在對(duì)話框里加一個(gè)Button控件OnButtonPsd將上述生成的.c文件的頭文件加入到工程的.cpp中，且將#ifdef_cplusplus
extern “c”{
#end if
c代碼聲明加入Vc的包含文件和生成的.C的包含文件之間
將#ifdef_cplusplus
}
#end if加入.cpp文件未尾
為了簡潔且便于處理將生成的c函數(shù)稍改動(dòng)，給出部分代碼如下： void CTestpsdwinDlg::OnButtonPsd(){mxArray* x_rhs_;//指向采集數(shù)據(jù)存放區(qū)
Fs=23510;//數(shù)據(jù)采集的頻率 nfft=1024;//1024點(diǎn)的fft
double datax[1024]//采集的數(shù)據(jù)
x_rhs_mxCreateDoubleMatrix(1,1024,mxReal);
memcpy(mxGetPr(x_rhs_),datax,1024*sizeof(double));
noverlap=512;
……………….
……………….
mccCopy(&Pxx,&Spec);
mccCopy(&f,&frevgg_vector);
for(int j=0;j<(int)(nfft/2+1);j++)
{
datap[j]=mccGetRealVectorElement(&Pxx, (j+1));//功率譜密度存于datap[]數(shù)組
dataf[j]=mccGetRealVectorElement(&f, (j+1));//相應(yīng)頻率存于數(shù)組dataf[]中
}
mccFreeMatrix(&Pxx);
……………….
SendMessageBox(WM_PAINT,0,0);//利用Vc下的圖形函數(shù)畫圖
Return;
}
如上生成的程序可讀性不太好，而生成的c++代碼則可讀性較好，但千萬注意只能用 Matlab的MATH庫，不可用c++的MATH庫，否則編譯會(huì)出錯(cuò)，限于篇幅在此不述。
3）利用Matcom調(diào)用工具箱中的函數(shù)
Matcom編譯M文件，先將M文件按照與Matcom的cpp庫的對(duì)應(yīng)關(guān)系翻譯為cpp源代碼，然后用對(duì)應(yīng)版本的c編譯器將cpp文件編譯成相應(yīng)的exe或dll文件，所以第一次運(yùn)行要指定c編譯器的路徑，否則無法編譯，指定好的編譯信息就寫在Matcombinmatcom.ini文件中，不過這一步按裝matcom時(shí)，它自動(dòng)尋找編譯器并將其寫入matcom.ini文件中，matcom4.5版中使用TeeChart3.0 OCX控件，因而它支持圖形操作。
我們依然用上述的testwin.m文件，不要將圖形函數(shù)注釋掉，利用Mideva來生成可被Vc調(diào)用的信號(hào)處理程序。
運(yùn)行Mideva在主界面上直打開M文件，在菜單中選擇compile to dll,輸入testwin..在Matcom debug目錄下可以找到這樣的幾個(gè)文件，testwin.c ,testwin.h,testwin.cpp,testwin.lib,testwin.dll,testwin.exp等。 將上述testwin.cpp和testwin.h加入工程中，project/add to project/files并且在相應(yīng)的文件中加入”stdafx.h” 加連接庫：Toolsoptiondirectory ， 選include選項(xiàng),加入e:matcom45lib (包含matcom.h)library選項(xiàng)，加入e:matcom45lib
4) projectadd to projectfiles 文件類型選項(xiàng)選（.lib）將e:matcom45libv4501.lib加入工程中編譯運(yùn)行。相應(yīng)代碼如下：
void CtestmatcomDlg::OnpsdButton(){
double datap[512],dataf[512];
initM(MATCOM_VERSION);//初始化matcom庫
Mm Fs,nfft,noverlap;//創(chuàng)建矩陣
Mm x=zeros(1,1024);
Fs=1024;nfft=1024;noverlap=128;
dMm(Pxx_o);dMm(f_o);//創(chuàng)建并命名矩陣
datax[];//數(shù)據(jù)采集的數(shù)據(jù)存于此數(shù)組中
for(int i=1;i<=1024;i++)
{
x.r(1,i)=datax[i+1];//給x陣賦值
}
testwin(Fs,nfft,noverlap,x,i_o,Pxx,f_o);//matcom生成的函數(shù)
for(i=0;i<513;i++){//取出功率譜密度分析結(jié)果
dataf[i]=f_o.r(i+1,1);
datap[i]=Pxx_o.r(i+1,1);}
exitM();
return;
}
可見利用Matcom進(jìn)行M文件轉(zhuǎn)換非常的容易，生成的代碼可讀性很好，以上的轉(zhuǎn)換同時(shí)生成了可供Vc調(diào)用的動(dòng)態(tài)連接庫，其使用和一般的動(dòng)態(tài)庫一樣使用。同時(shí)需指明Matcom不僅可轉(zhuǎn)換獨(dú)立的不依賴于其它M文件的M文件，同時(shí)可轉(zhuǎn)換調(diào)用其它M文件的M文件嵌套。條件是這此M文件在同一個(gè)目錄下面，如前所述的psd.m可直接用上述方法轉(zhuǎn)換，生成了多個(gè)重載形式的psd函數(shù)
結(jié)論： 利用Mtlab引擎調(diào)用工具箱中的函數(shù)可節(jié)省大量的系統(tǒng)資源，應(yīng)用程序整體性能較好，但不可脫離Matlab 的環(huán)境運(yùn)行。用Matlab編譯器進(jìn)行工具箱函數(shù)的調(diào)用，須轉(zhuǎn)換相應(yīng)的Ｍ文件使其成為獨(dú)立的Ｍ文件，且不支持圖形函數(shù)，轉(zhuǎn)換的代碼可讀性不太好。用Matcom 進(jìn)行轉(zhuǎn)換非常方便，生成的代碼可讀性很好，支持圖形函數(shù)，且代碼執(zhí)行的速度比不轉(zhuǎn)換平均要快１.5倍以上。以上程序在Vc++ 6.0,Matlab5.2,Matcom4.5中調(diào)試通過,以上方法在工程實(shí)踐中已得到很好的運(yùn)用。
