矩陣運算
矩陣的基本生成
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
-x24
方程1:ax^2+bx+c=0
syms a b c x
solve(a*x^2 + b*x + c, x)
(-b+b2-4 a c2 a-b-b2-4 a c2 a)
方程2:\(2\mathrm{x}-\mathrm{x}^2=\mathrm{e}^{-\mathrm{x}}\)
syms x;
y = 2*x - x^2 -exp(-x);
solve(y, x)
0.41640296239270530567959997618172
方程組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
5a b-b+1
res.y
5 a-1a b-b+1
方程組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
0.35324661959671746608371888721268
res.x2
0.60608173664146473530299588999127
方程和方程組的數值解(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 =
ex μ2-μ2+42 μ+μ2+42 μ2+4-ex μ2+μ2+42 μ-μ2+42 μ2+4
方程組1:
y
1
=
y
2
y
2
=
-
y
1
,初值
y
1
0
=
1
y
2
0
=
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 cos(x-π4)
2 cos(x+π4)
方程組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