4.11.1 一维插值
【例4.11.1-1】已知一组原始数据，确定它们所代表函数穿越y?0.95线的时刻。 （1）
t=linspace(0,5,100);y=1-cos(3*t).*exp(-t);
（2）
plot(t,y,'b');grid;hold on,plot(t,0.95*ones(size(t)),'r');hold off
1.41.210.80.60.40.20012345t_ginput=ginput(1) t_ginput = 0.4965 0.9500
图 4.11-1 （3）
it=min(find(y>0.95)); T=(it-3):(it+3); t_nearst=interp1(y(T),t(T),0.95,'nearst'); t_linear=interp1(y(T),t(T),0.95); t_cubic=interp1(y(T),t(T),0.95,'cubic'); t_spline=interp1(y(T),t(T),0.95,'spline'); disp([' t_nearst t_linear t_cubic t_spline']) disp([t_nearst t_linear t_cubic t_spline]) t_nearst t_linear t_cubic t_spline 0.5051 0.4965 0.4962 0.4962
% % % % %
<2> <3> <4> <5> <6>
（4）
t_zero=fzero('1-cos(3*x).*exp(-x)-0.95',0.5) t_zero =
0.4962
4.11.2 高维函数的插值
【例4.11.2-1】假设有一组海底深度测量数据，采用插值方式绘制海底形状图。 （1）
randn('state',2)
x=-5:5;y=-5:5;[X,Y]=meshgrid(x,y); %
zz=1.2*exp(-((X-1).^2+(Y-2).^2))-0.7*exp(-((X+2).^2+(Y+1).^2)); Z=-500+zz+randn(size(X))*0.05; surf(X,Y,Z);view(-25,25)
26
-498.5-499-499.5-500-500.5-50150-50-55图 4.11-2 （2）
xi=linspace(-5,5,50);yi=linspace(-5,5,50);[XI,YI]=meshgrid(xi,yi); ZI=interp2(X,Y,Z,XI,YI,'*cubic'); surf(XI,YI,ZI),view(-25,25)
图 4.11-3
4.12 样条函数及其应用
4.12.1 样条插值
【例4.12.1-1】根据连续时间函数w(t)?e检查重构误差。
?t的采样数据，利用spline重构该连续函数，并
t=-5:0.5:5;w=exp(-abs(t));
N0=length(t);tt=linspace(t(1),t(end),10*N0); ww=spline(t,w,tt);
error=max(abs(ww-exp(-abs(tt)))) plot(tt,ww,'b');hold on
stem(t,w,'filled','r');hold off error = 0.0840 10.90.80.70.60.50.40.30.20.10-505 27
图 4.12-1
【例4.12.1-2】用样条插值产生长、短轴分别在45度、135度线上的椭圆。
theta=[0:0.5:2]*pi;
y=[-0.5 1 -0.5 -1 0.5 1 -0.5;0.5 1 0.5 -1 -0.5 1 0.5]; %<3> theta2=linspace(theta(1),theta(end),50*length(theta)); yy=spline(theta,y,theta2);
plot(yy(1,:),yy(2,:),'b');hold on
plot(y(1,:),y(2,:),'or');hold off,axis('image') 10.80.60.40.20-0.2-0.4-0.6-0.8-1-1-0.500.51图 4.12-2
4.12.2 样条函数用于数值积分和微分
【例4.12.2-1】对于函数（1）
x=(0:0.1:1)*2*pi;y=sin(x); pp=spline(x,y); int_pp=fnint(pp); der_pp=fnder(pp); %
xx=(0:0.01:1)*2*pi;
err_yy=max(abs(ppval(pp,xx)-sin(xx)))
err_int=max(abs(ppval(int_pp,xx)-(1-cos(xx)))) err_der=max(abs(ppval(der_pp,xx)-cos(xx))) err_yy = 0.0026 err_int = 0.0010 err_der = 0.0253
y??cosx。本例将借此演示样条函数求数值不定积分、导函数的能力。
y?sinx，很容易求得S(x)??0sinxdx?1?cosx，
x（2）
%
DefiniteIntegral.bySpline=ppval(int_pp,[1,2])*[-1;1]; DefiniteIntegral.byTheory=(1-cos(2))-(1-cos(1)); %
Derivative.bySpline=fnval(der_pp,3); Derivative.byTheory=cos(3);
Derivative.byDiference=(sin(3.01)-sin(3))/0.01; DefiniteIntegral,Derivative DefiniteIntegral = bySpline: 0.9563 byTheory: 0.9564
% <2>
28
Derivative =
bySpline: -0.9895 byTheory: -0.9900 byDiference: -0.9907
（3）
fnplt(pp,'b-');hold on
fnplt(int_pp,'m:'),fnplt(der_pp,'r--');hold off legend('y(x)','S(x)','dy/dx') 2.521.510.50-0.5-1y(x) S(x) dy/dx01234567图 4.12-3
4.13 Fourier分析
4.13.1 快速Fourier变换和逆变换指令 4.13.2 连续时间函数的Fourier级数展开 4.13.2.1 展开系数的积分求取法
4.13.2.2 Fourier级数与DFT之间的数学联系 4.13.2.3 MATLAB算法实现
【例4.13.2.3-1】已知时间函数w(t)???t?0.5?00.5?t?1.5，运用符号法求该函数的
elseFourier级数展开系数。 （1）
[fzzysym.m]
function [A_sym,B_sym]=fzzysym(T,Nf,Nn) %
syms ttt n
if nargin<2;Nf=6;end if nargin<3;Nn=32;end yy=time_fun_s(ttt); A0=int(yy,ttt,0,T)/T;
As=int(yy*cos(2*pi*n*ttt/T),ttt,0,T); Bs=int(yy*sin(2*pi*n*ttt/T),ttt,0,T);
29
A_sym(1)=double(vpa(A0,Nn)); for k=1:Nf
A_sym(k+1)=double(vpa(subs(As,n,k),Nn)); B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn)); end
%------------------------------------------- function yy=time_fun_s(ttt) %
y1=sym('Heaviside(ttt-0.5)')*(ttt-0.5);
yy=y1-sym('Heaviside(ttt-1.5)')*((ttt-1.5)+1);
（2）
[A_sym,B_sym]=fzzysym(2) A_sym =
0.2500 -0.3183 0.0000 0.1061 -0.0000 -0.0637 0.0000 B_sym =
0 -0.2026 0.1592 0.0225 -0.0796 -0.0081 0.0531
【例4.13.2.3-2】运用数值积分法，按式(4.13.2.1-4)和(4.13.2.1-5)，求上例时间函数的Fourier级数展开系数。 （1）
[fzzyquad.m]
function [t,y,S,A,B]=fzzyquad(a,T,Nf,K) % %
if nargin<4;K=200;end
if (nargin<3|isempty(Nf));Nf=15;end k=1:K;
t=a+k*T/length(k); y=time_fun(t,T); S=zeros(Nf+1,length(k)); a0=mean(y); S(1,:)=a0;
A=zeros(1,Nf+1);B=zeros(1,Nf+1); A(1)=a0;
for n = 1:Nf
A(n+1) = quadl(@cos_y,a,a+T,[],[],n,T)/T*2; B(n+1) = quadl(@sin_y,a,a+T,[],[],n,T)/T*2;
S(n+1,:)=S(n,:)+A(n+1)*cos(2*n*pi*t/T)+B(n+1)*sin(2*n*pi*t/T); end
[cos_y.m]
function wcos=cos_y(t,n,T) %
y=time_fun(t,T);wcos=cos(2*n*pi*t/T).*y;
[sin_y.m]
function wsin=sin_y(t,n,T) %
y=time_fun(t,T);wsin=sin(2*n*pi*t/T).*y;
30