矩阵运算
矩阵的基本生成
m1 = 1:5 % 生成行矩阵[1,2,3,4,5]
m2 = 1:2:10 % 起点:步长:终点 [1,3,5,7,9]
linspace(x1,x2,n) % 生成 n 个点。这些点的间距为 (x2-x1)/(n-1)。
m3 = linspace(0,5,11) % [0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5]
m4 = 1:3;
m5 = 4:6;
m6 = [m4, m5] % [1,2,3,4,5,6] 矩阵组合
m7 = [m4; m5] % [1,2,3]
% [4,5,6]
特殊矩阵
eye(n) % 单位矩阵
zeros(n) % n*n的全零矩阵
zeros(m,n) % m*n的全零矩阵
zeros([m,n])
ones(n) & n*n的全一矩阵
ones(m,n) % m*n的全一矩阵
ones([m,n])
随机数矩阵
rand % 生成0~1随机数,符合均匀分布
rand(n) % 生成n*n的随机数矩阵
rand([m,n]) % 生成m*n的随机数矩阵
randi(max) % 1~max的随机整数
randi(max,m) % 1~max的随机m阶整数矩阵
randi(max,[m,n]) % 1~max的随机m*n整数矩阵
randn(n) % 生成n*n的随机数矩阵(服从正态分布)
randn([m,n]) % 生成m*n的随机数矩阵(服从正态分布)
获取矩阵的行列数
[row,col] = size(m) % 返回矩阵m的行数和列数的矩阵
矩阵转置、逆矩阵
m\' % 矩阵m的转置
inv(m) % 矩阵m的逆矩阵
特征值、特征向量
[V, D] = eig(m) % V:特征向量矩阵,D:特征值矩阵,对角元素为特征值
加减乘除、乘方
a = [1,2,3]
b = [4,5,6]
a + b
a - b
a * b % 行列元素对应相乘求和返回矩阵
a .* b % 仅仅是相同位置元素相乘
a / b % 等价于a * inv(b)
a ./ b % 仅仅是相同位置元素相除
a \ b % 等价于inv(a) * b 求 Ax = b的解
a ^ n % a的次方
a .^ n % 每个元素的n次方
广播机制
a = [1,2;4,5]
b = 1
a + b % 把b广播成[1,1;1,1]
逻辑运算
a = [1,2,3]
b = [1,3,2]
a == b % 返回比较矩阵,对应元素进行判断,相同为1,不同为0
a > b
a < b
a == 1 % 也满足广播机制
% 保留a大于2的元素,将小于等于2的元素置为0
a .* (a > 2)
矩阵的索引
m = [1,2,3;4,5,6]
m(3) = 2 % 按列检索
m(row,col) % 返回第row行,第col列的值
m([1,2],[1,2]) % 前两行,前两列的元素
m(1:2, 1:2)
MATLAB基础语法
变量类型和转化
- 数值类型
默认是double类型,可进行加减乘除等运算
- 字符串类型
s1 = "hello";
s2 = \'world\';
[s1,s2] % 字符串类型矩阵 [hello,world]
s1 + s2 % helloworld 用+实现字符串拼接
- 字符串和数值的转化
str2num() % 字符串转数字
str2num("5")+4 % 9
num2str() % 数值转字符串
num2str(num,n) % 把数字num转化成n位小数
num2str(1/3,2) % 0.33
输入和输出
- 输入语句
value = input("请输入一个值:") % 输入数值、矩阵
string = input("请输入一个字符串:",\'s\') % 输入字符串,需要加第二个参数\'s\'
- 输出语句
disp() % 输出多个字段时,需要将多个字段转化成字符串矩阵,再进行输出
disp(["hello",2]) % ["hello","2"]
disp(["1/3=",num2str(1/3,2)]) % "1/3=" "0.33"
运算符和if else控制语句
- 关系运算符
1 < 2 % 1表示真,非0都是真
1 > 2 % 0表示假
- 逻辑运算符
% & 与
% | 或
% ~ 非
(1 < 2) & (1 > 2) % 0
(1 < 2) | (1 > 2) % 1
~((1 < 2) & (1 > 2)) % 1
- if else 语句
value = input("请输入数字:")
if(value == 1)
disp("1既不是质数也不是合数")
elseif(isprime(value))
disp("输入是质数")
else
disp("输入是合数")
end
for循环和while循环
- for循环
% 求1加到100
sum = 0;
for i = [1:100]
sum = sum + i;
end
% 二重for循环
for i = 1:3 % 执行3次
for j = 1:3 % 执行3次
disp([i,j]) % 一共执行9次
end
end
% continue语句(结束当前循环)
for i = 1:3
for j = 1:3
if(i == 2)
continue;
end
disp("i = "+num2str(i)+", j = "+num2str(j));
end
end
% i = 1, j = 1
% i = 1, j = 2
% i = 1, j = 3
% i = 3, j = 1
% i = 3, j = 2
% i = 3, j = 3
% break语句
for i = 1:3
for j = 1:3
if(j == 2)
break;
end
disp("i = "+num2str(i)+", j = "+num2str(j));
end
end
% i = 1, j = 1
% i = 2, j = 1
% i = 3, j = 1
% return语句
for i = 1:3
for j = 1:3
if(j == 2)
return;
end
disp("i = "+num2str(i)+", j = "+num2str(j));
end
end
% i = 1, j = 1
- while循环
n = 1
while(n<5)
disp(n)
n = n + 1;
end
练习题
生成[min,max]之间的随机(整数)数矩阵
min和max为用户输入;行列数为用户输入;
小数或整数由用户指定;最后输出结果(如何使每次输出的结果相同?)
相关知识点:input、if……else、rand、randi、disp
val = input("请输入矩阵元素范围[min,max]:");
size = input("请输入矩阵行列数[row,col]:");
isInt = input("请输入元素类型 0)小数 1)整数:");
rand(\'seed\',0); % 加入随机种子使得每次获得的结果一致
if isInt == 0
res = (val(2) - val(1)) * rand(size) + val(1);
else
res = randi(val,size);
end
disp(res)
模拟圆周率pi
参考思路:蒙特卡罗法;点数为用户输入值
相关知识点:input、if……else、for循环或while循环、rand、disp
allPoints = input("请输入生成的总点数:");
count = 0;
for i = 1:allPoints
if rand^2 + rand^2 <= 1
count = count + 1;
end
end
res = count / allPoints * 4;
disp(res);
斐波那契数列求和
斐波那契数列的第n项和前n项和,n为用户输入值1,1,2,3,5,8,13,21,34,55......
递推公式:F[n] = F[n-1] +F[n-2];
相关知识点: input、for循环或 while循环、disp
N = input("请输入数列的项数:");
res = [1,1];
for i = 3:N
element = res(length(res)) + res(length(res) - 1);
res = [res,element];
end
disp(res); % 输出斐波那契数列
圆中四只鸭子在同一个半圆的概率
参考思路:蒙特卡罗法(进行N次试验,每次试验生成4个随机点,统计四点在同一个半圆的个数)
相关知识点: input、if...else、for循环或while循环、rand、disp
函数、匿名函数和脚本
函数
% function [输入参数] = myfun(输入参数)
% 函数体
% end
x = area(2);
function s = area(r) % 圆的面积
s = pi*r.^2;
end
匿名函数
% f = @(输入参数) 函数体
% f:函数句柄
f = @(x) x.^2;
f(2) % 4
f([2,3]) % [4,9]
f1 = @(x, y) x.^2 + y.^2 + 2*x.*y;
f1(2,3) % 25
f1([1,2,3], 6) % [49, 64, 81]
匿名函数和函数的转化
f2 = @fun; % 方式一
f2(5) % 25
f3 = @(x)fun(x); % 方式二
f3(6)
function y = fun(x)
y = x.^2;
end
使用场景
% 函数体复杂时,使用函数
% 函数体简单时,使用匿名函数
f4 = @fun1;
f4(-5);
f5 = @(x)fun1(x);
f5(-4)
f6 = @fun2;
f6(5,2)
function y = fun1(x)
if x>=0
y = x;
else
y = -x;
end
end
function y = fun2(x,a)
if x>=a
disp("x大于等于a")
else
disp("x小于a")
end
end
脚本
实际上就是后缀.m的文件;
当文件里只有函数时,就成为函数脚本文件或函数文件;
函数文件可被其他脚本调用(需要在同一文件目录下),也可在命令行调用。
value = input("请输入矩阵元素范围[min,max]:");
sz = input("请输入矩阵行列数[row,col]:");
isInt = input("请指定元素类型 0) 小数 1) 整数:");
m = GetMatrix(isInt, value, sz)
disp(m)
GetMatrix.m
function y = GetMatrix(isInt, val, size )
rand(\'seed\',0); % 加入随机种子使得每次获得的结果一致
if isInt == 0
y = (val(2) - val(1)) * rand(size) + val(1);
else
y = randi(val,size);
end
end
MATLAB绘图操作
基本绘图命令plot
- 二维绘图命令plot
绘制 y = sin(x)
x = linspace(0,2*pi,100);
y = sin(x);
plot(x,y);
- 一幅图绘制多条曲线
% 方式一:hold
% hold on 开启图形保持功能,绘制我个图形对象
% hold off 关闭图形保持功能
x = linspace(0,2*pi,100);
y = sin(x);
y2 = cos(x);
plot(x,y);
hold on
plot(x,y2);
hold off
% 方式二:plot(x1, y1, x2, y2, ……)
plot(x,y,x,y2);
- 添加坐标轴标签(lable)、标题(title)、图例(legend)
plot(x,y,x,y2);
xlabel(\'x\');
ylabel(\'y = sin(x), y2 = cos(x)\');
title(\'y = sin(x), y2 = cos(x)\');
legend(\'y = sin(x)\', \'y = cos(x)\');
- 绘制多幅图 figure
x = linspace(0,2*pi,100);
y = sin(x);
y2 = cos(x);
figure; plot(x,y);
figure; plot(x,y2);
- 绘制多个子图 subplot
% subplot(m,n,i):m行n列 第i个图形
x = linspace(0,2*pi,100);
y = sin(x);
y2 = cos(x);
y3 = x.^2;
y4 = x.^0.5;
subplot(2,2,1);plot(x,y);
subplot(2,2,2);plot(x,y2);
subplot(2,2,3);plot(x,y3);
subplot(2,2,4);plot(x,y4);
绘图修饰
% plot(x, y, \'- r o\') 线型、颜色、描点类型(顺序可变)
%线型; -实线、--虚线、:点线、-.点画线
%描点类型: .点、。圆、x叉号、+加号、*星号
% <、>、^、v、(三角形的朝向)
% s方形、d菱形、p五角星、h六角星
%颜色: r红色、g绿色、b蓝色、y黄色、w白色、k黑色
% plot(x, y, \'- r o\') 线性、颜色、描点类型(顺序可变)
x = linspace(0, 2*pi, 30);
y = sin(x);
plot(x,y, \'- r o\');
% plot( x, y, \'- b p\' ,...
% \'linewidth\', 1,... %线宽
% "MarkerEdgecolor\', \'r\' ,... %描点边框颜色
% "MarkerFacecolor \', \'y \' ,... %描点内部填充颜色
% \' Markersize\', 10) %描点大小
x = linspace(0, 2*pi, 30);
y = sin(x);
plot(x, y, \'- b p\', ...
\'lineWidth\',1, ...
\'MarkerEdgeColor\', \'r\', ...
\'MarkerFaceColor\', \'y\', ...
\'MarkerSize\', 10)
% grid on 添加网格
% grid off 取消网格
% axis on 显示坐标轴、刻度线和坐标轴标签
% axis off 关闭坐标轴、刻度线和坐标轴标签
% axis ( [xmin, xmax,ymin,ymax]) 设置x轴和y轴的显示范围
% axis equal 沿每个坐标轴使用相同的数据单位长度。
% axis square 使用相同长度的坐标轴线。相应调整数据单位之间的增量。
x = linspace(0,2*pi,100);
y = sin(x);
y2 = cos(x);
y3 = x.^2;
y4 = x.^0.5;
subplot(2,2,1);plot(x,y); grid on;
subplot(2,2,2);plot(x,y2); axis off;
subplot(2,2,3);plot(x,y3); axis ([-10, 10, 0, 100]);
subplot(2,2,4);plot(x,y4); axis equal
绘制gif动图
for i = [500, 1000, 2000, 5000, 10000]
x1 = linspace(0,1,10000);
y1 = (1 - x1.^2).^0.5;
x2 = rand([1,i]);
y2 = rand([1,i]);
count = 0;
for j = 1:i
if x2(j).^2 + y2(j).^2 <= 1
count = count + 1;
end
end
plot(x1,y1,\'k.\',x2,y2,".");
title(num2str(count)+" / 500 * 4 = "+num2str(count/i*4))
axis square;
frame = getframe(gcf); %捕获坐标区或图窗作为影片帧
I = frame2im(frame); %返回与影片帧关联的图像数据
[I,map] = rgb2ind(I,256); %将RGB图像转换为索引图像I,关联颜色图为map
if i == 500
imwrite(I,map,\'test.gif\',\'gif\',\'Loopcount\',inf,\'DelayTime\',0.2);
else
imwrite(I,map,\'test.gif\',\'gif\',\'WriteMode\',\'append\',\'DelayTime\',0.2);
end
end
更多二维绘图命令
errorbar 含误差图的线条
回归分析,拟合曲线的时候常用
% errorbar:含误差条的线图
% errorbar(x,y,err,["both"|"horizontal"|"vertical"])
% 绘制y对x的图,并在每个数据点处绘制一个垂直误差条,可设置误差条方向
x = linspace(0,100,10);
y = x.^0.5;
err = rand(size(x));
errorbar(x,y,err,"both");
histogram 直方图
% histogram:直方图
% histogram(x,n) 基于x创建直方图,n为区间数量
x = randn([1,10000]);
n = 100;
histogram(x,n);
scatter 散点图
% scatter:散点图
% scatter(x,y) 在向量x和y指定的位置创建一个包含圆形的散点图
% 用法类似plot
x = linspace(0, 2*pi, 30);
y = sin(x);
scatter(x, y, \'o\', \'MarkerEdgeColor\',\'b\',"MarkerFaceColor",\'r\');
bar 柱状图
% bar(y) 创建一个条形图,y中的每一个元素对应一个条形
% 如果y是mxn矩阵,则bar创建每组包含n个条形的m个组
y = [2,3; 6,11; 23,26];
bar(y);
饼图
% pie(X,explode) 使用x中的数据绘制饼图。饼图的每个扇区代表X中的一个元素;
% explode将扇区从饼图偏移一定的位置
X = [1 3 0.5 2.5 2];
explode = [0 1 0 1 0];
pie(X, explode);
三维绘图
plot3 三维曲线
% plot3(x1,y1,z1,LineSpec1,……,xn,yn,zn,LineSpec)
x = linspace(0, 6*pi, 200);
y = sin(x);
z = cos(x);
xlabel(\'x\');
ylabel(\'sin(x)\');
zlabel(\'cos(x)\');
title(\'y = sin(x), z = cos(x)\');
plot3(y, z, x);
scatter3 三维散点图
% scatter3 用法类似scatter
x = linspace(0, 6*pi, 300);
y = sin(x);
z = cos(x);
scatter3(x, y, z, \'o\', \'MarkerEdgeColor\',\'b\',"MarkerFaceColor",\'r\');
mesh、surf 三维曲面
meshgrid 用法
x = 1:3;
y = 1:5;
[X,Y] = meshgrid(x,y)
% 绘制 z = x*e^(-(x^2+y^2))
% mesh 绘制的是网格
[x,y] = meshgrid(-10:1:10, -10:1:10);
z = x.*exp(-x.^2-y.^2);
mesh(x,y,z);
% 绘制 z = x*e^(-(x^2+y^2))
% surf 绘制的是曲面
[x,y] = meshgrid(-10:1:10, -10:1:10);
z = x.*exp(-x.^2-y.^2);
surf(x,y,z);
MATLAB文件读写
将数据写入文件
% writetable(m, filename):将m写入名为filename的文件
% 支持的文件扩展名:.txt .csv .xls .xlsm 或 .xlsx
m = rand(4) + 1;
m = round(m, 2, \'decimals\'); % decimals:小数位数; significant:有效数字位数
t = table(m)
m
____________________________
1.51 1.38 1.94 1.59
1.82 1.81 1.88 1.21
1.79 1.53 1.55 1.3
1.64 1.35 1.62 1.47
writetable(t, \'m.txt\',"Delimiter","\t","WriteVariableNames",false); % 相对路径,设置数据分隔方式为空格
% Delimiter(指定分隔符):默认","、"\t"、";"、"|"
% WriteVariableNames:是否显示列名
writetable(t,\'E:\Project\MATLABProject\process\m.txt\'); % 绝对路径
writetable(t,\'E:\Project\MATLABProject\process\m.csv\');
writetable(t,\'E:\Project\MATLABProject\process\m.xls\');
writetable(t,\'E:\Project\MATLABProject\process\m.xlsm\');
writetable(t,\'E:\Project\MATLABProject\process\m.xlsx\');
type m.txt % 显示内容
1.51 1.38 1.94 1.59
1.82 1.81 1.88 1.21
1.79 1.53 1.55 1.3
1.64 1.35 1.62 1.47
% 将多个矩阵保存在同一个文件 "append"
t2 = table(eye(4))
Var1
________________
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
% writeMode(写入模式):\'overwrite\'(覆盖,默认)、\'append\'追加
writetable(t2, \'m.txt\',"Delimiter","\t","WriteVariableNames",false,"WriteMode","append");
type m.txt;
1.51 1.38 1.94 1.59
1.82 1.81 1.88 1.21
1.79 1.53 1.55 1.3
1.64 1.35 1.62 1.47
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
从文件读取数据
% t = readtable(filename) 从filename 文件中读取数据
% 支持的拓展名:.txt .csv .xls .xlsb .xlsm .xlsx .xltm .xltx
t = readtable("m.txt")
Var1 Var2 Var3 Var4
____ ____ ____ ____
1.51 1.38 1.94 1.59
1.82 1.81 1.88 1.21
1.79 1.53 1.55 1.3
1.64 1.35 1.62 1.47
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
% table转化成数组
m = table2array(t)
列 1 至 3
1.5100 1.3800 1.9400
1.8200 1.8100 1.8800
1.7900 1.5300 1.5500
1.6400 1.3500 1.6200
1.0000 0 0
0 1.0000 0
0 0 1.0000
0 0 0
列 4
1.5900
1.2100
1.3000
1.4700
0
0
0
1.0000
% 读取excel表单
t = readtable(\'E:\Project\MATLABProject\process\m.xls\')
m_1 m_2 m_3 m_4
____ ____ ____ ____
1.51 1.38 1.94 1.59
1.82 1.81 1.88 1.21
1.79 1.53 1.55 1.3
1.64 1.35 1.62 1.47
% 通过名称获取excel表单
t_grade = readtable("student.xls","Sheet","grade") % 指定从student.xls的grade表单读取数据
Var1 ID Chinese Math English
____________ ____ _______ ____ _______
{\'zhangsan\'} 1001 98 94 95
{\'lisi\' } 1002 94 99 98
{\'wangwu\' } 1003 95 95 97
% 通过顺序获取excel表单
t_grade = readtable("student.xls","Sheet",1) % 指定从student.xls的第一个表单读取数据
Var1 ID Chinese Math English
____________ ____ _______ ____ _______
{\'zhangsan\'} 1001 98 94 95
{\'lisi\' } 1002 94 99 98
{\'wangwu\' } 1003 95 95 97
t_info = readtable("student.xls","Sheet",2)
Var1 ID Height Weight
____________ ____ ______ ______
{\'zhangsan\'} 1001 126 54
{\'lisi\' } 1002 128 56
{\'wangwu\' } 1003 135 55
% 获取student.xls所有表单名称
sheets = sheetname("student.xls")
"grade"
"info"
% 获取表单个数
length(sheets)
ans = 2
% 指定单元格范围获取
readtable("student.xls","Range","B2:E4")
Var1 Var2 Var3 Var4
____ ____ ____ ____
1001 98 94 95
1002 94 99 98
1003 95 95 97
table的更多用法
% table的构造
Names = {\'zhangsan\';\'lisi\';\'wangwu\'};
ID = {1001;1002;1003};
Chinese = {98;94;95};
Math = {94;99;95};
English = {95;98;97};
table(Names,ID,Chinese,Math,English)
Names ID Chinese Math English
____________ ________ _______ ______ _______
{\'zhangsan\'} {[1001]} {[98]} {[94]} {[95]}
{\'lisi\' } {[1002]} {[94]} {[99]} {[98]}
{\'wangwu\' } {[1003]} {[95]} {[95]} {[97]}
table(ID,Chinese,Math,English,\'RowNames\',Names)
ID Chinese Math English
________ _______ ______ _______
zhangsan {[1001]} {[98]} {[94]} {[95]}
lisi {[1002]} {[94]} {[99]} {[98]}
wangwu {[1003]} {[95]} {[95]} {[97]}
% 访问表格元素
% 1、通过索引(和矩阵一致)
t_grade(1,2)
t_grade(1,1:5) ;
t_grade(1,:); % 获取第一行所有列
t_grade(:,[1,3]) % 获取所有行和第1、3列;
% 2、通过列名获取
t_grade(:,"ID")
ID
____
1001
1002
1003
t_grade(:,{\'Var1\',\'Chinese\'})
Var1 Chinese
____________ _______
{\'zhangsan\'} 98
{\'lisi\' } 94
{\'wangwu\' } 95
% 修改列名
t_grade.Properties.VariableNames
列 1 至 3
{\'Var1\'} {\'ID\'} {\'Chinese\'}
列 4 至 5
{\'Math\'} {\'English\'}
t_grade.Properties.VariableNames(1) = {\'Name\'}
Name ID Chinese Math English
____________ ____ _______ ____ _______
{\'zhangsan\'} 1001 98 94 95
{\'lisi\' } 1002 94 99 98
{\'wangwu\' } 1003 95 95 97
% 增加行
t_grade(4,:) = {\'zhao\',1004,98,95,100}
Name ID Chinese Math English
____________ ____ _______ ____ _______
{\'zhangsan\'} 1001 98 94 95
{\'lisi\' } 1002 94 99 98
{\'wangwu\' } 1003 95 95 97
{\'zhao\' } 1004 98 95 100
% 增加列
t_grade.total = t_grade.Chinese + t_grade.Math + t_grade.English
Name ID Chinese Math English total
____________ ____ _______ ____ _______ _____
{\'zhangsan\'} 1001 98 94 95 287
{\'lisi\' } 1002 94 99 98 291
{\'wangwu\' } 1003 95 95 97 287
{\'zhao\' } 1004 98 95 100 293
% 合并表格
t_grade = readtable("student.xls","Sheet","grade")
Var1 ID Chinese Math English
____________ ____ _______ ____ _______
{\'zhangsan\'} 1001 98 94 95
{\'lisi\' } 1002 94 99 98
{\'wangwu\' } 1003 95 95 97
t_info = readtable("student.xls","Sheet","info")
Var1 ID Height Weight
____________ ____ ______ ______
{\'zhangsan\'} 1001 126 54
{\'lisi\' } 1002 128 56
{\'wangwu\' } 1003 135 55
join(t_grade,t_info)
Var1 ID Chinese Math English Height Weight
____________ ____ _______ ____ _______ ______ ______
{\'zhangsan\'} 1001 98 94 95 126 54
{\'lisi\' } 1002 94 99 98 128 56
{\'wangwu\' } 1003 95 95 97 135 55
% 创建student表单,并把数据写入
writetable(t_student, \'student.xls\',"Sheet", \'student\')
图像处理
图片的读写和显示
% 像素点的个数对应矩阵的大小
% 矩阵元素值的范围:0:255,0:黑色,255:白色
% uint8:unsigned int 8:无符号整型
% 彩色图片是 m*n*3 的RGB三通道矩阵
pic_bw = imread("JLU.jpg")
部分矩阵:
% 显示图片
imshow(pic)
% 显示部分图片
tmp = pic(60:100,60:100,1:3);
imshow(tmp)
% 构造图像
m = randi([0,255],[400,400]);
m = uint8(m); % double转换成uint8
imshow(m)
% 生成一个渐变图像
bw = zeros([256,400]);
for i = 1:256
for j = 1:400
bw(i,j) = i-1;
end
end
bw = uint8(bw);
imshow(bw)
% 保存图片
imwrite(bw,"bw.jpg");
彩色图、灰色图和二值化
RGB分离与合并
pepper = imread("peppers.png"); % matlab自带图片
imshow(pepper);
% RGB的分离
R = pepper(:,:,1);
G = pepper(:,:,2);
B = pepper(:,:,3);
subplot(2,2,1);
imshow(pepper);
title("original");
subplot(2,2,2);
imshow(R);
title("R");
subplot(2,2,3);
imshow(G);
title("G");
subplot(2,2,4);
imshow(B);
title("B");
% RGB的合并
subplot(1,1,1);
rgb(:,:,1) = R;
rgb(:,:,2) = G;
rgb(:,:,3) = B;
imshow(rgb);
彩色图转灰度图
% rgb2gray 彩色图转灰度图
pepper_gray = rgb2gray(pepper);
imshow(pepper_gray)
二值化
[row, col] = size(pepper_gray);
for i = 1:row
for j = 1:col
if pepper(i,j) > 128
pepper_gray(i,j) = 1;
else
pepper_gray(i,j) = 0;
end
end
end
figure;
pepper_bw = logical(pepper_gray);
imshow(pepper_bw);
% imbinarize 系统自带的二值化方式
% method -用于二值化图像的方法: \'global\'(默认)| \'adaptive\'
% \'Sensitivity\' -自适应阈值的敏感度因子: 0.50(默认)| [0,1]范围内的数值
% \'ForegroundPolarity\' -确定哪些像素被视为前景像素: \'bright\'(默认)| \'dark \'
% \'bright\': 前景比背景亮
% \'dark\' : 前景比背景暗
pepper_gray = rgb2gray(pepper);
bw = imbinarize(pepper_gray, "adaptive","ForegroundPolarity","dark","Sensitivity",0.4);
imshow(bw)
二值化的应用:突出文字
I = imread(\'printedtext.png\');
imshow(I);
bw = imbinarize(I,"adaptive","ForegroundPolarity","dark","Sensitivity",0.4);
imshow(bw)
图像处理相关函数
调整图片大小
% imresize 调整图片大小
% I = imresize(pic, scale); scale为缩放倍数
% I = imresize(pic, [row col]); 调整大小为row*col
I = imread("peppers.png");
imshow(I)
J = imresize(I, 0.5);
imshow(J)
K = imresize(I, [200 200]);
imshow(K)
旋转图像
% I = imrotate(pic,angle); angle为旋转的角度
J = imrotate(I, 45);
imshow(J)
图像的加减乘除
% imadd() 两幅图像相加时,要求大小一致
% imsubtract() 矩阵的减法
% immultiply() 矩阵的点乘
% imdivide() 矩阵的点除
I = imread(\'rice.png\');
imshow(I)
J = imread(\'cameraman.tif\');
imshow(J)
K = imadd(I,J);
imshow(K)
L = imsubtract(I, J);
imshow(L)
M = immultiply(I, 2);
imshow(M)
直方图和直方图均衡化
% imhisteq 直方图均衡化
% imhist 直方图
I = imread(\'tire.tif\');
imshow(I)
imhist(I)
J = histeq(I);
imshow(J);
imhist(J)
标注连通分量
I = imread(\'coins.png\');
imshow(I)
J = imsubtract(I, imopen(I, strel(\'disk\',50))); % 开运算,先腐蚀后膨胀
% bwlabel 先二值化再标注连通分量
bw = imbinarize(J);
imshow(bw)
L = bwlabel(bw, 8);
max(max(L))
ans = 10
方程求解
方程和方程组的解析解(solve)
solve(方程1,方程2,……,变量1,变量2……)
% 多项式合并: (x + 3x - 5x)x/4
syms x
(x+3*x-5*x)*x/4
方程1:ax^2+bx+c=0
syms a b c x
solve(a*x^2 + b*x + c, x)
方程2:\(2\mathrm{x}-\mathrm{x}^2=\mathrm{e}^{-\mathrm{x}}\)
syms x;
y = 2*x - x^2 -exp(-x);
solve(y, x)
方程组1:\( \begin{cases} \mathrm{x}+\mathrm{by}=5\\ \mathrm{ax}-\mathrm{y}=\mathrm{x}\\ \end{cases} \)
syms a b x y
y1 = x + b*y -5;
y2 = a*x -y -x;
res = solve(y1, y2, x, y);
res.x
res.y
方程组2:\( \begin{cases} \mathrm{e}^{-\mathrm{e}^{-\mathrm{x}_1-\mathrm{x}_2}}=\mathrm{x}_2\left( 1+{\mathrm{x}_1}^2 \right)\\ \mathrm{x}_1\cos \left( \mathrm{x}_2 \right) +\mathrm{x}_2\sin \left( \mathrm{x}_1 \right) =\frac{1}{2}\\ \end{cases} \)
syms x1 x2
y1 = exp(-exp(-x1-x2)) - x2*(1+x1^2);
y2 = x1*cos(x2) + x2*sin(x1)-1/2;
res = solve(y1,y2,x1,x2);
res.x1
res.x2
方程和方程组的数值解(fsolve)
fsolve(函数句柄,初值)
% 初值一般根据经验给出
方程:\( 2\mathrm{x}-\mathrm{x}^2=\mathrm{e}^{-\mathrm{x}} \)
f = @(x)2*x-x^2-exp(-x);
fsolve(f,0)
% ans = 0.4164
方程组1:\( \begin{cases} \mathrm{x}+\mathrm{by}=5\\ \mathrm{ax}-\mathrm{y}=\mathrm{x}\\ \end{cases} \)
a = 3;
b = 5;
f = @(x)funs(x,a,b);
fsolve(f,[0,0])
function y = funs(x,a,b) % x = [x,y]
y(1) = x(1) + b*x(2) - 5;
y(2) = a*x(1) - x(2)-x(1);
end
ans = 0.4545 0.9091
方程组2:\( \begin{cases} \mathrm{e}^{-\mathrm{e}^{-\mathrm{x}_1-\mathrm{x}_2}}=\mathrm{x}_2\left( 1+{\mathrm{x}_1}^2 \right)\\ \mathrm{x}_1\cos \left( \mathrm{x}_2 \right) +\mathrm{x}_2\sin \left( \mathrm{x}_1 \right) =\frac{1}{2}\\ \end{cases} \)
f = @fun;
fsolve(f,[0,0])
function y = fun(x)
y(1) = exp(-exp(-x(1)-x(2))) - x(2)*(1+x(1)^2);
y(2) = x(1)*cos(x(2)) + x(2)*sin(x(1)) - 1/2;
end
function y = funs(x,a,b) % x = [x,y]
y(1) = x(1) + b*x(2) - 5;
y(2) = a*x(1) - x(2)-x(1);
end
ans = 0.3532 0.6061
常微分方程和常微分方程组的解析解
dsolve(方程1,方程2,……初值1,初值2……)
方程1:\( \mathrm{y}\'=2\text{x,初值y}\left( 0 \right) =1 \)
syms y(x) % 声明y是x的函数
eqn = diff(y) == 2*x; % diff表示导数
cond = y(0) == 1; % 初值
dsolve(eqn, cond)
% ans = x^2 + 1
方程2:\( \mathrm{y}\'\'=\mathrm{\mu y}\'+\text{y,初值y}\left( 0 \right) =1,\mathrm{y}\'\left( 0 \right) =0 \)
syms y(x) mu;
eqn = diff(y,2) == mu*diff(y)+y;
cond1 = y(0) == 1;
Dy = diff(y);
cond2 = Dy(0) == 0;
dsolve(eqn,cond1, cond2)
ans =
方程组1:
syms y1(x) y2(x)
eqn1 = diff(y1) == y2;
eqn2 = diff(y2) == -y1;
cond1 = y1(0) == 1;
cond2 = y2(0) == 1;
res = dsolve(eqn1,eqn2,cond1,cond2);
res.y1
方程组2:\( \begin{cases} \mathrm{y}_1=\mathrm{a}-\left( \mathrm{b}+1 \right) \mathrm{y}_1+{\mathrm{y}_1}^2\mathrm{y}_2\\ \mathrm{y}_2\'=\mathrm{by}_1-{\mathrm{y}_1}^2\mathrm{y}_2\\ \end{cases}\text{,初值}\begin{cases} \mathrm{y}_1\left( 0 \right) =3\\ \mathrm{y}_2\left( 0 \right) =4\\ \end{cases} \)
syms y1(x) y2(x) a b
eqn1 = diff(y1) == a-(b+1)*y1+y1^2*y2;
eqn2 = diff(y2) == b*y1 - y1^2*y2;
dsolve(eqn1,eqn2)
% 警告: Unable to find symbolic solution.
比较复杂的方程无法求出解析解,此时需要寻求数值解
常微分方程与常微分方程组的数值解
ode45(函数句柄,积分区间,初值)
基本思想:
例如y\' =2x,可以转化为(y[n]一y[n-1])/△x=2x→ y[n]= y[n-1]+2x△x,然后通过迭代的方式来求解y。
也因此,数值解法必须提供初值。
常微分方程(组)求解solver = ode45, ode23, ode113, ode15s, ode23s
- ode45: 4-5阶Runge-Kutta法
- ode23: 2-3阶Runge-Kutta法
- ode113: Adams-Bashforth-Moutlon PECE算法
- ode15s: 后向差分
- ode23s: 修正的二阶Rosenbrock公式
方程1: \( y\'=2x,初值y(0)=10,积分区间[0,10] \)
% 匿名函数必须同时接受两个输入(x,y),计数其中一个输入未使用也是如此
f = @(x,y)2*x;
tspan = [0,10]; % 积分区间
y0 = 10;
[x,y] = ode45(f,tspan,y0);
plot(x,y)
方程2:\(
y\'\'=u(1-y^2)y\'+y,初值y(0)=1,y\'(0)=0,积分区间[0,20]
\)
\(
首先转化为MATLAB标准求解格式:令y_1=y,y_2=y\'\text{,则转化为}
\\
\begin{cases}
\mathrm{y}_1\'=\mathrm{y}_2\\
\mathrm{y}_2\'=\mathrm{\mu}\left( 1-{\mathrm{y}_1}^2 \right) \mathrm{y}_2+\mathrm{y}_1\\
\end{cases}
\)
mu = 1;
f = @(x,y)fun(y,mu);
tspan = [0 20];
y0 = [1 0];
[x,y] = ode45(f,tspan,y0);
plot(x,y(:,1),"r",x,y(:,2),"g")
function ydot = fun(y,mu)
ydot = [y(2);mu*(1-y(1)^2*y(2)+y(1))];
end