数学实验与数学软件课件

数学实验与数学软件内容整理，参考教材：MATLAB教程

Lecture 1 MATLAB软件入门

1.1 MATLAB特点

MATLAB:MATrix LABoratory，是一种解释性程序语言，无需编译，大都按行逐句解释运行，直到发生错误或中断。

1.2 变量与特殊数值

... 是续行输入法

变量名区分大小写，首字母必须为字母，不可与关键词重复

realmax('classname') 返回某类型数据的实数最大值，同理 realmin

intmax('classname'='int32') 返回指定整数类型的最大值，同理 intmin

eps(2)等价于$2\epsilon$

1.3 运算符

A/B 和 A\B 区别： 2/3=0.66，3/2=1.5，在矩阵下定义为 A/B=A*inv(B) 和 A\B=inv(A)*B

A*B矩阵乘法，A.*B矩阵对应位置的乘法，同理A/B和A./B，矩阵乘除法需要满足行列关系否则报错

1.4 复数

i 以及 j 是预留的叙述单位，real() 取实部，imag() 取虚部，angle() 取幅角，可用 abs() 计算复数的模长，conj()计算复数共轭

以上函数均可直接在矩阵上使用，考虑向量化变成提高代码效率。

z = 4 + 3i
real_z = real(z)
image_z = imag(z)
magnitude_z = abs(z)
angle_z_radian = angle(z)
angle_z_degree = angle(z)*180/pi

绘图操作：

z1=4+3i;
z2=1+2i;
z12=z1+z2;
clf
hold on
plot([0,z1,z12],'-b','LineWidth',3)
plot([0,z12],'-r','LineWidth',3)
plot([z1,z12],'ob','MarkerSize',8)
hold off
grid on
axis equal
axis([0,6,0,6])
text(3.5,2.3,'z1')
text(5,4.5,'z2')
text(2.5,3.5,'z12')
xlabel('real')
ylabel('image')

复数相关计算，$\sqrt[3]{-8}$为例

a=-8;
r_a = a^(1/3)%计算主根
r_n = nthroot(a,3)%计算实根

p=[1,0,0,-a]; %p(r)=r^3-a的多项式系数
R=roots(p) %返回所有根

1.5 命令窗功能

clf 清空绘图窗口，clc 清空指令窗，clear 清除工作空间变量或clear name清除指定工作变量，close all 关闭所有弹出窗口

help name 显示帮助信息，doc name 在浏览器搜索指定内容， helpwin name 弹出窗口显示帮助信息， lookfor name 搜索词条

cd 改变目录，dir 当前目录下所有文件，exit 和 quit 退出MATLAB，diary 记录指令窗输入记录

Lecture 2 MATLAB编辑器与符号计算（一）

2.1 运算结果显示形式

format class 恢复预设值，默认为short

format short (e)/long (e) 5/15位有效数字（科学计数）

其它可选字段有rat近似有理数、hex十六进制、+大矩阵、bank金融，compact紧凑和loose稀疏

a=1/3;
% 默认显示精度为4，实际存储精度为双精度
a=0.3333

2.2 搜索路径原则

可以在设置搜索路径添加搜索路径，或者可以使用path或addpath命令

1. 首先检查内存（工作区），观察cont是否为变量
2. 若不是，检查cont是否为内建函数或常量
3. 在当前文件夹检查是否有名为cont的M文件，若有则需运行对应的M文件（可能为子函数或子程序）
4. 在搜索路径的其他文件夹查找是否有对应的M文件
5. 如果都没有，则会报错

2.3 符号对象

a = sym(0.3,'flag')则定义了一个符号变量，具体以3/10存储，flag参数指定存储类型，可选r有理形，d十进制浮点，e机器误差有理格式，f合并的有理形式

字符串小数转换为符号对象会被识别为32位精度的近似小数，计算会有误差。通常先定义符号对象再计算可以减小误差。

syms para1 para2 Flag则定义多个符号变量，同在Flag指定域下定义。可选real,positive等

syms F(x,y) 则定义符号函数，可以在后续定义具体形式。g(z)=u*w^2+z*w==v 可以定义符号恒等式函数

symvar(EXPR) 按字典序显示所有基本符号变量

symvar(EXPR,20) 按照ASCII码与’x’距离显示最近的20个

disp([‘XXXXX’,str,…])可依次输出多条字符串，非字符串可以用自带变量类型输出，也可用int2str转换，disp此时可等价换为display，多目标输出小括号内中括号不可省略

2.4 符号表达式

syms u v w z %先定义好基本符号变量
Eq=u*w^2+z*w-v %然后定义符号表达式
g(z)=u*w^2+z*w==v %定义符号恒等式函数

solve(Eq,var) 缺省变量则默认解的自由变量是x或离x最近的符号

vpa(symb,num)将符号对象或传统数化为有指定位数有效数字的结果

subs(symb,x,4)将符号对象中指定符号替换为新的对象

2.5 符号微积分

limit(f,x,a,'right') 计算f(x)在x=a处极限，x缺省则为默认变量，a缺省则默认为0，'right'也可选'left'，缺省则为两侧极限。

使用该函数计算左右导数。例如limit((subs(f_p,x,x+d)-f_p)/d,d,0)

g = diff(f,v,n) 计算f(x)关于v的n阶导数，v缺省则为默认变量，n缺省则为一阶

fjac=jacobian(f,v)对矩阵表达式求其雅可比矩阵，f与v需要定义成函数向量与自变量向量，若雅可比矩阵fjac为方阵，可使用detjac=det(fjac)来获得雅可比行列式的抽象表达式（符号型表达式）

r=taylor(g,x,a,'Order',10)可获得函数g在x=a点的9阶泰勒展开表达式，x缺省则为默认变量，a缺省则默认为0，不指定阶数则默认为5阶展开。

symsum(f,v,a,b) 对f(v)从a到b求级数和

Lecture 3 MATLAB符号计算（二）

$$(arctan(x))'=\frac{1}{1+x^2}$$ $$(arctanh(x))'=\frac{1}{1-x^2}$$ $$(arcsin(x))'=\frac{1}{\sqrt{1+x^2}}$$ $$(arcsinh(x))'=\frac{1}{\sqrt{1-x^2}}$$

3.1 符号对象的识别

size()返回矩阵行列数

class()显示变量类型

isa(matrix,'char')返回1或0，确定输入是否具有指定数据类型

whos para1 para2列表显示所有内存变量类别

3.2 符号转双精度或符号型小数

reset(symengine) 符号运算设置重置

vpa(x,n) 对x截取n位有效数字（仍是sym），n缺省则按当前digits位数

digit(n) 设置vpa截断位数n，可以用digit显示当前截断位数

eval(x) sym型转化成命令计算双精度结果

3.3 符号表达式的基本操作

collect(equation, variable) 合并同类项

expand(equation) 多项式展开

simplify(equation)简化函数符号表达式，可选参数step调整简化步数。使用时注意有时需要设置'IgnoreAnalyticConstraints',true 忽视复数根

[N,D]=numden(equation) (Numerator and denominator of a symbolic expression)通分为分子N和分母D

pretty(equation) 美化（并不美观）

[r,s]=subexpr(equation) 提取公子式（$r(s)$=equation）

% 特征向量分解
syms a b c d
A = [a b;c d];
[V,D] = eig(A);
[RVD,w] = subexpr([V;D])

subs(eqation,old,new) variable可以由tuple输入完成批量代换，对应返回高维数组，数据类型相同时可以由vector输入

int(function,variable,lower,upper) 符号函数计算定积分，上下限缺省时为不定积分

h1 = fplot(x(t),y(t),t_range) 例如t_range=[0,2pi]传入上下界的向量，h1是*绘图句柄，可以用set(h1,'Color','r','LineWidth',5)设置绘图具体参数

h1 = ezplot(y(1),[x_range,y_range])

3.4 微分方程组的符号解法

声明符号变量时需要指明是符号函数syms x y(x)

S=dsolve([eq1,eq2],condition)常微分方程理论解，返回解可能有常数C1,C2等。具体每一个解可以通过$S.name$得到。在对常数进行赋值时，可以借助symvar(S)来获得对应的常数符号。

3.5 傅里叶变换

进行时频分析，决定信号在什么频率上更强。

$$F(\omega)=\int_{-\infty}^{+\infty} f(t)e^{-j\omega t}dt$$

逆变换定义

$$f(t)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(\omega)e^{j\omega t}d\omega$$

Dirac function

在$\omega\not ={0}$处$F(\omega)\equiv 0$，而在$\omega=0$处“无穷大”

$$\delta(x)=\lim_{a\rightarrow 0^+}\frac{1}{a\sqrt{\pi}}e^{\frac{-x^2}{a^2}}$$

fourier(function,time,w) 时间域到频率域

ifourier(function,w,time) 频率域到时间域

傅里叶变换的一些性质

$$F(\omega)=\int_{-\infty}^{+\infty}f(t)e^{-j\omega t}dt$$

平移性质:

对于有平移混淆的信号（信号相似，有时间差），可以借助傅里叶变换平移性质，即模相等的性质取消平易混淆

$$|\mathscr{L}(f(t-t_0))|=|\mathscr{L}(f(t))|$$

导数性质:

$$\mathscr{L}(f'(t))=i\omega F(\omega)$$

其中$\mathscr{L}(f(t))=F(\omega)$

卷积性质:

傅里叶变换下对卷积计算的复杂性降低（从$O(n^4)$到$O(n^2logn)$），用于解决反卷积问题

Lecture 4 MATLAB符号计算（三）

4.1 MATLAB的符号限制性假设

在使用sym或syms定义符号变量时，MATLAB启动MuPAD引擎并提供内存工作空间执行符号运算。若不带限定性假设，默认为复数。

syms a positive 定义正实数变量

assume(a,'positive') 对已定义符号变量更改限制性假设，也可以添加表达式的限制性假设如assume(imag(x)==0)等形式

assumptions(x) 查看变量x当前的限定性假设

clear all 或 syms a clears 清除全部或符号变量a的限制性假设

assumeAlso() 对符号变量增加限制性假设（不覆盖原有假设）

4.2 Laplace变换

$$F(s)=\int_0^{\infty}f(t)e^{-st}dt$$

逆变换

$$f(t)=\frac{1}{2\pi j}\int_{c-j\infty}^{c+j\infty}F(s)e^{st}ds$$

laplace(ft,time,w) 从时间域到频率域

ilaplace(Fs,w,time) 从频率域到时间域

4.3 Z变换

$$F(z)=\sum_{n=0}^\infty f(n)z^{-n}$$

逆变换

$$f(n)=\frac{1}{2\pi j}\oint_\Gamma F(z)z^{n-1}dz$$

Z变换描述离散序列的复数频率域信息，可以把差分方程转化为代数方程。

ztrans(fn,time,z) 从时间域到频率域

iztrans(FZ,z,time) 从频率域到时间域

4.4 符号卷积

$$(u*h)(t)=\int_0^t u(\tau)h(t-\tau)d\tau$$

可以通过Laplace变换和逆变换的性质快速求出：

$$(u*h)(t)=\mathscr{L}^{-1}(\mathscr{L}(u(t))·\mathscr{L}(h(t)))$$

clear all
syms T t tao s
ut=exp(-t);ht=exp(-t/T)/T;
uh_tao=subs(ut,t,tao)*subs(ht,t,t-tao);
% 直接法
yt1=int(uh_tao,tao,0,t)
% 间接法
yt21=ilaplace(laplace(ut,t,s)*laplace(ht,t,s),s,t)
yt22=collect(yt21,(T-1))%合并同类项

4.5 符号矩阵分析与计算（数值矩阵）

det(A) 求行列式

diag(A) 取A对角元或利用向量形成对角矩阵

expm(A) 计算A的指数矩阵，需要将A特征值取指数复原

inv(A) 求逆矩阵

rank(A) 求A的秩

tril(A) 化A为下三角矩阵，也有triu(A)化为上三角矩阵

[V D] = eig(A) 计算A的特征值保存到D，特征向量保存到V

[V J] = jordan(A) 计算A的Jordan分解，AV=VJ

[U S V] = svd(A) 计算A的奇异值分解，A=USV’

[L U] = lu(A) 计算A的LU分解，A=LU

[Q R] = qr(A) 计算A的QR分解，A=QR

4.6 一般方程组符号解

[v1,v2,..] = solve(exp1,exp2,..,v1,v2,..) 对任何代数方程组的符号解

4.7 符号计算结果可视化

• ezcountour(fun) 画等位线
• ezcontourf(fun) 画填色等位线
• ezplot(fun,[xmin,xmax]) 二维曲线绘图，可接受参数方程：ezplot(funx,funy,[tmin,tmax])
• ezplot3(funx,funy,funz,[tmin,tmax]) 三维曲线绘图
• ezpolar(fun,[a,b]) 极坐标曲线
• ezsurf(funx,funy,funz,[smin,smax,tmin,tmax]) 空间曲面图，也可以对$z=z(x,y),a<x<b,c<y<d$的形式：ezsurf(z(x,y),[a,b,c,d])以及$z=z(x,y),a<x,y<b$的形式：ezsurf(z(x,y),[a,b])
• ezsurfc(funx,funy,funz,[smin,smax,tmin,tmax]) 带等位线的曲面图
• ezmesh(funx,funy,funz,[smin,smax,tmin,tmax]) 三维网格绘图函数，用法同ezsurf()
• ezmeshc(funx,funy,funz,[smin,smax,tmin,tmax]) 带等位线的网格图，用法同ezsurfc()

积分互补性质： Week4 PPT 39

Lecture 5 MATLAB数值数组与数组化编程

5.1 一维数组（向量）

A = start:step:end 创建行向量，while(start<end): 添加start, start+=step

linspace(start,end,n)创建一位数组，等价于start:(end-start)/(n-1):end

（步长为实数时，近似值向零取整，可以确保得到希望的步数，例如0:pi/4:pi一定能得到4个值）

logspace(a,b,n)创建等比数列，也即10.^linspace(a,b,n)，等价于10.^(a:(b-a)/(n-1):b)

ones(1,n)生成1行n列的全为1的行向量，one(n)生成元素全为1的n阶方阵

zeros(1,n), rand(1,n), randn(1,n)生成全为零，0-1均匀分布随机，标准正态分布随机的行向量

rng default使随机状态回复到Matlab初始值（伪随机）

A.'转置，conj(A)共轭，A'共轭转置

ndims(A)矩阵维数，size(A,nd)矩阵第nd维的规模（缺省nd则返回每一维的大小，类似shape）

length(A)矩阵最大维度规模（或向量长度）

numel(A)矩阵所含元素的总数目

5.2 矩阵（二维数组）

load A 读取当前文件夹下A.mat的矩阵

ones(m,n), zeros(m,n), rand(m,n), randn(m,n)常用zeros(m,n)申请存储空间，缺省一维则生成方阵

eye(n)生成单位方阵

magic(n)生成n阶幻方（行列和相等）

diag(A)A为方阵：提取方阵对角元为列向量 / A为1维向量：利用向量生成对角阵

random(gm,n), randsrc(m,n,[set;prob]), gallery(name)一般随机数字矩阵、指定字符集随机数组、特殊功能的测试数组

rng(seed) 设置随机种子

speye(n) 生成稀疏的n阶单位矩阵

B=sparse(A)正常定义的矩阵转化为稀疏矩阵（由三元组方式进行定义）

reshape(A,m,n)将矩阵A变形成m行n列的矩阵，matlab中矩阵存储可以理解为列优先的数组

repmat(A,m,n) 按行数复制m份，列数复制n份

flipud(A), fliplr(A) 上下、左右对称翻转

rot90(A,n) 矩阵A逆时针旋转n*90度

5.3 二维数组编址

A(r,c) 全下标编址，可以计算A(i,j)=A((j-1)*M+i), j=ceil(k/M), i=k-(j-1)*M, A(k)=A(i,j)，其中r和c可以接受向量。

计算时只需注意矩阵是列优先，按列存储的。

A(ind)访问ind包含的元素，得到规模与ind相同的向量，A(:)默认为列向量（ind类似定义为r=[2,3]，或L=A<=0逻辑方式索引）

注：A([1,end])是传入ind为[1,end]的两个单序号寻址，对应矩阵数组第一个和最后一个值

[rowsub, colsub] = ind2sub([M N],IND)得到IND对应的行列指标

IND = sub2ind([M N], rowsub, colsub)得到单序号IND

5.4 数组的运算

矩阵运算可以减少for循环，但使用时需要注意是矩阵乘除法还是数组乘除法。

bsxfun(@fun,A,B) 对矩阵A和B使用fun

all(vector) 返回0或1，是否全部非零

any(vector) 返回0或1，是否存在非零

注：对向量t用逻辑方法Lt = (t==0)会得到一个逻辑向量，常用于判断是否满足条件并进行定位操作

ceil() 向$+\infty$取整

fix() 向0取整

floor() 向$-\infty$取整

mod() 取模，rem() 取余，x和y的符号不同时，rem函数结果的符号和x的一样，而mod和y一样。这是由于这两个函数的生成机制不同，rem函数采用fix函数，而mod函数采用了floor函数

round() 四舍五入

sign() 返回1或-1或0，符号函数

[azimuth,elevation,r] = cart2sph(x,y,z) 直角坐标转球坐标，反向有sph2cart()

[theta,rho] = cart2pol(x,y) 直角坐标转极坐标，反向有pol2cart()

[theta,rho,z] = cart2pol(x,y,z) 直角坐标转柱坐标，反向有pol2cart()

Lecture 6 MATLAB矩阵函数与程序设计初步

6.1 矩阵的乘除运算

A*B要求满足矩阵乘法的行列要求

A\B=inv(A)*B或A/B=A*inv(B)可以解线性方程组$Ax=b$也就是x=A\b，矩阵方程$XA=B$可以用X=B/A来解出。超定方程（A列满秩）可以提供最小二乘解

关于数组运算，$$构成复数集合$\mathscr{C}^S$上的数组域

关于矩阵运算，$$构成复数集合$\mathscr{C}^{(M\times N)}$上的阿贝尔群

6.2 矩阵的幂运算与矩阵函数

若$A=Q·diag(\lambda_1,…,\lambda_n)·Q^{-1}$且特征值各异：

b^A为$Q·diag(b^{\lambda_1},..,b^{\lambda_n})·Q^{-1}$

A^b为$Q·diag(\lambda_1^{b},..,\lambda_n^{b})·Q^{-1}$

expm(A)为矩阵指数函数$Q·diag(e^{\lambda_1},…,e^{\lambda_n})·Q^{-1}$

logm(A)为矩阵对数函数$Q·diag(\ln \lambda_1,…,\ln \lambda_n)·Q^{-1}$

sqrtm(A)为矩阵平方根函数$Q·diag(\sqrt{\lambda_1},…,\sqrt{\lambda_n})·Q^{-1}$

funm(A,Hfun)为矩阵函数$Q·diag(f(\lambda_1),…,f(\lambda_n))·Q^{-1}$，Hfun为函数句柄

6.3 MATLAB控制流

if-elseif-else-end end不可缺省

switch-case value-otherwise-end value常用胞元数组

for-end以及while-end循环

break和continue同C++，死循环时可以在命令行用Ctrl+C停止运行

pause程序在该行暂停，可以用pause(n)暂停n秒。keyboard进入调试模式，键入dbcont恢复代码段的运行

a = input('input a number') 屏幕输出，并将输入值赋给a

return结束当前函数并返回上一级函数，若在主程序中则程序运行结束

a = inputdlg('input a number) 弹窗输入，返回值是胞元。另外有listdlg和questdlg进行列表和多选弹窗。

fminbnd(fx,a,b)可以计算某个区间内的一个极小值，对于波动较大的函数，极小值不一定是最小值

倍增区间极值搜索法，将分割区间的长度不断变小，从而可以获得“更多的极小值”，每一次分割都将获取最小的一个极小值。在单调区间个数有限的假设下，分割越细，就越有可能获得最小值。

最终当进一步分层后，获取的最小值在连续几次倍增后不变（或变化量极小），就退出循环，结束区间倍增，然后输出最终得到的最小值、最小值点以及最终倍增区间的个数。

6.4 胞元

B(1,2) 援引，返回内容标签 如[3x3 double]

B{1,2} 内容援引，返回具体内容

Lecture 7 MATLAB数值微积分

7.1 差分方法

dx = diff(X) 为向量或矩阵的向前差分运算，等价于dx = X(2:n,:) - X(1:n-1,:)，结果比X少一行

FX = gradient(F) 为向量或矩阵的中心差分运算，等价于FX(1)=F(2)-F(1),FX(end)=F(end)-F(end-1),FX(2:end-1)=0.5*(F(3:end)-F(1:end-2))，结果与F长度相同。矩阵的一元梯度仅包括横向

[FX,FY] = gradient(F)二元梯度运算，FX为横向梯度（每一行右减左），FY为纵向梯度（每一列下减上）

注：用realmin近似计算函数极限时，若分子分母都趋于$\epsilon$且阶数较高时，MATLAB不能正常返回结果而是会生成NaN。近似导数计算时步长太小会出现毛刺现象。

7.2 数值求和与累积和

sum(X,dim) 矩阵求和，dim默认为1即按列方向求和，返回行向量

Scs = cumsum(X) 矩阵按列方向求累积和，返回同型矩阵

7.3 数值积分

将定积分$\int_a^bf(x)dx$积分区间划分为$a=x_0<x_1<x_2<…<x_n=b$，对每一个节点分段$[x_i,x_{i+1}]$，分段积分$\int_{x_i}^{x_{i+1}}f(x)dx$的近似

中矩形公式

$$\int_{x_i}^{x_{i+1}}f(x)dx\approx f(\frac{x_i+x_{i+1}}{2})(x_{i+1}-x_i)$$

复合梯形公式

$$\int_{x_i}^{x_{i+1}}f(x)dx\approx \frac{f(x_i)+f(x_{i+1})}{2}(x_{i+1}-x_i)$$

trapz(y) 使用梯形公式方法对函数矩阵y按列方向求和，首尾元素求和时乘以权值0.5

d = pi/8; t = 0:d:pi/2; y=0.2+sin(t);
s_ta = d*trapz(y);

trap(x,y) 函数矩阵y按列方向求关于x的和，x为函数采样序列y对应的自变量列表，可以接受非均匀采样。

cumtrapz(y) 函数矩阵累积式按列方向求和，返回同型矩阵

cumtrapz(x,y) 函数矩阵累积式按列方向求和，x为指定自变量列表

Y = [1 4 9 16 25];
Q = cumtrapz(Y);
%Q = 1×5
%         0    2.5000    9.0000   21.5000   42.0000

辛普森公式

$$\int_{x_i}^{x_{i+1}}f(x)dx\approx \frac{f(x_i)+4f(\frac{x_i+x_{i+1}}{2})+f(x_{i+1})}{6}(x_{i+1}-x_i)$$

Monte Carlo法

$$\int_{x_i}^{x_{i+1}}f(x)dx\approx \frac{\sum_{j=1}^Kf(x_{i,j})}{K}(x_{i+1}-x_i)$$

其中$\forall 1\leq j\leq K,x_{i,j}\in [x_i,x_{i+1}]$为子区间随机采样点，一般用均匀分布。

Monte Carlo法除去计算数值积分外，也可以用于估计某些实数的值。例如$\sqrt{2}$可以通过在$X\sim U[0,2]$上多次取样本值，计算$x^2\leq 2$的占比，理论值为$\frac{\sqrt{2}}{2}$

7.4 匿名函数

F = @(var) sin(var).*var 一句话函数，F为函数句柄，括号内为参数

7.5 内建数值积分函数

Q = integral(fun,xmin,xmax) 匿名函数句柄以及积分上下限。可选参数'AbsTol'控制绝对误差（注意需要调小'RelTol'保证精度） 'ArrayValued'被积函数是否为阵列函数，'Waypoints'指定积分路径

Lecture 8 MATLAB微分方程数值解法

8.1 差分方程

利用固定间隔的函数值构成的递推关系称为差分方程。

$$f(x)=F(f(x-1),f(x-2),...,f(x-n))$$

常系数线性差分方程的形式如下

$$f(x)=a_1f(x-1)+a_2f(x-2)+...+a_nf(x-n)+g(x)$$

若$g(x)\equiv 0$则为常系数齐次线性差分方程，解法如下：

$$f(x)=a_1f(x-1)+a_2f(x-2)+...+a_nf(x-n)$$

解其特征方程$\lambda^n-a_1\lambda^{n-1}-a_2\lambda^{n-2}-…-a_{n-1}\lambda-a_n=0$

若得到n个两两不同的实特征根$\lambda_1<\lambda_2<…<\lambda_n$，则差分方程的通解为

$$f(x)=c_1\lambda_1^x+c_2\lambda_2^x+...+c_n\lambda_n^x$$

8.2 差分方程与微分方程

差分方程$f(x+d)=(1+d)·f(x)$等价于

$$f(x+d)-f(x)=d·f(x)$$

则当$d\leftarrow 0$时，问题等价于微分方程

$$f'(x)=\lim_{d\rightarrow 0}\frac{f(x+d)-f(x)}{d}=f(x)$$

结论： 小步长差分方程是微分方程的数值近似。

$$f(x+d)=(1+d)·f(x)=f(x)+d·f(x)$$

就是$f’(x)=f(x)$的一种简单的数值近似求解方法

8.3 欧拉法

假设微分方程为$y’(x)=f(x,y(x))$，初值条件$y(x_0)=y^*$

则从初值条件$y_0=f(x_0)=y^*$开始逐步迭代

$$y_{n+1}=y_n+d·y'(x_n) = y_n+d·f(x_n,y_n)$$

8.4 改进欧拉法

假设微分方程为$y’(x)=f(x,y(x))$，初值条件$y(x_0)=y^*$

改进的欧拉法对每一步迭代所导入的导数值进行预估和矫正：

$$\begin{cases} \overline{y_{n+1}}=y_n+d\cdot f(x_n,y_n)\\ y_{n+1}=y_n+\frac{d}{2}\cdot [f(x_n,y_n)+f(x_{n+1},\overline{y_{n+1}})] \end{cases}$$

改进欧拉法同样也可以解一阶的微分方程组，只需要将未知函与导函数看成向量$\vec{y}’(x)=\vec{f}(x,\vec{y}(x))=\vec{f}(x,y_1(x),y_2(x),…,y_n(x)),\vec{y}(x_0)=\vec{y^*}$

8.5 龙格-库塔法

假设微分方程为$y’(x)=f(x,y(x))$，初值条件$y(x_0)=y^*$

使用不同的增量位置代入并加权平均$\sum_{i=1}^nc_i=1$

$$y_{n+1}=y_n+d\cdot \sum_{i=1}^rc_if(x_n,\lambda_id,y(x_n+\lambda_id))$$

二阶龙格-库塔法的中点公式

$$y_{n+1}=y_n+d\cdot f(x_n+\frac{d}{2},y_n+\frac{d}{2}f(x_n,y_n))$$

8.6 一阶微分方程组演化实例（见PPT）

设矩阵$A$有实特征值$\lambda_1,\lambda_2$，并假设$A$有两个线性无关的特征向量$v_1,v_2$.过原点按$v_1,v_2$方向做直线，即构成了两条演化轨迹（trajectory）的渐近线。若特征值$\lambda>0$，则在对应的特征向量方向上趋于无穷；若特征值$\lambda=0$，则在对应方向上不作位移（静止不动或平行于另一条渐近线移动）；若特征值$\lambda<0$，则在对应方向趋于原点。

8.7 高阶微分方程

[t,y] = ode45(odefun,tspan,y0) 4阶龙格库塔法，t和y返回时间与函数序列，odefun为导函数句柄，tspan为数值解法的时间区间，y0为初值列向量

高阶微分方程可以转化成一阶微分方程组。

Lecture 9 MATLAB数值线性代数

9.1 矩阵的范数

norm(A) 计算矩阵（或向量）A的2-范数，也称为谱范数

$$||A||_2=\max_{||x||_2=1}||Ax||_2=\sqrt{\lambda_{\max}(A^TA)}$$

其中向量$x$的2-范数是欧式空间长度，$\lambda_{\max}$为$A^TA$的最大特征值，因此谱范数也等于$A$的最大奇异值，但不一定是最大特征值（$A^TA$半正定方阵，特征值非负）

norm(A,'fro') 计算矩阵$A$的Frobenius范数，即拉成长向量再计算2-向量范数，常用于矩阵误差分析

$$||A||_F=\sqrt{\sum_i(A^TA)_{i,i}}=\sqrt{trace(A^TA)}$$

norm(A,1) 计算矩阵（或向量）A的1-范数，也称为列范数

$$||A||_1=\max_{||x||_1=1}||Ax||_1=\max_{1\leq j\leq n}\sum_{i=1}^m|a_{ij}|$$

其中向量$x$的1-范数为向量$x$所有元素绝对值的和

norm(A,inf) 计算矩阵（或向量）A的$\infty$-范数，也称为行范数

$$||A||_{\infty}=\max_{||x||_\infty=1}||Ax||_\infty=\max_{1\leq i\leq m}\sum_{j=1}^n|a_{ij}|$$

“由向量p-范数诱导得到的矩阵p-范数”

sum(sum(A~=0.0)) 计算矩阵A的0-范数（不满足$||kA||=|k|\cdot||A||$故不是传统意义下的范数），含义为矩阵A的非零元素的个数

$$||A||_0=\# (A_{ij}\not ={0})$$

0-范数是对矩阵稀疏性的一种直接和准确的度量，可以用于最优化问题中求稀疏解（但是算法多数不收敛）：

$$\min\frac{1}{2}||Ax\cdot b||^2+\lambda||x||_0$$

[U S V] = svd(A); norm(diag(S),1) 计算矩阵A的核范数，定义为$A$所有奇异值的和。对核范数的约束是低秩逼近与低秩求解最常用的方法。（对奇异值的0-范数换位1-范数）

$$||A||_*=\sum_{i=1}^n\sqrt{\lambda_i(A^TA)}=trace(A)$$

9.2 矩阵变换与特征值分解

[R,ci]=rref(A) 实现初等变换将$A$化为阶梯型矩阵，其中R存储了行阶梯型的矩阵，ci为阶梯元素（线性独立）的列指标

null(A) 得到$A$的核空间的一组标准正交基

orth(A) 得到$A$的像空间的一组标准正交基

[VR,DR]=cdf2rdf(V,D) 将复数对角型转换为实数块对角型(VR*DR/VR=A，了解)

9.3 MATLAB解线性方程组

A\b=inv(A)*B 可解得$x=A^{-1}b$的解，反斜杠效率高于inv()，无解问题或超定问题提供最小二乘解，不满秩且有解则随机提供一个解

对于$A\cdot x=b$

• 线性正问题：已知$A,x$求$b$，复杂度为$O(n^2)$
• 线性反问题：已知$A,b$求$x$，复杂度为$O(n^3)$

9.3.1 数值迭代法

从线性方程组$A\cdot x=b$出发，令$A=M-N$，则

$$x=M^{-1}Nx+M^{-1}b=Bx+f$$

其中$B=M^{-1}N,f=M^{-1}b$

则可以由迭代法求解：

$$x^{(k+1)}=Bx^{(k)}+f$$

这种算法在$||B||<1$的情况下可以保证收敛性

9.3.2 雅克比迭代法

采用矩阵分裂记号，令$M=D,N=L+U$

$$A=D-L-U$$

D=diag(A);L=-tril(A)+D;U=-triu(A)+D

则有

$$x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b$$

9.3.3 高斯-赛德尔迭代法

令$M=D-L,N=U$则有

$$x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b$$

9.3.4 基于变分法解线性方程组

若矩阵$A$对称正定，定义

$$\phi(x)=\frac{1}{2}<Ax,x>-<b,x>=\frac{1}{2}x^TAx-b^Tx$$

则线性方程组$Ax=b$的解集等于$\underset{x\in \mathbb{R}^n}{\argmin}\phi(x)$的解集，将线性方程组求解问题转化为最优化问题的方法称为变分法，当$A$正定时，最优化问题是一个凸问题，仅有一个全局最优解

9.3.5 变分-最速梯度下降法法

光滑函数$\phi(x)$在$x=x^{(k)}$下降最快的方向为其负梯度方向

$$-\nabla\phi(x^{(k)})=-(Ax^{(k)}-b)=p^{(k)}$$

则只需通过迭代$x^{(k+1)}=x^{(k)}+\lambda p^{(k)}$直到收敛。最速下降法定义为

$$\lambda^{(k)}=\underset{\lambda}{\argmin}\phi(x^{(k)}+\lambda p^{(k)})$$

解得最速步长$\lambda_k=\frac{}{}$

9.3.5 变分-共轭梯度法（了解）

与梯度下降法的主要区别在于每一步的移动方向$p^{(k)}$会在$span\left\{p^{(0)},…,p^{(k-1)}\right\}$中进行最优化的选取

x=pcg(A,b) Preconditioned conjugate gradients 实现带有简单病态预处理的共轭梯度法求解$Ax=b$的问题，矩阵$A$需要是对称正定的矩阵

Lecture 10 MATLAB进阶程序设计与问题求解

Recap:

M脚本：代码的组合，一次运行可以将一串代码按照顺序依次运行

M函数：单个返回值可以省略[]，只有一个函数的M函数可以不写end，函数申明行和函数体必须包含在内

edit function_name 查看函数源码

10.1 多种输入/输出格式的函数设计

H1后，空行表示分级，help只显示H1级的内容

function [S,L]=exm060201(N,R,str)
% 下面的一行称为H1行，接下来的%层均help可见
% exm060201.m    The area and perimeter of a regular polygon (正多边形的面积和周长)
%       N       The number of sides
%       R       The circumradius
%       str     A line specification to determine line type/color
%       S       The area of the regular polygon
%       L       The perimeter of the regular polygon
% exm060201    %所有可行的用法
% exm060201(N)
% exm060201(N,R)
% exm060201(N,R,str)
% S=exm060201(...)   %返回值如果少一个，默认保留第一个
% [S,L]=exm060201(...)

%  Zhang Zhiyong 修改于 2013-12-20

nargin为输入参数列表的个数，nargout为输出变量个数

switch nargin%检查输入变量（参数列表的个数，如果不同，补值或报错）
    case 0,N=100;R=1;str='-b';
    case 1,R=1;str='-b';
    case 2,str='-b';
    case 3,;%不进行任何操作并跳出switch分支,分号可省略
    otherwise,error('输入量太多。');%error意味弹出错误提示并中止程序运行
end;%分号可省略
t=0:2*pi/N:2*pi;x=R*sin(t);y=R*cos(t);
if nargout==0%无输出变量则只画空心圆图
    plot(x,y,str);
elseif nargout>2%输出变量过多则报错
    error('输出量太多。');
else %返回值如果少一个，默认保留第一个
    S=N*R*R*sin(2*pi/N)/2;
    L=2*N*R*sin(pi/N);
    fill(x,y,str)       %画填色图
end

10.2 MATLAB主函数与子函数

一个M文件中，可能会有多个函数，其中第一个称为主函数（与m文件名相同），后面的所有函数称为子函数。所有的子函数都可以被M文件内的脚本或主函数调用，但无法被其他M文件或命令行直接调用。因此，子函数是一种减少M文件数量，封装外部脚本不直接调用的函数的好方法。

同名函数调用顺序：子函数->MATLAB内建函数->其它M文件主函数

10.3 函数句柄

fun_handle=@fun 创建函数句柄，用于函数传参，或作函数返回值（返回封装的子函数的句柄）

fun=@(var1,var2) expression 创建匿名函数，并赋值函数句柄

函数句柄类似指针，既可以通过匿名函数定义，也可以由内建函数或自定义函数创建。调用子函数时，函数句柄可以“完全地代替”本身子函数的函数名，格式与子函数直接调用相同。但利用函数句柄多次调用可以大大节省调用的时间。

which('fun_name') 检查文件夹内有无函数定义

isa(A,dataType) 确定输入是否具有指定数据类型

10.4 函数极值问题

[x,fval,exitflag]=fminbnd(fun,x1,x2) 求一元函数fun在[x1,x2]的一个极小值，fun可以用字符串，匿名函数或函数句柄。返回值列表x为极小值点，fval为函数极小值，exitflag>0代表此函数成功找到了一个极值点。

[x,fval,exitflag]=fminsearch(fun,x0) 利用无导数方法（单纯形法）从x0（由列向量可定义为多个不同起点）出发，求多元函数fun在多维空间的一个极小值，fun建议使用多元匿名函数或多元函数句柄（输入值为向量）。返回值列表x为极小值点向量，fval为函数极小值

ff=@(x) 100*(x(2)-x(1)^2)^2 多元匿名函数的一个定义示例，即输入为向量

注意上面两个函数都只能求局部最优解，不一定能够得到全局最优解。

fminunc(ff,x0,options) 多元函数局部无约束最小值函数（拟牛顿法），仅接受一种初值计算

10.5 单峰函数

单峰函数 $\argmin_{x\in [a,b]}f(x)=\bar{x},\forall x\in [a,b]$，且$f(x)$在$x<\bar{x}$单调递减，在$x>\bar{x}$单调递增

单峰函数先减后增（可以不连续不可导），凸函数是单峰函数，而单峰函数不一定是凸函数（羊角型）。

从单峰函数极值点的两侧点$[a_0,b_0]$不断向内搜索，并且不断缩小左右端点的距离的方法被称为试探法（或三分法）

10.5.1 黄金分割法（要求）

对于向内的试探法，从$[a_k,b_k]$开始，当$a_1\approx 0.618a_k+0.382b_k,a_2\approx 0.382a_k+0.618b_k$时，探索的效率最高，这样的试探法就成为黄金分割法（0.618法）

例如寻找最小值时，两个新节点选择函数值更大一个的和另一侧端点（想象一下画面理解：同一侧坡面和不同侧坡面的情形）

10.5.2 一元问题牛顿迭代法（要求）

原理是一阶泰勒展开近似$f(x)\approx f(x_0)+f’(x_0)(x-x_0)$

若$\exist \bar{x},f(\bar{x})=0$，给定任意一点$x$有

$$f(x)\approx f(\bar{x})+f'(\bar{x})(x-\bar{x})=f'(\bar{x})(x-\bar{x})$$

移项后得到$\bar{x}\approx x-\frac{f(x)}{f’(x)}$，通过下式迭代可以计算$f(x)=0$的解：

$$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}$$

对于求局部最小值问题，因为$\argmin_x f(x)\Leftrightarrow f’(x)=0$，因此对应的牛顿迭代法只需改为导数方程问题：

$$x_{n+1}=x_n-\frac{f'(x_n)}{f''(x_n)}$$

10.5.3 最速梯度下降法

将一元最优化问题推广到多元问题$\underset{x\in\mathbb{R}^n}{\argmin}\phi(x)$，其中$\phi(x)$满足一定光滑性且其梯度除极值点外不为0向量。对于凸问题，只需迭代计算下降最快方向（负梯度）计算即可接近理论最小值：

$$x^{(k+1)}=x^{(k)}+\lambda p^{(k)}$$

多元问题下的最速下降法，最佳步长计算较困难，常用固定步长代替。判断多元问题是否为凸问题，即是要计算$\phi(x)$对应的Hessian矩阵是否半正定

10.5.4 牛顿法解多元最小值问题（了解）

定义Hessian矩阵

$$\nabla^2f(x,y)= \left[ \begin{matrix} \frac{\partial^2f}{\partial x^2} & \frac{\partial^2f}{\partial x\partial y}\\ \frac{\partial^2f}{\partial x\partial y} & \frac{\partial^2f}{\partial y^2} \end{matrix} \right]$$

多元问题在使用最速梯度下降法求解时，由于相邻迭代的梯度方向正交，移动方向呈现锯齿形，收敛效果差，对于病态矩阵无健壮性。因此用下式迭代（Hessian必须可逆）：

$$x^{(k+1)}=x^{(k)}-(\nabla^2f(x,y))^{-1}\nabla f(x,y)|_{x^{(k)}}$$

10.6 MATLAB解非线性方程（组）

[x,fval] = fzero(fun,x0) 以x0为初值，尝试寻找函数句柄或匿名函数的一个零点，fval为对应函数值。此解是fun(x)变号的位置，搜索算法并不能保证得到最近的零点。

[x,fval]=fsolve(fun,x0) 以x0为初值(向量)，尝试寻找函数向量fun的一个零点，fval为对应函数向量值

Lecture 11 MATLAB概率统计与曲线拟合

11.1 常见分布的概率和密度函数

11.1.1 二项分布

$$P(X=k)=\binom{N}{k}p^kq^{(N-k)},k=0,1,...,N$$

即独立重复N次伯努利试验后，总发生次数为k的概率

pk=binopdf(k,N,p) 即$P(X=k)$

dk=binocdf(k,N,p) 即$P(X\leq k)$

R=binornd(N,p,m,n) 生成元素服从二项分布$Bin(N,p)$的规模mxn的随机矩阵

11.1.2 正态分布

服从正态分布$N(\mu,\sigma^2)$的概率分布函数为

$$\phi(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{\frac{-(x-\mu)^2}{2\sigma^2}\right\},x\in (-\infty,+\infty)$$

对应的密度函数为$\Phi(x)=\int_{-\infty}^x \phi(t)dt$

px=normpdf(x,Mu,Sigma) 即$\phi(x)$

px=normcdf(x,Mu,Sigma) 即$\Phi(x)$

R=normrnd(Mu,Sigma,m,n) 生成元素服从正态分布$N(\mu,\sigma^2)$的规模mxn的随机矩阵

R=randn(m,n) 生成元素服从标准正态分布的规模mxn的随机矩阵

11.1.3 其它分布

pd=makedist('Poisson',lambda) 建立一个特定的概率分布(ProbDist)对象，后可调用pdf(pd,0:m)获取$k=0,1,…,m$的概率值或cdf(pd,0:m)同理

11.2 随机流的控制

rng default 初始默认的随机流状态（伪随机）

rng shuffle 根据时间变量设置随机流（真随机）

rnd(seed) 以非负整数seed(0~65535)为种子生成随机数

rng(sd,generator) 以特定的生成器生成随机数

st=rng 获得随机生成器构架(struct)并可用rng(st)再现

rand(m,n) 生成(0,1)间均匀分布的随机矩阵 randn(m.n)生成标准正态分布的随机矩阵

randi([imin,imax],m,n) 生成[imin,imax]内所有整数均匀分布的随机矩阵

11.3 统计分析命令

min(X);max(X) 计算矩阵各列的最大值或最小值，若要计算整个矩阵的最值可以用max(max(X))或max(X(:))

mean(X);median(X) 计算矩阵各列的均值和中位数

std(X);var(X) 计算矩阵各列的样本标准差和样本方差即ddof=1，若要ddof=0可以通过std(X,1)得到

cov(X) 计算矩阵各列所组成列向量计算出的协方差矩阵

corrcoef(A) 计算矩阵A各列所组成列向量的相关系数，或corrcoef(A,B)计算两个向量之间的相关矩阵

11.4 多项式拟合

$$y=a_1x^n+a_2x^{n-1}+...+a_nx+a_{n+1}$$

利用函数曲线上采样点$(x_i,y_i)$拟合确定系数$a_i$的问题，一般用逼近或插值方法拟合

p=polyfit(x,y,n) 通过数组x和y的数据进行拟合，阶数为n，返回多项式系数p（降幂）

yy=polyval(p,x) 多项式系数回代，一般用2-范数（平分残差），1-范数（绝对值残差）或无穷范数（一致逼近残差）分析拟合情况

polyfit函数的方法也即解最小二乘问题

$$\min_{a_i}||y-(a_1x^n+...+a_nx+a_{n+1})||_2$$

构造设计阵$X$即可将问题写作$Y=Xa$，其解为a=X\Y，注意一般设定维数小于样本个数

Ridge Regression 2-范数惩罚系数，使参数尽可能小

$$\min_{a_i}||y-(a_1x^n+...+a_nx+a_{n+1})||_2^2+\lambda ||a||_2^2$$

LASSO 1-范数惩罚，易得稀疏解

$$\min_{a_i}\frac{1}{2}||y-(a_1x^n+...+a_nx+a_{n+1})||_2^2+\lambda ||a||_1$$

LAR 最小绝对残差模型，对离群值处理更优

$$\min_{a_i}||y-(a_1x^n+...+a_nx+a_{n+1})||_1$$

11.5 cftool工具箱的使用

权值的选择，相当于解$\underset{a}{\min}||W(Xa-b)||_2$

傅里叶（三角）函数拟合 $y\approx a_0+\sum_k a_kcos(wx)+b_ksin(wx)$

高斯函数拟合 $y\approx \sum_k a_k \exp\left\{-(\frac{x-b_k}{c_k})^2\right\}$

有理式拟合 $y\approx\frac{p_1x^m+…+p_mx+p_{m+1}}{x^n+…+q_{n-1}x+q_{n}}$

非参数拟合 Cubic Spline

对于噪声，可以人工去除离群点，或用smooth五点均值滤波，也可以使用工具箱的Robust方法，如LAR最小绝对残差

Lecture 12 MATLAB数字信号与声音处理

12.1 函数信息与数字信号

数字信号是离散定义的，如声音信号有每秒若干次的采样，强度也有一个固定精度的取值。数字化信号一般需要抽样、量化和编码。

波形信号通过拟合、插值，往往代表这一种具有光滑性与周期性的函数。

声音信号则往往有多种频段的信息，对应的连续傅里叶变换在部分区段或位点出现更大的系数

一维数字信号还可能是分片光滑函数或分片常函数，也可能是符合某种特定分布的随机函数

含噪信号定义为$y=f(x)+\epsilon$，噪声有多种，如高斯加性噪声则满足$\epsilon\sim N(0,\sigma^2)$

平均平方误差(MSE) Mean Square Error

$$\sum[y-f(x)]^2$$

信噪比 Signal to Noise Ratio

$$10\cdot\log_{10}(\frac{\sum f^2(x)}{\sum[y-f(x)]^2}) (dB)$$

snr(s,n) s为信号，n为误差（也即y-s）

12.2 均值滤波

对于一些理想的0均值噪声，其幅度明显小于信号本身的强度（SNR>10），故采用均值方法对带噪信号进行光滑化处理就可以得到近似的恢复

f1=smooth(y)五点均值滤波(MA5)，返回的是列向量，或f1=conv(y,ones(1,5)*1/5, 'same')

f2=conv(y,[1/16,1/4,3/8,1/4,1/16], 'same')五点加权滤波

更大的滤波范围会导致图像过度平滑

12.3 离散傅里叶变换

一维离散傅里叶变换定义

$$Y(k)=\sum_{i=1}^nX(i)e^{\frac{(-2\pi j)(i-1)(k-1)}{n}},1\leq k\leq n$$

特别地，当$k=1$时，$Y(1)=\sum_{i=1}^nX(i)$称为离散傅里叶变换的低频系数，其余系数称为高频傅里叶系数。低频系数表达了函数的总和，高频系数在不同程度刻画函数的变化情况。

离散傅里叶变换实际上是将序列$X(i)$假想为以n为周期的无限长周期序列，进而可以利用有限个特定频率采样完成其频率表示

fft(X) 通过上式计算向量X的离散傅里叶变换

离散傅里叶逆变换定义

$$X(i)=\frac{1}{n}\sum_{k=1}^nY(k)e^{\frac{(2\pi j)(i-1)(k-1)}{n}},1\leq i\leq n$$

ifft(Y)通过上式计算离散傅里叶逆变换

离散傅里叶变换可以刻画傅里叶级数的系数？正交基只对某维度影响？例题$f(x)=sin(x)$傅里叶系数是稀疏的，在中间部分可以直接设为0达到去噪的效果。

12.4 硬阈值和软阈值

硬阈值算子

$$T_\lambda(y)= \begin{cases} y,|y|\geq\lambda\\ 0,|y|<\lambda \end{cases}$$

对y的傅里叶系数做硬阈值算子处理相当于求解问题$\underset{x}{\min}\frac{1}{2}||x-y||_2^2+\frac{1}{2}\lambda^2||\hat{x}||_0$

软阈值算子

$$\Tau_\lambda(y)=sign(y)\cdot\max(|y|-\lambda,0)$$

对y的傅里叶系数做软阈值算子处理相当于求解问题$\underset{x}{\min}\frac{1}{2}||x-y||_2^2+\frac{1}{2}\lambda^2||\hat{x}||_1$

硬阈值算子对分片光滑函数效果较差（只处理了判定域内的部分），软阈值方法更好，但仍然不如五点均值，是因为函数不具备全局光滑性或周期性，导致傅里叶系数的稀疏度不足

12.5 数字信号去噪问题

去噪方法通常基于噪声类型和强度以及真实信号值的函数特性。常见的噪声除高斯噪声外，还有冲击噪声（随机的部分采样值出现与真实信号无关的错误函数值），可以通过中值滤波器进行判断和调整

泊松噪声是一种比较难处理的非加性噪声，在每点设$\lambda=k\cdot real$并基于参数为$\lambda$的泊松分布下获取一个整数值$Y$，带噪信号对应的信号值为$\frac{Y}{k}$

傅里叶变换对于三角函数去噪效果极佳，因为三角函数在此变换系数满足完美稀疏，对于分片光滑函数效果则一般，根据目标函数特性选择合适变换及约束（正则化）方法很关键

12.6 声音信号处理

info=audioinfo('filepath') 获取声音文件的主要特性

[y,fs]=audioread('filepath') 读取声音文件的数字信号和采样率。采样率即一秒内有多少个采样点组成信号，y常有两列分别对应左右声道，取值在[-1,1]之间

ys=y(1:k*fs,1) 截取前k秒的左声道信号

sound(ys,fs) 按采样率fs播放ys声音信号

audiowrite('filename',ys,fs) 将声音信号ys按采样率fs写入文件

clear sound 结束播放

响度、振幅、音色等等，乐理中A4定义为440Hz，半音之间频率差异为$2^{1/12}$倍，一般可以通过短时傅里叶变化最大模所在的频率对音乐进行音调分析

ub = length(ys)-4410; % 窗口个数
height = zeros(ub,1);
for i = 1:ub
    temp = ys(i:i+4410);
    tempf = fft(temp);
    [maxv, maxp]=max(abs(tempf));
    height(i) = maxp; % 最大模的频率
end
plot(1:ub,height);

音叉音生成 使用特定的包络线函数控制振幅变化

mod = sin(t/44100*pi) 半周期三角函数包络线，用mod.*ys即可控制振幅

Lecture 13 MATLAB多项式运算与数据可视化（一）

13.1 多项式表示和运算

降幂排序$a(x)=a_1x^n+a_2x^{n-1}+…+a_nx+a_{n+1}$存储为系数向量a=[a1,a2,...,an+1]

多项式乘法和离散卷积的联系：

$(x^2+1)(x-1)=x^3-x^2+x-1$相当于计算系数的离散卷积$[1,0,1]*[1,-1]=[1,-1,1,-1]$

一种卷积计算方式的理解（错车）

1. 反转向量
2. 另一个向量逐格移动，上下相乘并加和得到卷积结果

如$[1,0,2]*[1,-1]$对应$[2,0,1]$和$[1,-1]$的“错车”，求得$[1,-1,2,-2]$

c=conv(a,b) 获取多项式相乘系数

[q,r]=deconv(b,a) 多项式带余除法，q为商r为余式

[r,p,k]=residue(b,a) 即为$\frac{b(x)}{a(x)}=\sum\frac{r_i}{x-p_i}+k(x)$，其中r为留数向量，p为极点向量，k为余式系数行向量

r=roots(a) 求多项式$a(x)$的所有根（含重根和复根），返回列向量

a=poly(r) 当r为行向量时返回根为r的多项式；当r为方阵时返回r的特征多项式（用于计算特征值）

V=polyval(p,X) p为多项式的稀疏行向量，X为任意规模矩阵，将X每一个元素代入并获得数值矩阵p(X)

V=polyvalm(p,X) 以矩阵X为整体的多项式代入结果的矩阵

poly2str(a,'s') 将以a向量元素为系数，s为自变量生成多项式字符串，按降幂排列

快速构造有形如$p(x)=x^3+ax^2+bx+c$特征多项式的方阵：$\begin{pmatrix}-a&-b&-c\\1&0&0\\0&1&0\end{pmatrix}$

diag(v,k) 将向量v的元素放置在第k条对角线上，k>0时为上方

13.2 有限长向量卷积

已知$a=[a_1,a_2,…,a_n],b=[b_1,b_2,…,b_m]$

$$(a*b)_k=\sum_{i=1}^{k-1}a_ib_{k-i},2\leq k\leq m+n$$

这里补充定义$\forall i>n,a_i=0,\forall j>m,b_j=0$

conv(a,b,'same') 返回长度与a相同的卷积结果，即取卷积结果中间n个元素的结果（默认选后面的）

conv(a,b,'valid') 返回长度等于n-m-1的卷积结果，即仅包含b的全部元素参与卷积的中间项

[q,r] = deconv(u,v) 反卷积函数，即u = conv(v,q)+r

13.3 离散数据的可视化

hist(x,nbins) 直方图绘制

stem(X,Y) 将数据序列Y绘制为从沿x轴的基线延伸的针状图

plot(x,y,'r.','Name','Value') 以x为横坐标序列，y为纵坐标序列，’s’表示线性与线色的字符串；允许X或Y为矩阵，对应绘出多条曲线；允许x缺省，以y的下标作为x绘图；允许一句命令内绘制多条曲线

线型：-实线；:虚点线；--虚划线；o空心圆圈；.实心点；+十字符；*星字符

线色：b蓝色；g绿色；r红色；c青色；m品红；y黄色；k黑色；w白色

可控的调用格式：MarkerSize点大小；Color颜色；LineStyle线型；LineWidth线宽；Marker点形状；MarkerEdgeColor点边界色彩；MarkerFaceColor点域色彩

坐标轴控制： axis auto

equal横纵轴等长刻度；off取消轴背景；image等长刻度且坐标框紧贴数据周围；on使用轴背景

grid on/off 开启或关闭网格线

hold on/off 开启/关闭多层叠绘

box on/off 右侧和上侧是否有实刻线包裹

title(S);xlabel(S);ylabel(S) 标题，x轴和y轴标签名

legend(S1,S2,...) 绘图图例字段名

text(xt,yt,S) 在(xt,yt)坐标插入文本

xticks;yticks 设置横坐标与纵坐标刻度线位置或刻度间距

字符串设置：\fontname{Courier New}设置字体；\fontsize{14}设置字体大小；\bf粗体；\it斜体；\sl另外一种斜体；\rm罗马？

Latex标识（略）

plotyy(x,y,x,s,'stem','plot') 双纵坐标视图

Lecture 14 MATLAB数字图像处理初步

14.1 数字图像与矩阵

灰度图像对应一个二维矩阵$f(x,y),1\leq x\leq M, 1\leq y\leq N$，其中(M,N)为图像尺寸，而$f(x,y)$为图像的灰度。坐标类似矩阵定义，如左下(256,1),右上(1,256)

彩色图像存储方式为RGB存储模式，红绿蓝分别对应一种亮度（灰度）值，记为$f(x,y,z),1\leq x\leq M, 1\leq y\leq N,z\in \left\{1,2,3\right\}$

info = imfinfo('filename') 获取灰度或彩色图像文件的基本信息

A = imread('filename') 获取灰度或彩色图像的信息并存储在矩阵中

imshow(A) 弹出Figure并显示图像A，默认最小值为0，最大值为1，0~255的double图像可以标准化图示如imshow(A/255)，或者指定imshow(A,[low,high])

数字图像的默认存储格式为uint8，但在进行数值计算时会被自动转化成double，而在图示前矩阵强制转换为uint8。

14.2 彩色图像通道分离与图像存储

RGB格式彩色图像矩阵第三维顺序对应红色、绿色、蓝色通道的灰度。

C = rgb2gray(A) 根据彩色图像A的整体亮度均匀转化为灰度图像C

imwrite(A,'filename') 将任意灰度或彩色图像矩阵存入对应文件，注意double型矩阵存储时范围为0~1

14.3 彩色图像的颜色编码更换

色调(Hue),饱和度(Saturation),明度(Value)模型HSV是另一种存储方式。

B = rgb2hsv(A) 将RGB模型的矩阵A化为HSV型的矩阵B

C = hsv2rgb(B) 将HSV模型的矩阵转化为RGB模型的矩阵

若A为uint8型的数据就会自动转为double并除以255，再进行对应的色彩空间变换。HSV的指标范围为0~1

14.4 数字图像的尺寸处理

imresize(RGB, [numrows,numcols]) 若numrows或numcols设为NaN则为等比例放大或缩小。默认使用双三次差值的方法放大图像

RGB(rows,cols,:) 按rows和cols的索引位置截取图像。可以使用start:-1:end实现翻转的效果，或start:2:end实现下采样

[x,y,z] = sphere(100)获取单位球面坐标，warp(x,y,z,RGB) displays the image on the surface (X,Y,Z)

14.5 数字图像的亮度和对比度

若int8型的灰度图像二维矩阵A，设置常整数$c\in(-255,255)$，用A+c完成亮度调整。A+c超过255时自动设置为255，小于0则自动设置为0

设uint8型的灰度图像二维矩阵A，设置正常数$c>0$和0~255之间的整数k，则$c\times(A-k)+k$表示为以k为灰度中心，对比放大/缩小c倍。（通常对比度定义为最大最小灰度比值，k=0），结果仍为uint8型数据

亮度和对比度类似“均值”移动和“方差”改变

14.6 图像的直方图显示与均衡化

imhist(A) 绘制图像的灰度直方图，统计各灰度出现频率

B = histeq(A) 矩阵A进行灰度均衡化（增大对比度）

14.7 图像背景提取与计算（了解）

imopen 基于图像的开运算算法，核心是对图像的腐蚀再膨胀的运算（提取邻域内最暗点并代替领域，提取邻域内最亮点并代替领域）

imsubstract 图像相减

imabsdiff 前景值不论非负都保留

immultiply 图像乘法

14.8 二维离散傅里叶变换

二维离散傅里叶变换可以理解为一维离散傅里叶变换的张量积，或理解按照两个方向分别作一次一维傅里叶变换。

$$F(u,v)=\sum\sum f(x,y)\exp\left\{\frac{-2\pi j(x-1)(u-1)}{m}+\frac{-2\pi j(y-1)(v-1)}{n}\right\}$$

$F(1,1)$称为离散傅里叶变换的低频系数，其余系数为高频傅里叶系数。低频系数表达了函数的区域总和，高频系数在不同程度刻画函数的变化情况

fft2(X) 计算矩阵或图像X的二维傅里叶变换

ifft2(Y) 计算离散傅里叶变换逆变换

fftshift(FA) 将零频分量移到频谱中心（平移傅里叶变换函数将大值置于中部）

傅里叶变换对于自然图像往往不具有稀疏性。对于小块亮度函数，其频谱线会集中在坐标轴，并且保持旋转不变性。

14.9 二维离散余弦变换

二维离散余弦变换与二维离散傅里叶变换类似，但不同点为余弦变换是实变换，并且更贴近于真实图像的频谱分析需求（因为真实图像在二维空间往往并不具有周期性）

$$C(u,v)=a(u)b(v)\sum_{x=1}^m\sum_{y=1}^n f(x,y)cos\frac{(2x-1)(u-1)\pi}{2m}cos\frac{(2y-1)(v-1)\pi}{2n}$$

其中$a(u)=\sqrt{2/m},u\neq 1,a(1)=\sqrt{1/m}$以及$b(v)=\sqrt{2/n},v\neq 1,b(1)=\sqrt{1/n}$

$C(1,1)$称为离散余弦变换的低频系数，其余系数为高频余弦系数。

dct2(X) 计算矩阵或图像X的二维离散余弦变换

idct2(Y) 离散余弦逆变换

14.10 图像去噪声问题

与信号去噪问题相似，图像去噪问题即图像对应的目标函数受到了一些未知的退化变换，得到了一幅带有噪点的不准确图像。去噪问题即利用图像的内在性质与噪声分布的特点，完成目标清晰图像的估计与获取的过程。

imnoise(A,'gaussian',mu,sigma) 对图像A添加高斯噪声，均值$\mu$，方差$255^2\sigma$，imnoise可以自动将结果化为uint8型

14.10.1 均值滤波

imfilter(B,filter,'symmetric','same') 需要先定义filter向量，如filter = 1/25 * ones(5)

峰值信噪比(PSNR) 基于图像函数误差定义信噪比，可以对图像的恢复质量进行数值上的分析

$$PSNR=10\cdot\lg(\frac{255^2\cdot\#\text{points}}{\sum [y-f(x)]^2})$$

psnr(pic,realpic) 计算PSNR

14.10.2 中值滤波

imnoise(A,'salt & pepper',0.1) 设置10%的随机像点为椒盐噪声

认为在椒盐噪声噪点周围，大多数信息是准确的，可以在噪点附近一定范围内取灰度值的中位数来估计灰度。

C = medfilt2(B,[3,3]) 3x3的中值滤波

14.10.3 图像区域填充问题

图像整块缺失时，利用逐步填充的思想，利用已知或未擦除区域对剩余区域进行估计。反复进行中值滤波，真实值会渐渐渗透到填充区内部，从而达到图像填充的理想效果

14.10.4 图像模糊的退化模型

MATLAB内置有一系列预设好的filter，可以通过fspecial函数创建模糊核

kernel = fspecial('disk',10) 对焦不准模糊核，半径为10

kernel = fspecial('motion',15,90) 移动模糊核，向上平移15格

imfilter(A,kernel,'replicate') 复制型边界条件

imfilter(A,kernel,'circular') 周期型边界条件

deconvreg(B,kernel) 图像反卷积函数

deconvblind(B,kernel) 盲反卷积函数，此处kernel为初值

Lecture 15 MATLAB数据可视化（二）

15.1 绘图参数设置

subplot(m,n,k) 设定的mxn子图中的第k个位置绘图（行优先）

subplot('position',[left,bottom,width,height]) 指定位置绘图，0~1，左下为0，右上为1

[x,y] = ginput(n) 采用鼠标点击的方法获取二维图形坐标，返回n个，以列向量方式存储

plot3(X,Y,Z,'s') 三维坐标绘图

view([longtitude,latitude]) 0度，东经为正，北纬为正，或通过view([x,y,z])控制观察点位置

[X,Y] = meshigrid(x,y) 生成以x,y为横纵坐标的网格，分别存储在X和Y中（矩阵）

surf(X,Y,Z,C) 三维曲面，若X与Y缺省则以Z的两个维度下标代替，C为颜色函数，默认色彩按曲面高度决定每个点颜色，C指定为大小与Z相同的矩阵，或指定为RGB三元组的m×n×3数组

mesh(X,Y,Z,C) 三维网格线，若X与Y缺省则以Z的两个维度下标代替

colormap (h1,parula) 对子图h1=subplot(1,2,1)选定预设的色图

shading options options可选参数faceted一致颜色并显示边缘，flat一致颜色，interp用双线性插值使四顶点颜色连续变化

alpha(v) 设置透明度，0为完全透明，1完全不透明

light('Name',value) 灯光设置：'color',[1,0,0]色光设置，'style'可选'infinite'平行光或'local'点光源，'position'光源位置

lighting options 照明模式，options可选参数flat均匀分布，gouraud先对顶点颜色插补再对面色插补，phong顶点法线插值，再计算各像素反光，none关闭

material[ka,kd,ks,n,ns] 光反射材质的环境/漫反射/镜面反射强度/镜面反射指数/镜面反射颜色反射率

material options options可选参数shiny使对象比较明亮，镜面反射额较大，dull使对象比较暗淡，漫反射额较大，metal使对象有金属光泽，defalut缺省值

sphere(n) 生成每一维度n等分的三维单位球面

hidden on/off 控制网格叠压图形的透视

pcolor(X,Y,Z) 伪彩图绘制，二维图片用颜色表现第三维

contour(X,Y,Z) 创建一个包含矩阵Z的等值线的等高线图

[C,h] = contourf(X,Y,Z,n,'k:') 带颜色填充的等位线图，n为分级个数

clabel(C,h) 对contour对象添加标识

四维信息：三维曲面+颜色

slice(X,Y,Z,V,sx,sy,sz) 绘制切片图，X,Y,Z,V为同型三维矩阵，表示一组点的各维度坐标和颜色函数值。切片方向与坐标面平行，每一个sx,sy,sz的取值都会带来一个新的切片绘制。

comet(x,y,p)或comet3(x,y,z,p) 彗星线，其长度为p*length(y)或p*length(z)

spinmap(t,inc) 色图变幻，t为变化时间，inc标识变幻的速度

15.2 动画

M(i)=getframe 从当前绘图窗口保存图形到结构体数组，i为帧序号

M(i)=struct('cdata',[],'colormap',[]) 包含cdata和colormap两个域

movie(M,k,fps) 将画面按顺序以每秒fps帧（默认12）的速度播放，重复k次

drawnow 将动态的绘图显示在绘图窗口，用pause(n)可以辅助控制两次绘图之间的时间差

imwrite(I,map,'shuxueshiyan.gif','gif','Loopcount',inf,'DelayTime',0.5) 保存gif，第一帧需要设置'Loopcount',inf，其后只需要设置'WriteMode','append'

附录

倍增区间极值搜索法

fminbnd(fx,a,b)可以计算某个区间内的一个极小值，对于波动较大的函数，极小值不一定是最小值

倍增区间极值搜索法，将分割区间的长度不断变小，从而可以获得“更多的极小值”，每一次分割都将获取最小的一个极小值。在单调区间个数有限的假设下，分割越细，就越有可能获得最小值。

最终当进一步分层后，获取的最小值在连续几次倍增后不变（或变化量极小），就退出循环，结束区间倍增，然后输出最终得到的最小值、最小值点以及最终倍增区间的个数。

function [xmin,fmin,n]=exm060101(fx,a,b,Nt)
%fx为函数,a,b为左右端点，Nt为子区间细分的过程中极小值点重复超过多少次即可停止细分
%xmin和ymin代表最小值所在的横坐标和纵坐标，n为获得最小值时分割的子区间数

%求fx在(a,b)一个极小值，f0记录最小值，第一个返回值（最小值点）舍去
[~,f0]=fminbnd(fx,a,b);
n=1;%初始区间没有分割，共1个区间
jj=1;%同一最小值的重复出现次数
while 1 %无限的循环（需要带有条件的break跳出循环）
n=2*n;%循环区间数X2
d=(b-a)/n;%区间长度/2
    x=a:d:b;
ii=0;%本轮循环（本次精度）极小值的数量
xc=zeros(1,n);fc=xc;%本次循环极值点与极值的点的数组
    for k=1:n%
        [w,f,eflag]=fminbnd(fx,x(k),x(k+1));
        %子区间取极值点，w为点，f为值，flag为是否找到
        if eflag>0
            ii=ii+1;%极小值数量+1
            xc(ii)=w;
            fc(ii)=f;
        end
    end
    xc = xc(1:ii);fc = fc(1:ii);
    [fmin,kk]=min(fc);%最小值fmin，位置kk
    xmin=xc(kk);
    if abs(f0-fmin)<1e-6 %足够小的最小值
        jj=jj+1;
        if jj>Nt
            break
        end
    elseif f0-fmin>1e-6
        f0=fmin;%确立新最小值
        jj=1;%reset重复次数
    end
end

Hilbert矩阵的创建

N = 1000 %阶数
n = repmat(1:N,N,1);
m = n';
A = 1./(n+m-1);