matlab笔记
鉴于命令老是会忘记, 这里就分类对常用的matlab命令及相应的知识点进行一个总结了, 方便使用, 本文基于6大模块, 插值, 线代, 优化, 积分, 微分, 差分,主要内容是基于《高等应用数学问题的 Matlab 求解》这本书
插值
插值分类如下
- 包括一维, 二维, 多维插值, 基于插值求解数值积分和离散数据的最优化问题
- 常用的两种样条插值, 三次分段多项式插值, B样条插值, 基于样条插值的数值微积分运算, 其结果较1中更好
插值与数据礼盒
一维插值
y1 = interp1(x,y,x1,’spline’), 其中x,y为两组向量, x1为自定义的坐标, 除了’spline’,还包括’linear’,’nearest’,’cubic’, ‘pchip’, 一般建议用’spline’
对拉格朗日插值有公式, $${\phi}(x)=\sum^N_{i=1}yi\prod^N{j=1,j{\neq}i}\frac{x-x_j}{(x_i-x_j)}$$
所以可以写出lagrange插值函数为
function y = lagrange(x0,y0,x)
len = 1:length(x0);
y = zeros(size(x));
for i = len
ij = find(len~=i)
y1 = 1;
for j = 1:length(ij)
y1 = y1.*(x-x0(ij(j)));
end
y = y + y1*y0(i)/prod(x0(i)-x0(ij));
end
那么对于离散数据的定积分来说, 可以先进行插值后再积分, 结果更准确, 如
function y = quadspln(x0,y0,a,b)
f = @(x) interp1(x0,y0,x,'spline')
y = integral(f,a,b)
end
如此可以直接通过原来的值构造插值函数, 直接对函数进行积分
而对现成方法, 有trapz, 求解定积分提醒方法
二维插值
z1 = interp2(x0,y0,z0,x1,y1,’spline’)
这里的x0,y0,z0,一般用meshgrid生成, 需要为二维数组, 这个只能处理网格数据
对非网格数据, 有
z = griddata(x0,y0,z0,x,y,’v4’), ‘v4’是matlab提供的算法
对三维图像, 可以用 plot3,plot,surf, contour
注意这里和r一样, 可以用line进行直线添加
用find进行查找满足条件的索引
高维数据
1-3维数据可以用meshgrid生成
[x,y,z] = meshgrid(x1,y1,z1)
对n维数据 一般使用ndgrid
[x1,x2,..] = ndgrid(v1..vn)
两者的区别在于, 前者只能用于2-3维生成; 此外,两者生成数据的防线是相反的
所以高纬插值, 可以进一步用interp3或者interpn, 散则使用griddata
3维数据的可视化可以用slice, vol_visual4d
在最优化问题的利用
由于插值可以用于生成可能的函数, 那么将其函数作为最优化的目标函数, 就转换成了最优化问题
如
x = -3+6*rand(200,1);
y = -2+4*rand(200,1);
z = (x.^2-2*x).*exp(-x.^2-y.^2-x.*y);
f = @(p) griddata(x,y,z,p(1),p(2),'v4');
x = fminunc(f,[0,0]);
[x0,y0] = meshgrid(-3:0.1:2,-2:0.1:2);
z0 = (x0.^2-2*x0).exp(-x0.^2-y0.^2-x0.*y0);
contour(x0,y0,z0,30)
注意对最简单的无约束优化有fminunc,有约束用fmincon
样条插值
样条插值也是一种函数逼近的方法, 一般有三次样条插值和B样条插值,
三次样条插值
三次样条插值定义为, 满足一下三个条件
- 该函数必然经过样本点
- 在每个子区间上$[xi,x{i+1}]$为三次多项式
- 在整个区间$[x_1,x_n]$桑有连续的一阶与二阶导数
样条插值工具箱里有
S = csapi(x,y), 其中x,y就是样本点
S 为一三次样条函数对象
结果可以用fnplt画出来, 也可以fnval进行计算
fnplt(S), fnval(S,xp)
系数可以用如S.coefs(1:4,:)调出, 即S.breaks矩阵和S.coefs
x0 = [0,0.4,1,2,pi];
y0 = sin(x);
sp = csapi(x0,y0);
fnplt(sp,':');
hold on
ezplot('sin(t)',[0,pi]);
plot(x0,y0,'o')
注意这里ezplot的调用和mathematica一致
对样条插值,
对多变量的csapi,有
S = csapi({x1,x2,…,xn},z)
B样条插值
调用为spapi, S = spapi(k,x,y)
S = spapi({k1,k2},{x1,x2},y)
k为自定义的B样条阶次, 一般使用k=4,5, 提高阶次能改善插值效果
同样用fnval,fnplt进行绘制
基于样条插值的数值微积分
- 微分
可以用fnder函数, 对样条函数进行微分运算, 一般的, 对符号运算, 可以直接对符号用diff进行微分, 然后ezplot
syms x;
f = (x^2-3*x+5)*exp(-5*x)sin(x);
ezplot(diff(f),[0,1])
syms x y;
z = (x^2-2*x)*exp(-x^2-y^2-x*y);
ezsurf(diff(diff(z,x),y),[-3 3],[-2 2])
Sd = fnder(S,k) 求S的k阶导数
Sd = fnder(S,[k1,…,kn])
- 积分
可以用fnint
如a = fnint(s,1)可以返回积分函数,
对积分的值可以用
xx = fnval(a,[0,pi]);
xx(2)-xx(1);
即积分函数进行相加减
多项式拟合
直接调用polyfit
p = polyfit(x,y,n)
n为阶次, 返回多项式的系数, 之后可以用polyval进行求解, 求中可以用poly2sym转换成多项式形式, 这里可以把其看成符号表达式
如
syms x
p = [1,2,3,4,5,6,7];
a = poly2sym(p);
ezplot(a,[1,5])
可以绘制出图形
多项式拟合实际上相当于对已知函数用taylor级数进行展开, 但是由于函数是未知的
syms x
y = (x^2-3*x+5)*exp(-5*x)*sin(x);
p = taylor(y,'Order',9);
talor进行符号运算展开!, 这里可以选择展开的阶数
学会运用符号运算!
拓展-函数线性组合的曲线拟合方法
$$ g(x)=c_1f_1(x)+c_2f_2(x)+c_3f_3(x)+\ldots+c_nf_n(x)$$
对这个函数进行线性拟合, 在已经知道$f_n(x)$的情况下, 可以直接转换成最小二乘进行求解
$$Ac=y$$
那么我们可以把多项式拟合看成这次拟合的特例, $f_i(x)=x^{n+1-i},i=1,2,\ldots,n$, 所以可以直接用最小二乘对polyfit进行求解
x = [0:.1:2]';
y = (x.^2-3*x+5).*exp(-5*x).*sin(x);
n = 7;
A = zeros(length(x),n+1);
for i = 1:(n+1)
A(:,i) = x.^(n-i+1);
end
c = A\y
p = polyfit(x,y,n)
norm(c'-p)
最小二乘曲线拟合
[a,Jm] = lqscurvefit(Fun,a0,x,y,options)
f = @(a,x) a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x);
[xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y)
x1 = 0:0.01:10;
y1 = f(xx,x1)
plot(x1,y1,x1,y1'o')
可以看到的是返回拟合的参数,
这里需要注意的是需要提供样本点, 和参数的初始值和拟合函数, 一般用匿名函数进行, 对多变量的函数, 可以用x(:,1)代表, 传入数据x,y的时候 注意对应
特殊函数
- gamma函数, 可以看成广义阶乘, y = gamma(x)
- beta函数, y = beta(m,x), 有一个参数m
- bessel函数, 一类besselj(lambda,x), lambda作为阶次, 二类bessely, 三类besselh
- y = legendre(n,x)
- mittag-leffler可以看成指数函数的拓展, f = mittag_leffler(a,z)
信号分析与数字信号处理
相关性
对相关系数, 有corr自相关系数
xcorr可以给出自相关函数, cxx = xcorr(x,N)
cxy = xcorr(x,y,N), N表示相关点
stem用来画上点图,
滤波
线性滤波器
滤波器其实可以按照传函的形式进行书写, 所以a,b向量可以确定一个传递函数
- FIR滤波器
- IIR滤波器
- ARMA滤波器
matlab中可以调用y = filter(b,a,x)进行滤波,
一般滤波器分为低通, 高通和带通滤波器, 实际上就是放大器, 效果为, 在非通过频率放大倍数接近0, 在通过频率段放大倍数为接近1
如果滤波器已知道, 用freqz可以进行放大倍数分析
[h,w] = freqz(b,a,N)
综上, 可以用filter对一个线性滤波器进行构造, freqz对一个线性滤波器进行分析
滤波器设计
常用的有两种,
butter(n,wn)
cheby1(n,r,wn)
cheby2(n,r,wn)
n为阶次
wn为归一化频率,