通信系统仿真(一些简单的MATLAB函数学习)（多项式插值与拟合）(一)

最新推荐文章于 2024-04-08 09:57:52 发布

daijingxin

最新推荐文章于 2024-04-08 09:57:52 发布

阅读量999

点赞数 1

分类专栏：通信系统

本文链接：https://blog.csdn.net/daijingxin/article/details/103482265

版权

通信系统专栏收录该内容

28 篇文章 223 订阅

订阅专栏

一、创建向量的函数

1.线性增量序列模式

基本格式
linspace(startValue, endValue)
linspace(startValue, endValue, nElements)
注解: 当 nElements缺省时，nElements=100

例如：

u=linspace(0,0.25,5)
u=
0 0.0625 0.1250 0.1850 0.2500

2.指数增量序列模式

基本格式
logspace(startValue, endValue)
logspace(startValue, endValue, nElements)
注解: 从 10^startValue 到 10^endValue

u=logspace(1,4,4)
u=
10 100 1000 10000

3.使用冒号创建向量

基本格式
startValue: endValue
startValue: increment: endValue
注解：缺省增量为单位值1

s=1:4
s =
1 2 3 4

t=0:0.1:0.4
t=
0 0.1 0.2 0.3 0.4

%计算函数值z=...,x和y取值区间为[-8,8]，等间隔取10个点
x = linspace(-8,8,10);
y = linspace(-8,8,10);
z = (sin(sqrt(x.^2 + y.^2)))./(sqrt(x.^2 + y.^2))
%计算函数y=...的数值，t取值区间为[0,4*pi],取值间隔为0.1pi
t = (0:0.1:4)*pi;
y = exp(t./3).*sin(t);

二、多项式与插值

1.多项式的表达

数学表达式
y=c1x+c2x+c3x+c4x+…+cn*x
求值计算

c=[c1,c2,c3,…,cn,];
y=polyval(c,x)

2.多项式的插值

(1).一维插值

yi=interp1(x,y,xi,’method’)
x,y为行向量
xi为插值范围内任意点的x坐标
yi为插值运算后的对应y坐标
method— linear (默认)/ nearest/ spline/ previous / next
/ pchip/ cubic
在这里插入图片描述
例如：
x0=0:0.12:1;
y0=(x0.^2-3*x0+5).exp(-5x0).*sin(x0);
plot(x0,y0,’ro-’)
x1=0:0.02:1;
y1=interp1(x0,y0,x1,’linear’);
y2=interp1(x0,y0,x1,’cubic’);
y3=interp1(x0,y0,x1,’spline’);
y4=interp1(x0,y0,x1,’nearest’);
plot(x1,y1,‘k’);
.
.
.
plot(x0,y0)
legend(‘linear’,’cubic’,’spline’,’nearest’,’样本点’，‘原函数’ );
在这里插入图片描述

(2).二维插值

zi= griddata(x,y,z,xi,yi,method);
例：
zi= griddata(x,y,z,xi,yi,‘nearest’);
mesh(xi,yi,zi)
hold on
plot3(x,y,z,’mo’)
在这里插入图片描述
interp2能较好的进行二维运算，但只能处理以网格形式给出的数据，通过使用mesh(xi,yi,zi)函数先进行处理。

(3).多维插值

三维插值：
interp3(x,y,z,v,xi,yi,zi,method)
n维插值：
interpn(x1,x2,x3…,v,y1,y2,y3,…,method)
n维空间栅格
[X1,X2,X3,…]=ndgrid(x1,x2,x3,…)

二、多项式拟合

问题：给定一批数据点，需确定满足特定要求的曲线或曲面
解决方案：
• 若要求所求曲线（面）通过所给所有数据点，就是插值问题；
• 若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，就是数据拟合。

1.拟合的基本原理

曲线拟合问题最常用的解法——线性最小二乘法的基本思路

2.用MATLAB解决拟合问题

(1).线性最小二乘拟合

用一个多项式来逼近一组给定数据：
p=polyfit(x,y,n)
x、 y表示数据点的横、纵坐标
n表示用来拟合的多项式的阶次
p为拟合的多项式，为n+1个系数的行向量

确定n的方法：

for i=1:8
	y2=polyfit(x0,y0,i);
	Y=polyval(y2,x0);%计算拟合函数在x处的值
	if sum((Y-y0).^2)<0.1 %误差平方和精度范围
		c=i；
		break;
	end
end

评估拟合效果的方法：
• 差值求和， sum(abs(nihe-shiji)) 取得最小即可；
• 和方差： sum((nihe-shiji)2) ，越接近0，效果越好

(2).非线性最小二乘拟合

Matlab提供的求非线性最小二乘拟合的函数有：
lsqcurvefit:
[x,resnorm]=lsqcurvefit(fun,x0,xdata,ydata))
lsqnonlin:
x=lsqnonlin(fun,x0,lb,ub)
nlinfit:
[beta,r,J]=nlinfit(x,y,fun,beta0)

• 先定义函数，通常用inline函数定义。
fun=inline(‘f(x)’,‘参变量’,‘x’)
• 然后利用nlinfit函数如下命令：
[beta,r,J]=nlinfit(x,y,fun,beta0)
x， y为原始数据， fun是定义好的函数， beta0是函数中参数的初始值
例子：

x=2:2:20;
y=[zeros(1,5) ones(1,5)];
plot(x,y,'o');
axis([0 20 -0.2 1.2])
myfun=inline('1./(1+exp(-a*(x-10)))','a','x')
a=5;
beta = nlinfit(x,y,myfun,a)
y1=1./(1+exp(-beta*(x-10)));
plot(x,y,'o',xp,y1,'r',xp,y5,'g',xp,ys,'k')
legend('data','一阶多项式拟合','五阶多项式拟合','s型曲线拟合')

在这里插入图片描述

% 已知某信号电平的测量值为：[0.9845, 0.7676, 
% 0.1897, -0.4542, -0.8871, -0.9896, -0.8076, -
% 0.4797, -0.1496, 0.0888, 0.2055 0.2189, 0.1699, 
% 0.0998, 0.0371, -0.0048, -0.0247, -0.0280, -0.0222, 
% -0.0135, -0.0057] mV，以第一个数据测量点作为时间零点，每隔0.5秒测量一次。用三次样条插值的方法完成如下操作：
% (1) 画出每隔1毫秒的插值曲线图形，并将原始的数据画在同一图上。
% (2) 估算6.6秒时的信号电平。
y0 = [0.9845,0.7676,0.1897,-0.4542,-0.8871,-0.9896,-0.8076,-0.4797,-0.1496,0.0888,0.2055,0.2189,0.1699,0.0988,0.0371,-0.0048,-0.0247,-0.0280,0.1699,0.0988,0.0371,-0.0048,-0.0247,-0.0280,-0.0222,-0.0135,-0.0057];
x0 = 0:0.5:10;
x1 = 0:0.001:10;
y1 = interp1(x0,y0,x1,'spline');
plot(x1,y1,'b-',x0,y0,'r-');
legend('spline','样本点')
%根据生成的数据进行插值处理，得出较为平缓的曲线
x0 = 0:0.12:1;
y0 = (x0.^2 - 3*x0 + 5).*exp(-5*x0).*sin(x0);
figure(1)
plot(x0,y0,'ro-')
x1 = 0:0.02:1;
y = (x1.^2 - 3*x1 + 5).*exp(-5*1).*sin(x1);
y1 = interp1(x0,y0,x1,'linear');
y2 = interp1(x0,y0,x1,'cubic');
y3 = interp1(x0,y0,x1,'spline');
y4 = interp1(x0,y0,x1,'nearest');
plot(x1,y1,'kx-',x1,y2,'rs-',x1,y3,'go-',x1,y4,'m^-',x0,y0,'yo-',x1,y,'b-')
legend('linear','cubic','spline','nearest','样本点','原函数')
figure(2)

在这里插入图片描述

% 已知平板表面分为5*3的网格，各处温度值为:
% [82 81 80 82 84;79 63 61 65 81;84 84 82 85 86], 做出平板表面温度分布曲面，在x,y方
% 向上每隔0.2个单位进行插值，画出插值后的图形。
[x,y] = meshgrid(0:4,0:2);
z = [82 81 80 82 84;79 63 61 65 81;84 84 82 85 86];
mesh(x,y,z) %mesh函数是绘制三维曲面的函数
[x1,y1] = meshgrid(0:0.2:4,0:0.2:2)
z1 = interp2(x,y,z,x1,y1,'cubic');
mesh(x1,y1,z1)
% 随机数据的二维插值
x=-3+6*rand(1,200);
y=-2+4*rand(1,200);
plot(x,y,'*')
z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);
z1=griddata(x,y,z,x1,y1,'cubic');
figure
surf(x1,y1,z1)

在这里插入图片描述

% 自变量x的取值为x=[129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5];
% 自变量y的取值为 y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5];
% 变量z的取值为Z=[4 8 6 8 6 8 8 9 9 8 8 9 4 9];
% 插值的取值范围：x为70:3:200; y为 -70:3:150；
x=[129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5];
y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5];
z=[4 8 6 8 6 8 8 9 9 8 8 9 4 9];
[x1,y1] = meshgrid(70:3:200,-70:3:150);
z1 = griddata(x,y,z,x1,y1,'cubic')
mesh(x1,y1,z1)

在这里插入图片描述

% 试用多项式拟合的方法在不同阶次下进行拟合，并观察效果，找到合适阶次。
x0=0:0.1:1;
y0=(x0.^2-3*x0+5).*exp(-5*x0).*sin(x0);
p3=polyfit(x0,y0,3);%p3为拟合多项式的系数
y3=polyval(p3,x0);
p4=polyfit(x0,y0,4);
y4=polyval(p4,x0);
p5=polyfit(x0,y0,5);
y5=polyval(p5,x0);
p8=polyfit(x0,y0,8);
y8=polyval(p8,x0);
p1=polyfit(x0,y0,1);
y1=polyval(p1,x0);
plot(x0,y3,'bo-',x0,y4,'g--',x0,y5,'m--',x0,y8,'b*-',x0,y0,'ro-')
legend('三次拟合','四次拟合','五次拟合','八次拟合','原函数曲线')

在这里插入图片描述

%某学生做材料实验，得到温度与长度的变化数据：温度 = 20:5:45; 长度= [81,82.3,84,86.9,89,90.8];
% 请用多项式拟合法，确定拟合阶数，确定L 和t的关系表达式(poly2str)，并计算和方差
x0 = 20:5:45;
y0 = [81,82.3,84,86.9,89,90.8];
for i=1:8
pi = polyfit(x0,y0,i);
yi = polyval(pi,x0)
if sum((yi-y0).^2)<1e-5
a = i
b = pi
c = sum((yi-y0).^2)
plot(x0,y0,'g--',x0,yi,'ro-')
legend('原函数','拟合函数')
break;
end
end
% x=1:20,y=x+3sin(x).
% 利用多项式拟合法，确定拟合阶数，并绘制出理论曲线与拟合曲线，对比拟合效果。```