第4章 MATLAB數(shù)據(jù)分析與多項式計算

第四章MATLAB數(shù)據(jù)分析與多項式計算4.1數(shù)據(jù)統(tǒng)計處理4.2數(shù)據(jù)插值4.3曲線擬合4.4離散傅立葉變換4.5多項式計算4.1數(shù)據(jù)統(tǒng)計處理4.1.1最大值和最小值MATLAB提供的求數(shù)據(jù)序列的最大值和最小值的函數(shù)分別為max和min，兩個函數(shù)的調(diào)用格式和操作過程類似。1．求向量的最大值和最小值求一個向量X的最大值的函數(shù)有兩種調(diào)用格式，分別是：(1)y=max(X)：返回向量X的最大值存入y，如果X中包含復(fù)數(shù)元素，則按模取最大值。(2)[y,I]=max(X)：返回向量X的最大值存入y，最大值的序號存入I，如果X中包含復(fù)數(shù)元素，則按模取最大值。求向量X的最小值的函數(shù)是min(X)，用法和max(X)完全相同。例4-1求向量x的最大值。命令如下：x=[-43,72,9,16,23,47];y=max(x)%求向量x中的最大值[y,l]=max(x)%求向量x中的最大值及其該元素的位置2．求矩陣的最大值和最小值求矩陣A的最大值的函數(shù)有3種調(diào)用格式，分別是：(1)max(A)：返回一個行向量，向量的第i個元素是矩陣A的第i列上的最大值。(2)[Y,U]=max(A)：返回行向量Y和U，Y向量記錄A的每列的最大值，U向量記錄每列最大值的行號。(3)max(A,[],dim)：dim取1或2。dim取1時，該函數(shù)和max(A)完全相同；dim取2時，該函數(shù)返回一個列向量，其第i個元素是A矩陣的第i行上的最大值。求最小值的函數(shù)是min，其用法和max完全相同。例4-2分別求3×4矩陣x中各列和各行元素中的最大值，并求整個矩陣的最大值.x=[16930;0124334;1495736];A=max(x,[],1);%各列元素的最大值B=max(x,[],2);%各行元素的最大值Y1=max(A);Y2=max(B);Z=max(y1,y2)%矩陣最大值3．兩個向量或矩陣對應(yīng)元素的比較函數(shù)max和min還能對兩個同型的向量或矩陣進(jìn)行比較，調(diào)用格式為：(1)U=max(A,B)：A,B是兩個同型的向量或矩陣，結(jié)果U是與A,B同型的向量或矩陣，U的每個元素等于A,B對應(yīng)元素的較大者。(2)U=max(A,n)：n是一個標(biāo)量，結(jié)果U是與A同型的向量或矩陣，U的每個元素等于A對應(yīng)元素和n中的較大者。min函數(shù)的用法和max完全相同。例4-3求兩個2×3矩陣x,y所有同一位置上的較大元素構(gòu)成的新矩陣p。4.1.2求和與求積數(shù)據(jù)序列求和與求積的函數(shù)是sum和prod，其使用方法類似。設(shè)X是一個向量，A是一個矩陣，函數(shù)的調(diào)用格式為：sum(X)：返回向量X各元素的和。prod(X)：返回向量X各元素的乘積。sum(A)：返回一個行向量，其第i個元素是A的第i列的元素和。prod(A)：返回一個行向量，其第i個元素是A的第i列的元素乘積。sum(A,dim)：當(dāng)dim為1時，該函數(shù)等同于sum(A)；當(dāng)dim為2時，返回一個列向量，其第i個元素是A的第i行的各元素之和。prod(A,dim)：當(dāng)dim為1時，該函數(shù)等同于prod(A)；當(dāng)dim為2時，返回一個列向量，其第i個元素是A的第i行的各元素乘積。例4-4求矩陣A的每行元素的乘積和全部元素的乘積。4.1.3平均值和中值求數(shù)據(jù)序列平均值的函數(shù)是mean，求數(shù)據(jù)序列中值的函數(shù)是median。兩個函數(shù)的調(diào)用格式為：mean(X)：返回向量X的算術(shù)平均值。median(X)：返回向量X的中值。mean(A)：返回一個行向量，其第i個元素是A的第i列的算術(shù)平均值。median(A)：返回一個行向量，其第i個元素是A的第i列的中值。mean(A,dim)：當(dāng)dim為1時，該函數(shù)等同于mean(A)；當(dāng)dim為2時，返回一個列向量，其第i個元素是A的第i行的算術(shù)平均值。median(A,dim)：當(dāng)dim為1時，該函數(shù)等同于median(A)；當(dāng)dim為2時，返回一個列向量，其第i個元素是A的第i行的中值。例4-5分別求向量x與y的平均值和中值。4.1.4累加和與累乘積在MATLAB中，使用cumsum和cumprod函數(shù)能方便地求得向量和矩陣元素的累加和與累乘積向量，函數(shù)的調(diào)用格式為：cumsum(X)：返回向量X累加和向量。cumprod(X)：返回向量X累乘積向量。cumsum(A)：返回一個矩陣，其第i列是A的第i列的累加和向量。cumprod(A)：返回一個矩陣，其第i列是A的第i列的累乘積向量。cumsum(A,dim)：當(dāng)dim為1時，該函數(shù)等同于cumsum(A)；當(dāng)dim為2時，返回一個矩陣，其第i行是A的第i行的累加和向量。cumprod(A,dim)：當(dāng)dim為1時，該函數(shù)等同于cumprod(A)；當(dāng)dim為2時，返回一個向量，其第i行是A的第i行的累乘積向量。例4-6求向量X的累加和向量。x=[1:8];

cumsum(x)結(jié)果顯示為ans=13610152128364.1.5標(biāo)準(zhǔn)方差與相關(guān)系數(shù)1．求標(biāo)準(zhǔn)方差在MATLAB中，提供了計算數(shù)據(jù)序列的標(biāo)準(zhǔn)方差的函數(shù)std。對于向量X，std(X)返回一個標(biāo)準(zhǔn)方差。對于矩陣A，std(A)返回一個行向量，它的各個元素便是矩陣A各列或各行的標(biāo)準(zhǔn)方差。std函數(shù)的一般調(diào)用格式為：Y=std(A,flag,dim)其中dim取1或2。當(dāng)dim=1時，求各列元素的標(biāo)準(zhǔn)方差；當(dāng)dim=2時，則求各行元素的標(biāo)準(zhǔn)方差。flag取0或1，當(dāng)flag=0時，按σ1所列公式計算標(biāo)準(zhǔn)方差，當(dāng)flag=1時，按σ2所列公式計算標(biāo)準(zhǔn)方差。缺省flag=0，dim=1。例4.7求標(biāo)準(zhǔn)方差。已知A=[248;136;589],v=Std(v)%求向量V的1型標(biāo)準(zhǔn)差Std(v,1)%求向量V的1型標(biāo)準(zhǔn)差Std(A,0,2)%矩陣沿行方向的1型標(biāo)準(zhǔn)差。ans=3.05512.51662.08172．相關(guān)系數(shù)MATLAB提供了corrcoef函數(shù)，可以求出數(shù)據(jù)的相關(guān)系數(shù)矩陣。corrcoef函數(shù)的調(diào)用格式為：corrcoef(X)：返回從矩陣X形成的一個相關(guān)系數(shù)矩陣。此相關(guān)系數(shù)矩陣的大小與矩陣X一樣。它把矩陣X的每列作為一個變量，然后求它們的相關(guān)系數(shù)。corrcoef(X,Y)：在這里，X,Y是向量，它們與corrcoef([X,Y])的作用一樣。例4-8生成滿足正態(tài)分布的10000×5隨機(jī)矩陣，然后求各列元素的均值和標(biāo)準(zhǔn)方差，再求這5列隨機(jī)數(shù)據(jù)的相關(guān)系數(shù)矩陣。命令如下：X=randn(10000,5);M=mean(X)D=std(X)R=corrcoef(X)4.1.6排序MATLAB中對向量X的排序函數(shù)是sort(X)，函數(shù)返回一個對X中的元素按升序排列的新向量。sort函數(shù)也可以對矩陣A的各列或各行重新排序，其調(diào)用格式為：[Y,I]=sort(A,dim)其中dim指明對A的列還是行進(jìn)行排序。若dim=1，則按列排；若dim=2，則按行排。Y是排序后的矩陣，而I記錄Y中的元素在A中位置。例4.9A=[248;136;589];V=;sort(v)%向量V作升序排列ans=12345667910sort(V,‘descend’)%向量V作降序排列ans=10976654321sort(A)%將矩陣A按列作升序排列ans=136248589sort(A,2)%將矩陣A按行作升序排列ans=248136589[B,I]=sort(A)%將矩陣按列排序并求出每個元素的索引B=136248589I=2221113334.2數(shù)據(jù)插值4.2.1一維數(shù)據(jù)插值在MATLAB中，實現(xiàn)這些插值的函數(shù)是interp1，其調(diào)用格式為：Y1=interp1(X,Y,X1,'method')函數(shù)根據(jù)X,Y的值，計算函數(shù)在X1處的值。X,Y是兩個等長的已知向量，分別描述采樣點和樣本值，X1是一個向量或標(biāo)量，描述欲插值的點，Y1是一個與X1等長的插值結(jié)果。method是插值方法，允許的取值有‘linear’（線性插值法）、‘nearest’（最近鄰插值法）

、‘cubic’（利用pchip函數(shù)進(jìn)行三次Hermite插值）、‘spline’（三次樣條插值）。Method=“nearest”,最近鄰插值法，將與待插值點最近的已知點的函數(shù)值作為待插值點函數(shù)的估計。Method=“l(fā)inear”,線性插值法，根據(jù)相鄰已知數(shù)據(jù)點間的線性函數(shù)估計落在該區(qū)域內(nèi)待插值數(shù)據(jù)點的函數(shù)值。Method=“spline”,三次樣條插值，這種方法在相鄰已知數(shù)據(jù)點間建立三次多項式函數(shù)，根據(jù)該多項式函數(shù)確定待插值數(shù)據(jù)點的函數(shù)值。Method=“cubic”,利用pchip函數(shù)分段地進(jìn)行三次Hermite插值。注意：X1的取值范圍不能超出X的給定范圍，否則，會給出“NaN”錯誤。例4-10用不同的插值方法計算在π/2點的值。MATLAB中有一個專門的3次樣條插值函數(shù)Y1=spline(X,Y,X1)，其功能及使用方法與函數(shù)Y1=interp1(X,Y,X1,‘spline’)完全相同。例4-11某觀測站測得某日6:00時至18:00時之間每隔2小時的室內(nèi)外溫度(℃)，用3次樣條插值分別求得該日室內(nèi)外6:30至17:30時之間每隔2小時各點的近似溫度(℃)。設(shè)時間變量h為一行向量，溫度變量t為一個兩列矩陣，其中第一列存放室內(nèi)溫度，第二列儲存室外溫度。命令如下：h=6:2:18;t=[18,20,22,25,30,28,24;15,19,24,28,34,32,30]';XI=6.5:2:17.5YI=interp1(h,t,XI,‘spline’)%用3次樣條插值計算4.2.2二維數(shù)據(jù)插值在MATLAB中，提供了解決二維插值問題的函數(shù)interp2，其調(diào)用格式為：Z1=interp2(X,Y,Z,X1,Y1,'method')其中X,Y是兩個向量，分別描述兩個參數(shù)的采樣點，Z是與參數(shù)采樣點對應(yīng)的函數(shù)值，X1,Y1是兩個向量或標(biāo)量，描述欲插值的點。Z1是根據(jù)相應(yīng)的插值方法得到的插值結(jié)果。method的取值與一維插值函數(shù)相同。X,Y,Z也可以是矩陣形式。同樣，X1,Y1的取值范圍不能超出X,Y的給定范圍，否則，會給出“NaN”錯誤。例4-12設(shè)z=x^2+y^2，對z函數(shù)在[0,1]×[0,2]區(qū)域內(nèi)進(jìn)行插值。例4-13某實驗對一根長10米的鋼軌進(jìn)行熱源的溫度傳播測試。用x表示測量點0:2.5:10(米)，用h表示測量時間0:30:60(秒)，用T表示測試所得各點的溫度(℃)。試用線性插值求出在一分鐘內(nèi)每隔20秒、鋼軌每隔1米處的溫度TI。命令如下：x=0:2.5:10;h=[0:30:60]';T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];xi=[0:10];hi=[0:20:60]';TI=interp2(x,h,T,xi,hi)4.3曲線擬合在MATLAB中，用polyfit函數(shù)來求得最小二乘擬合多項式的系數(shù)，再用polyval函數(shù)按所得的多項式計算所給出的點上的函數(shù)近似值。polyfit函數(shù)的調(diào)用格式為：[P,S]=polyfit(X,Y,m)函數(shù)根據(jù)采樣點X和采樣點函數(shù)值Y，產(chǎn)生一個m次多項式P及其在采樣點的誤差向量S。其中X,Y是兩個等長的向量，P是一個長度為m+1的向量，P的元素為多項式系數(shù)。polyval函數(shù)的功能是按多項式的系數(shù)計算x點多項式的值，將在4.5.3節(jié)中詳細(xì)介紹。例4-14已知數(shù)據(jù)表[t,y]，試求2次擬合多項式p(t)，然后求ti=1,1.5,2,2.5,…,9.5,10各點的函數(shù)近似值。t=[5:10];y=[21,23,24,65,96,76];ti=[1:0.5:10];[ps]=polyfit(t,y,2)polyval(p,ti)4.4離散傅立葉變換傅立葉變換將函數(shù)表示為不同頻率的正弦、余弦函數(shù)之和，對于離散數(shù)據(jù)，傅立葉變換為DFT，即離散傅立葉變換。4.4.1離散傅立葉變換的應(yīng)用（1）計算線性卷積X(n)、H(n)，Y(n)=X(n)＊H(n)X(K)、H(n)，Y(K)=X(K)×H(K)（2）對信號進(jìn)行譜分析為了利用DFT，要先對模擬信號進(jìn)行時域采樣，得到X(n)=Xa(nT),再對X(n)進(jìn)行DFT得到的X(K)是X(n)的傅立葉變換在頻域區(qū)間[0,2pi]上的N點等間隔采樣，這樣時域和頻域均為離散信號，便于計算機(jī)處理。4.4.2離散傅立葉變換的實現(xiàn)一維離散傅立葉變換函數(shù)，其調(diào)用格式與功能為：(1)fft(X)：返回向量X的離散傅立葉變換。設(shè)X的長度(即元素個數(shù))為N，若N為2的冪次，則為以2為基數(shù)的快速傅立葉變換，否則為運算速度很慢的非2冪次的算法。對于矩陣X，fft(X)應(yīng)用于矩陣的每一列。(2)fft(X,N)：計算N點離散傅立葉變換。它限定向量的長度為N，若X的長度小于N，則不足部分補(bǔ)上零；若大于N，則刪去超出N的那些元素。對于矩陣X，它同樣應(yīng)用于矩陣的每一列，只是限定了向量的長度為N。(3)fft(X,[],dim)或fft(X,N,dim)：這是對于矩陣而言的函數(shù)調(diào)用格式，前者的功能與fft(X)基本相同，而后者則與FFT(X,N)基本相同。只是當(dāng)參數(shù)dim=1時，該函數(shù)作用于X的每一列；當(dāng)dim=2時，則作用于X的每一行。值得一提的是，當(dāng)已知給出的樣本數(shù)N0不是2的冪次時，可以取一個N使它大于N0且是2的冪次，然后利用函數(shù)格式fft(X,N)或fft(X,N,dim)便可進(jìn)行快速傅立葉變換。這樣，計算速度將大大加快。相應(yīng)地，一維離散傅立葉逆變換函數(shù)是ifft。ifft(F)返回F的一維離散傅立葉逆變換；ifft(F,N)為N點逆變換；ifft(F,[],dim)或ifft(F,N,dim)則由N或dim確定逆變換的點數(shù)或操作方向。例4-15計算線性卷積clcclearX1=[122143613];X2=[2-31-524];X3=conv(X1,X2);Y1=fft(X1,16);%取16點.Y2=fft(X2,16);Y3=Y1.*Y2;%卷積Y4=real(ifft(Y3));subplot(2,1,1);stem(X3);subplot(2,1,2);stem(Y4);例4-15求傅立葉變換clcclearN=128;%取128個點t=linspace(0,3,N);%時間點f=2*exp(-3*t);%計算函數(shù)值Ts=t(2)-t(1);%采樣周期Ws=2*pi/Ts;%采樣角頻率F=fft(f);%計算fftFp=F(1:N/2+1)*Ts;%傅立葉頻譜歸一化W=Ws*(0:N/2)/N;%取1/2采樣頻率Fa=2./(3+j*W);%解析傅立葉變換plot(W,abs(Fa),W,abs(Fp),'+')%繪圖xlabel('Frequency,Rad/s'),ylabel('|F(w)|')%標(biāo)記例4-16給定數(shù)學(xué)函數(shù)X(t)=12sin(2π×10t+π/4)+5cos(2π×40t)取N=128，試對t從0~1秒采樣，用fft作快速傅立葉變換，繪制相應(yīng)的振幅-頻率圖。在0~1秒時間范圍內(nèi)采樣128點，從而可以確定采樣周期和采樣頻率。由于離散傅立葉變換時的下標(biāo)應(yīng)是從0到N-1，故在實際應(yīng)用時下標(biāo)應(yīng)該前移1。又考慮到對離散傅立葉變換來說，其振幅|F(k)|是關(guān)于N/2對稱的，故只須使k從0到N/2即可。程序如下：N=128;%采樣點數(shù)T=1;%采樣時間終點t=linspace(0,T,N);%給出N個采樣時間ti(I=1:N)x=12*sin(2*pi*10*t+pi/4)+5*cos(2*pi*40*t);%求各采樣點樣本值xdt=t(2)-t(1);%采樣周期fs=1/dt;%采樣頻率(Hz)X=fft(x);%計算x的快速傅立葉變換XF=X(1:N/2+1);%F(k)=X(k)(k=1:N/2+1)f=fs*(0:N/2)/N;%使頻率軸f從零開始plot(f,abs(F),'-*')%繪制振幅-頻率圖xlabel('Frequency');ylabel('|F(k)|')4.5數(shù)字濾波（1）常用函數(shù)為y=filter(B,A,x,[],dim)其中B為濾波器傳輸函數(shù)的分子多項式行向量，A為分母多項式行向量，x為輸入信號。Dim＝1表示對列向量濾波，dim＝2表示對行濾波。（2）窗函數(shù)法設(shè)計濾波器Y=fir1(N,Wn,’ftype’,window)其中N為濾波器的階數(shù)，Wn為截止頻率，其值在0－1之間，1對應(yīng)奈奎斯特頻率。ftype為高通或帶阻等類型，window為窗函數(shù)類型。(提一下freqz)例4-17設(shè)計一個低通濾波器，獲取10Hz信號clear,randn('state',1)ws=1024;t=0:1/ws:0.4;x=sin(2*pi*10*t)+cos(2*pi*100*t)+0.2*randn(size(t));wn=ws/2;[BA]=butter(10,30/wn);y=filter(B,A,x);yy=fft(y,512);subplot(2,1,1)f=yy(1:ws/2);ff=0:ws/2-1;plot(ff,abs(yy))gridonsubplot(2,1,2)plot(t,x,'b-',t,y,'r.','markersize',10);legend('input','output',0)用窗函數(shù)法實現(xiàn)clear,randn('state',1)ws=1024;t=0:1/ws:0.4;x=sin(2*pi*10*t)+cos(2*pi*100*t)+0.2*randn(size(t));y=fir1(64,0.06);%freqz(y);yy=conv(x,y);subplot(2,1,1);plot(t,x,'b-');subplot(2,1,2)plot(yy,'r.','markersize',10);legend('input','output',0)4-6系統(tǒng)分析常用函數(shù)S=ss(A,B,C,D)微分方程生成LTI模型S=zpk(z,p,k)零極點及增益生成LTI模型S=tf(num,den)傳輸函數(shù)生成LTI模型[y,t]=step(s_lti)計算系統(tǒng)的階躍響應(yīng)[y,t]=impulse(s_lti)計算系統(tǒng)的沖激響應(yīng)該函數(shù)實現(xiàn)模塊化處理，用simulink實現(xiàn)更方便4.7多項式計算多項式的表示：在MATLAB中，用向量表示多項式，向量的元素為多項式的系數(shù)按階數(shù)的降序排列。例如，對多項式p(x)=anxn+an-1xn-1+…+a1x+a0多項式表示為p=[an，an-1，…a1，a0]4.7.1多項式的四則運算1．多項式的加減運算2．多項式乘法運算函數(shù)conv(P1,P2)用于求多項式P1和P2的乘積。這里，P1、P2是兩個多項式系數(shù)向量。例4-16求多項式x4+8x3-10與多項式2x2-x+3的乘積。p1=[1800-10];p2=[2-13];conv(p1,p2)ans=215-524-2010-303．多項式除法函數(shù)[Q,r]=deconv(P1,P2)用于對多項式P1和P2作除法運算。其中Q返回多項式P1除以P2的商式，r返回P1除以P2的余式。這里，Q和r仍是多項式系數(shù)向量。deconv是conv的逆函數(shù)，即有P1=conv(P2,Q)+r。例4-17求多項式x4+8x3-10除以多項式2x2-x+3的結(jié)果。p1=[1800-10];p2=[2-13];[q,r]=deconv(p1,p2)q=0.50004.25001.3750r=000-11.3750-14.12504.7.2多項式的導(dǎo)函數(shù)對多項式求導(dǎo)數(shù)的函數(shù)是：p=polyder(P)：求多項式P的導(dǎo)函數(shù)p=polyder(P,Q)：求P·Q的導(dǎo)函數(shù)[p,q]=polyder(P,Q)：求P/Q的導(dǎo)函數(shù)，導(dǎo)函數(shù)的分子存入p，分母存入q。上述函數(shù)中，參數(shù)P,Q是多項式的向量表示，結(jié)果p,q也是多項式的向量表示。例4-18求有理分式的導(dǎo)數(shù)。命令如下：P=[1];Q=[1,0,5];[p,q]=polyder(P,Q)p=-20q=10100254.7.3多項式的求值MATLAB提供了兩種求多項式值的函數(shù)：polyval與polyvalm，它們的輸入?yún)?shù)均為多項式系數(shù)向量P和自變量x。兩者的區(qū)別在于前者是代數(shù)多項式求值，而后者是矩陣多項式求值。1．代數(shù)多項式求值polyval函數(shù)用來求代數(shù)多項式的值，其調(diào)用格式為：Y=polyval(P,x)若x為一數(shù)值，則求多項式在該點的值；若x為向量或矩陣，則對向量或矩陣中的每個元素求其多項式的值。例