第四章MATLAB數據分析與多項式計算4.1數據統(tǒng)計處理4.2數據插值4.3曲線擬合4.4離散傅立葉變換4.5多項式計算4.1數據統(tǒng)計處理4.1.1最大值和最小值MATLAB提供的求數據序列的最大值和最小值的函數分別為max和min，兩個函數的調用格式和操作過程類似。1．求向量的最大值和最小值求一個向量X的最大值的函數有兩種調用格式，分別是：(1)y=max(X)：返回向量X的最大值存入y，如果X中包含復數元素，則按模取最大值。(2)[y,I]=max(X)：返回向量X的最大值存入y，最大值的序號存入I，如果X中包含復數元素，則按模取最大值。求向量X的最小值的函數是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的最大值的函數有3種調用格式，分別是：(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時，該函數和max(A)完全相同；dim取2時，該函數返回一個列向量，其第i個元素是A矩陣的第i行上的最大值。求最小值的函數是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．兩個向量或矩陣對應元素的比較函數max和min還能對兩個同型的向量或矩陣進行比較，調用格式為：(1)U=max(A,B)：A,B是兩個同型的向量或矩陣，結果U是與A,B同型的向量或矩陣，U的每個元素等于A,B對應元素的較大者。(2)U=max(A,n)：n是一個標量，結果U是與A同型的向量或矩陣，U的每個元素等于A對應元素和n中的較大者。min函數的用法和max完全相同。例4-3求兩個2×3矩陣x,y所有同一位置上的較大元素構成的新矩陣p。4.1.2求和與求積數據序列求和與求積的函數是sum和prod，其使用方法類似。設X是一個向量，A是一個矩陣，函數的調用格式為：sum(X)：返回向量X各元素的和。prod(X)：返回向量X各元素的乘積。sum(A)：返回一個行向量，其第i個元素是A的第i列的元素和。prod(A)：返回一個行向量，其第i個元素是A的第i列的元素乘積。sum(A,dim)：當dim為1時，該函數等同于sum(A)；當dim為2時，返回一個列向量，其第i個元素是A的第i行的各元素之和。prod(A,dim)：當dim為1時，該函數等同于prod(A)；當dim為2時，返回一個列向量，其第i個元素是A的第i行的各元素乘積。例4-4求矩陣A的每行元素的乘積和全部元素的乘積。4.1.3平均值和中值求數據序列平均值的函數是mean，求數據序列中值的函數是median。兩個函數的調用格式為：mean(X)：返回向量X的算術平均值。median(X)：返回向量X的中值。mean(A)：返回一個行向量，其第i個元素是A的第i列的算術平均值。median(A)：返回一個行向量，其第i個元素是A的第i列的中值。mean(A,dim)：當dim為1時，該函數等同于mean(A)；當dim為2時，返回一個列向量，其第i個元素是A的第i行的算術平均值。median(A,dim)：當dim為1時，該函數等同于median(A)；當dim為2時，返回一個列向量，其第i個元素是A的第i行的中值。例4-5分別求向量x與y的平均值和中值。4.1.4累加和與累乘積在MATLAB中，使用cumsum和cumprod函數能方便地求得向量和矩陣元素的累加和與累乘積向量，函數的調用格式為：cumsum(X)：返回向量X累加和向量。cumprod(X)：返回向量X累乘積向量。cumsum(A)：返回一個矩陣，其第i列是A的第i列的累加和向量。cumprod(A)：返回一個矩陣，其第i列是A的第i列的累乘積向量。cumsum(A,dim)：當dim為1時，該函數等同于cumsum(A)；當dim為2時，返回一個矩陣，其第i行是A的第i行的累加和向量。cumprod(A,dim)：當dim為1時，該函數等同于cumprod(A)；當dim為2時，返回一個向量，其第i行是A的第i行的累乘積向量。例4-6求向量X的累加和向量。x=[1:8];

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

、‘cubic’（利用pchip函數進行三次Hermite插值）、‘spline’（三次樣條插值）。Method=“nearest”,最近鄰插值法，將與待插值點最近的已知點的函數值作為待插值點函數的估計。Method=“l(fā)inear”,線性插值法，根據相鄰已知數據點間的線性函數估計落在該區(qū)域內待插值數據點的函數值。Method=“spline”,三次樣條插值，這種方法在相鄰已知數據點間建立三次多項式函數，根據該多項式函數確定待插值數據點的函數值。Method=“cubic”,利用pchip函數分段地進行三次Hermite插值。注意：X1的取值范圍不能超出X的給定范圍，否則，會給出“NaN”錯誤。例4-10用不同的插值方法計算在π/2點的值。MATLAB中有一個專門的3次樣條插值函數Y1=spline(X,Y,X1)，其功能及使用方法與函數Y1=interp1(X,Y,X1,‘spline’)完全相同。例4-11某觀測站測得某日6:00時至18:00時之間每隔2小時的室內外溫度(℃)，用3次樣條插值分別求得該日室內外6:30至17:30時之間每隔2小時各點的近似溫度(℃)。設時間變量h為一行向量，溫度變量t為一個兩列矩陣，其中第一列存放室內溫度，第二列儲存室外溫度。命令如下：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二維數據插值在MATLAB中，提供了解決二維插值問題的函數interp2，其調用格式為：Z1=interp2(X,Y,Z,X1,Y1,'method')其中X,Y是兩個向量，分別描述兩個參數的采樣點，Z是與參數采樣點對應的函數值，X1,Y1是兩個向量或標量，描述欲插值的點。Z1是根據相應的插值方法得到的插值結果。method的取值與一維插值函數相同。X,Y,Z也可以是矩陣形式。同樣，X1,Y1的取值范圍不能超出X,Y的給定范圍，否則，會給出“NaN”錯誤。例4-12設z=x^2+y^2，對z函數在[0,1]×[0,2]區(qū)域內進行插值。例4-13某實驗對一根長10米的鋼軌進行熱源的溫度傳播測試。用x表示測量點0:2.5:10(米)，用h表示測量時間0:30:60(秒)，用T表示測試所得各點的溫度(℃)。試用線性插值求出在一分鐘內每隔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函數來求得最小二乘擬合多項式的系數，再用polyval函數按所得的多項式計算所給出的點上的函數近似值。polyfit函數的調用格式為：[P,S]=polyfit(X,Y,m)函數根據采樣點X和采樣點函數值Y，產生一個m次多項式P及其在采樣點的誤差向量S。其中X,Y是兩個等長的向量，P是一個長度為m+1的向量，P的元素為多項式系數。polyval函數的功能是按多項式的系數計算x點多項式的值，將在4.5.3節(jié)中詳細介紹。例4-14已知數據表[t,y]，試求2次擬合多項式p(t)，然后求ti=1,1.5,2,2.5,…,9.5,10各點的函數近似值。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離散傅立葉變換傅立葉變換將函數表示為不同頻率的正弦、余弦函數之和，對于離散數據，傅立葉變換為DFT，即離散傅立葉變換。4.4.1離散傅立葉變換的應用（1）計算線性卷積X(n)、H(n)，Y(n)=X(n)＊H(n)X(K)、H(n)，Y(K)=X(K)×H(K)（2）對信號進行譜分析為了利用DFT，要先對模擬信號進行時域采樣，得到X(n)=Xa(nT),再對X(n)進行DFT得到的X(K)是X(n)的傅立葉變換在頻域區(qū)間[0,2pi]上的N點等間隔采樣，這樣時域和頻域均為離散信號，便于計算機處理。4.4.2離散傅立葉變換的實現一維離散傅立葉變換函數，其調用格式與功能為：(1)fft(X)：返回向量X的離散傅立葉變換。設X的長度(即元素個數)為N，若N為2的冪次，則為以2為基數的快速傅立葉變換，否則為運算速度很慢的非2冪次的算法。對于矩陣X，fft(X)應用于矩陣的每一列。(2)fft(X,N)：計算N點離散傅立葉變換。它限定向量的長度為N，若X的長度小于N，則不足部分補上零；若大于N，則刪去超出N的那些元素。對于矩陣X，它同樣應用于矩陣的每一列，只是限定了向量的長度為N。(3)fft(X,[],dim)或fft(X,N,dim)：這是對于矩陣而言的函數調用格式，前者的功能與fft(X)基本相同，而后者則與FFT(X,N)基本相同。只是當參數dim=1時，該函數作用于X的每一列；當dim=2時，則作用于X的每一行。值得一提的是，當已知給出的樣本數N0不是2的冪次時，可以取一個N使它大于N0且是2的冪次，然后利用函數格式fft(X,N)或fft(X,N,dim)便可進行快速傅立葉變換。這樣，計算速度將大大加快。相應地，一維離散傅立葉逆變換函數是ifft。ifft(F)返回F的一維離散傅立葉逆變換；ifft(F,N)為N點逆變換；ifft(F,[],dim)或ifft(F,N,dim)則由N或dim確定逆變換的點數或操作方向。例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);%計算函數值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)|')%標記例4-16給定數學函數X(t)=12sin(2π×10t+π/4)+5cos(2π×40t)取N=128，試對t從0~1秒采樣，用fft作快速傅立葉變換，繪制相應的振幅-頻率圖。在0~1秒時間范圍內采樣128點，從而可以確定采樣周期和采樣頻率。由于離散傅立葉變換時的下標應是從0到N-1，故在實際應用時下標應該前移1。又考慮到對離散傅立葉變換來說，其振幅|F(k)|是關于N/2對稱的，故只須使k從0到N/2即可。程序如下：N=128;%采樣點數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數字濾波（1）常用函數為y=filter(B,A,x,[],dim)其中B為濾波器傳輸函數的分子多項式行向量，A為分母多項式行向量，x為輸入信號。Dim＝1表示對列向量濾波，dim＝2表示對行濾波。（2）窗函數法設計濾波器Y=fir1(N,Wn,’ftype’,window)其中N為濾波器的階數，Wn為截止頻率，其值在0－1之間，1對應奈奎斯特頻率。ftype為高通或帶阻等類型，window為窗函數類型。(提一下freqz)例4-17設計一個低通濾波器，獲取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)用窗函數法實現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)分析常用函數S=ss(A,B,C,D)微分方程生成LTI模型S=zpk(z,p,k)零極點及增益生成LTI模型S=tf(num,den)傳輸函數生成LTI模型[y,t]=step(s_lti)計算系統(tǒng)的階躍響應[y,t]=impulse(s_lti)計算系統(tǒng)的沖激響應該函數實現模塊化處理，用simulink實現更方便4.7多項式計算多項式的表示：在MATLAB中，用向量表示多項式，向量的元素為多項式的系數按階數的降序排列。例如，對多項式p(x)=anxn+an-1xn-1+…+a1x+a0多項式表示為p=[an，an-1，…a1，a0]4.7.1多項式的四則運算1．多項式的加減運算2．多項式乘法運算函數conv(P1,P2)用于求多項式P1和P2的乘積。這里，P1、P2是兩個多項式系數向量。例4-16求多項式x4+8x3-10與多項式2x2-x+3的乘積。p1=[1800-10];p2=[2-13];conv(p1,p2)ans=215-524-2010-303．多項式除法函數[Q,r]=deconv(P1,P2)用于對多項式P1和P2作除法運算。其中Q返回多項式P1除以P2的商式，r返回P1除以P2的余式。這里，Q和r仍是多項式系數向量。deconv是conv的逆函數，即有P1=conv(P2,Q)+r。例4-17求多項式x4+8x3-10除以多項式2x2-x+3的結果。p1=[1800-10];p2=[2-13];[q,r]=deconv(p1,p2)q=0.50004.25001.3750r=000-11.3750-14.12504.7.2多項式的導函數對多項式求導數的函數是：p=polyder(P)：求多項式P的導函數p=polyder(P,Q)：求P·Q的導函數[p,q]=polyder(P,Q)：求P/Q的導函數，導函數的分子存入p，分母存入q。上述函數中，參數P,Q是多項式的向量表示，結果p,q也是多項式的向量表示。例4-18求有理分式的導數。命令如下：P=[1];Q=[1,0,5];[p,q]=polyder(P,Q)p=-20q=10100254.7.3多項式的求值MATLAB提供了兩種求多項式值的函數：polyval與polyvalm，它們的輸入參數均為多項式系數向量P和自變量x。兩者的區(qū)別在于前者是代數多項式求值，而后者是矩陣多項式求值。1．代數多項式求值polyval函數用來求代數多項式的值，其調用格式為：Y=polyval(P,x)若x為一數值，則求多項式在該點的值；若x為向量或矩陣，則對向量或矩陣中的每個元素求其多項式的值。例