矩陣運算

矩陣的基本生成

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