Matlab 從入門到放棄: Application of MATLAB in Engineering 台大郭彥甫

本文轉載自查看原文 2021-01-02 23:03 347 Matlab

Lecture 01 學習導引：

Lecture 02 基本操作與矩陣操作：

Introduction

基本界面：

有兩種編程方式：

MATLAB as calculator

^:這里是次方號！

Note:

Command Window 如何清屏： clc 命令

遇到不熟悉的命令，可以使用谷歌，可以使用Matlab 自帶文檔，也可以在命令行中使用help ,

例如：在Command Window中：

e 如何表示: exp(1) ， e 的平方是 exp(2)

總結起來就是：遇到不會就online help ,

變量：

變量類型：

Matlab 中的關鍵字：

變量命名不規范容易導致下面的問題：

Note : clear + 變量名會將工作區的該變量給清除，如果只有clear 會清空工作區的所有變量。

format 表示

>> format rat
>> 3/13 + 4/14 + 5/15

ans =

232/273

>> format long
>> 3/13+4/14 + 5/15

ans =

0.849816849816850

分號的作用：

一些常用的指令：

Array operation

Array 可以分為矩陣（Matrix）和向量（Vector）。

如何輸入A矩陣呢？

數組索引：

大概是兩種方法：

1，(row ，column) 用逗號分開

2，(idx) 按列索引

小練習：

答案： A(1,2) = 76 ; A(3,2) = 76; A([1 2],[2 3]) = 0;

冒號操作：

練習：如何去掉下面矩陣的第三行？

1, A= A([1 2],:)

2, A(3,:) = []

數組拼接：

三種：

1，空格或逗號橫着拼接

2，分號豎着拼接

A =

1 2
3 4

B =

5 6
7 8

>> C = [A B] // 空格

C =

1 2 5 6
3 4 7 8

>> C = [A,B] // 逗號

C =

1 2 5 6
3 4 7 8

>> C = [A;B] // 分號相當於換行

C =

1 2
3 4
5 6
7 8

答案：

A =

1 2
3 4

B =

9 9
9 9

C =

5 6 7 8

D =

-2 -1 0 1

>> F = [A B;C;D]

F =

1 2 9 9
3 4 9 9
5 6 7 8
-2 -1 0 1

數組操作：

主要有

+ 加

- 減

* 乘

/ 除，要用逆矩陣

^ 次方

. 點乘或點左除或點右除或點次方 dot

' 轉置

一些特殊的矩陣：

注：

diag(A) 得到矩陣A的對角向量。

rand() 返回一個（0,1）的小數

與矩陣相關的常用函數：

注：

find(A>2) ,find 的參數一般跟一個條件表達式，

Lecture 03 結構化程序與自定義函數：

Script writing

Structured programming

注： matlab 中的！= 是 ~=

if 語句：

switch 語句：

while 循環語句：

for 循環

預先申請空間：

預申請空間比動態內存更有效率。

break 和 continue

有用的小技巧

Functions

Note: which Locate functions and files. edit Edit or create a file

User-defined function

函數名一般和文件名一致，一般一個函數自己一個腳本，然后其他腳本里面可以調用它，

有多個返回值：

函數內置變量：

等等，其他自行探索...

匿名函數：

總結：隨用隨定義，只能在定義匿名函數的文件中使用，

函數指針：

%% 函數 做為一個參數  
f1 = @(x) sin(x);
f2 = @(f,arg) f(arg); % 參數 f 為一個函數指針  

ret = f2(f1,pi/6);
disp(ret);

ret2 = f2(@cos,pi/6); % @cos 引用內置的函數 
disp(ret2);

input ，num2str ,isempty 的簡單用法：

Lecture 04 數據結構和文件訪問：

Variables:string ,structure,cell

字符串

字符串的拼接

1 space, 2 ；如下：

反轉字符串：

結構體：

結構體數組：

結構體相關的函數

嵌套的結構體：

Cell :

另一種形式的結構體，

訪問cell 的內容

使用（）訪問的是 pointer ,

使用 {} 訪問的才是具體的內容，

cell 相關的函數

num2cell Convert numeric array into cell array.

mat2cell Break matrix up into a cell array of matrices.

關於mat2cell :

多維數組

使用cat() 函數可以拼接兩個數組

其實，以上的操作也可以用cat 來完成，如下：

所以可以用 cat 來構建多維數組，

reshaper() 函數

判斷數據的類型

Data access

load() 和 save()

將 workspace 工作區的數據保存到文件中：

如何只存儲指定的變量呢？

可以參數指定要存儲的變量名，這里是存儲了a 變量

將 文件中的數據加載到workspace

xlsread() 和 xlswrite()

xlsread()

它會自動去掉Head的部分，只保留數字，

那么如何得到文字部分呢？

得到的head 是個cell

xlswrite()

底層文件操作：

向文件寫入：

Note： type 命令，List program file 。相當於Linux的cat 命令，把文件的內容顯示出來。

從文件中讀出：

這里使用的是 fscanf ,

最終得到的 res 為：

Lecture 05 基礎繪圖：

Basic plotting

plot() 函數：

x = 1:10;
y = sin(x);
plot(x,y);  
plot(x);  % a line. plot(x) will replace plot(x,y)!


%% 
x2 = 1:10;
y2 = sin(x2);
plot(x2,y2);
hold on % change mode, add plot ,not replace. default value is replace!
plot(x2);  % a line. plot(x) will be added!
plot(x2+y2);

注：這里展示的僅僅只是一部分，更多請查看文檔。

x2 = 1:10;
y2 = sin(x2);
plot(x2,y2,'.-r');  % marker is .  line is -  color is red

legend() 函數：

默認的legend 是加在右上角的，當然也可以改變它的位置，

legend(...,'Location',LOC) adds a legend in the specified
    location, LOC, with respect to the axes.  LOC may be either a
    1x4 position vector or one of the following strings:
        'North'              inside plot box near top
        'South'              inside bottom
        'East'               inside right
        'West'               inside left
        'NorthEast'          inside top right (default for 2-D plots)
        'NorthWest'           inside top left
        'SouthEast'          inside bottom right
        'SouthWest'          inside bottom left
        'NorthOutside'       outside plot box near top
        'SouthOutside'       outside bottom
        'EastOutside'        outside right
        'WestOutside'        outside left
        'NorthEastOutside'   outside top right (default for 3-D plots)
        'NorthWestOutside'   outside top left
        'SouthEastOutside'   outside bottom right
        'SouthWestOutside'   outside bottom left
 'Best'           least conflict with data in plot   這個不錯哦
        'BestOutside'        least unused space outside plot

title() 函數和 ?label：

? 可以是 x,y ,z

x3 = 1:0.3:10;
y3 = sin(x3);
hold on % add ,not replace 
plot(x3,y3,'--r');
plot(x3,x3,'--g');

legend('y=sin(x)','y=x') % add legend 

title('title example.Latex:e^{x}. t = 0 to 2\pi')% title , latex 
xlabel('x label test,t = 0 to 10.test:e^x 2\pi')
ylabel('y label test')

關於如何獲取latex 文本： https://latex.91maths.com/

text() and annotation() 函數

x3 = 1:0.3:10;
y3 = sin(x3);
hold on % add ,not replace 
plot(x3,y3,'--r');
plot(x3,x3,'--g');
line([1 10],[0 0]) % line from (1,0) to (10,0)
line([4 4], [0 sin(4)])  % line from (4,0) to (4 ,sin(4))  
line([5 5], [0 sin(5)])  % line from (5,0) to (5 ,sin(5))

legend('y=sin(x)','y=x') % add legend 

title('Just a test.')% title , latex 
xlabel('x label test')
ylabel('y label test')

t = '$$ \int_{4} ^{5} \sin(x) dx $$';   % latex format 
text(1.5,5,t,'Interpreter','latex') % pos is (1.5,5)

annotation('arrow',[0.25 0.45],[0.55 0.2]) % creates an annotation object  
                   %head pos: (0.25,0.55) .* (width,height)% 這里的坐標是 歸一化坐標！  
                   %tail pos: (0.45,0.2) .* (width,height)

當然，annotation 也可以畫其他的形狀，具體看幫助文檔，

'rectangle'
'ellipse'
'textbox'
'line'
'arrow'
'doublearrow' = two headed arrow
'textarrow' = arrow with text at tail end

Graphical object properties

Figure Adjustment

在知道如何調成第二個圖之前，我們首先要先知道一些知識：

1，Graphical Objects

而且它們還有一定的層次關系，

2 Figure Properties

如何進入這里：在圖片界面，View -> Property Editor

改變某個Object的屬性：

第一步，在創建的時候保存，例如 h = plot(x,y);

gca 和 gcf :

x4 = 1:0.3:10;
y4 = sin(x4);
h = plot(x4,y4);  
ret1 = gca;  % current axes  object   
ret2 = gcf;  % current figure object  
disp(ret1);
disp(ret2);

　　output:

allchild() 函數：

x4 = 1:0.3:10;
y4 = sin(x4);
h = plot(x4,y4);  
ret1 = allchild(gca);  
ret2 = allchild(gcf);  
disp(ret1); 
disp(ret2);

OUTPUT：
  Line with properties:

              Color: [0 0.4470 0.7410]
          LineStyle: '-'
          LineWidth: 0.5000
             Marker: 'none'
         MarkerSize: 6
    MarkerFaceColor: 'none'
              XData: [1x31 double]
              YData: [1x31 double]
              ZData: [1x0 double]

  Show all properties
=========================
  10x1 graphics array:

  Menu       (figMenuHelp)
  Menu       (figMenuWindow)
  Menu       (figMenuDesktop)
  Menu       (figMenuTools)
  Menu       (figMenuInsert)
  Menu       (figMenuView)
  Menu       (figMenuEdit)
  Menu       (figMenuFile)
  Toolbar    (FigureToolBar)
  Axes

第二步，使用 get() 和 set() 進行操作：

使用 get

x4 = 1:0.3:10;
y4 = sin(x4);
h = plot(x4,y4);  
disp('----line object-----')
get(h) % line object's properties 
disp('----axes object-----')
get(gca) % axes object's properties 
disp('----figure object-----')
get(gcf) % figure object's properties

----line object-----
    AlignVertexCenters: 'off'
            Annotation: [1x1 matlab.graphics.eventdata.Annotation]
          BeingDeleted: 'off'
            BusyAction: 'queue'
         ButtonDownFcn: ''
              Children: []
              Clipping: 'on'
                 Color: [0 0.4470 0.7410]
             CreateFcn: ''
             DeleteFcn: ''
           DisplayName: ''
      HandleVisibility: 'on'
               HitTest: 'on'
         Interruptible: 'on'
             LineStyle: '-'
             LineWidth: 0.5000
                Marker: 'none'
       MarkerEdgeColor: 'auto'
       MarkerFaceColor: 'none'
            MarkerSize: 6
                Parent: [1x1 Axes]
         PickableParts: 'visible'
              Selected: 'off'
    SelectionHighlight: 'on'
                   Tag: ''
                  Type: 'line'
         UIContextMenu: []
              UserData: []
               Visible: 'on'
                 XData: [1x31 double]
             XDataMode: 'manual'
           XDataSource: ''
                 YData: [1x31 double]
           YDataSource: ''
                 ZData: [1x0 double]
           ZDataSource: ''

----axes object-----
                       ALim: [0 1]
                   ALimMode: 'auto'
     ActivePositionProperty: 'outerposition'
          AmbientLightColor: [1 1 1]
               BeingDeleted: 'off'
                        Box: 'on'
                   BoxStyle: 'back'
                 BusyAction: 'queue'
              ButtonDownFcn: ''
                       CLim: [0 1]
                   CLimMode: 'auto'
             CameraPosition: [5.5000 0 17.3205]
         CameraPositionMode: 'auto'
               CameraTarget: [5.5000 0 0]
           CameraTargetMode: 'auto'
             CameraUpVector: [0 1 0]
         CameraUpVectorMode: 'auto'
            CameraViewAngle: 6.6086
        CameraViewAngleMode: 'auto'
                   Children: [1x1 Line]
                   Clipping: 'on'
              ClippingStyle: '3dbox'
                      Color: [1 1 1]
                 ColorOrder: [7x3 double]
            ColorOrderIndex: 2
                  CreateFcn: ''
               CurrentPoint: [2x3 double]
            DataAspectRatio: [4.5000 1 1]
        DataAspectRatioMode: 'auto'
                  DeleteFcn: ''
                  FontAngle: 'normal'
                   FontName: 'Helvetica'
                   FontSize: 10
              FontSmoothing: 'on'
                  FontUnits: 'points'
                 FontWeight: 'normal'
                  GridAlpha: 0.1500
              GridAlphaMode: 'auto'
                  GridColor: [0.1500 0.1500 0.1500]
              GridColorMode: 'auto'
              GridLineStyle: '-'
           HandleVisibility: 'on'
                    HitTest: 'on'
              Interruptible: 'on'
    LabelFontSizeMultiplier: 1.1000
                      Layer: 'bottom'
             LineStyleOrder: '-'
        LineStyleOrderIndex: 1
                  LineWidth: 0.5000
             MinorGridAlpha: 0.2500
         MinorGridAlphaMode: 'auto'
             MinorGridColor: [0.1000 0.1000 0.1000]
         MinorGridColorMode: 'auto'
         MinorGridLineStyle: ':'
                   NextPlot: 'replace'
              OuterPosition: [0 0 1 1]
                     Parent: [1x1 Figure]
              PickableParts: 'visible'
         PlotBoxAspectRatio: [1 0.7903 0.7903]
     PlotBoxAspectRatioMode: 'auto'
                   Position: [0.1300 0.1100 0.7750 0.8150]
                 Projection: 'orthographic'
                   Selected: 'off'
         SelectionHighlight: 'on'
                 SortMethod: 'childorder'
                        Tag: ''
                    TickDir: 'in'
                TickDirMode: 'auto'
       TickLabelInterpreter: 'tex'
                 TickLength: [0.0100 0.0250]
                 TightInset: [0.0506 0.0532 0.0134 0.0202]
                      Title: [1x1 Text]
    TitleFontSizeMultiplier: 1.1000
            TitleFontWeight: 'bold'
                       Type: 'axes'
              UIContextMenu: []
                      Units: 'normalized'
                   UserData: []
                       View: [0 90]
                    Visible: 'on'
              XAxisLocation: 'bottom'
                     XColor: [0.1500 0.1500 0.1500]
                 XColorMode: 'auto'
                       XDir: 'normal'
                      XGrid: 'off'
                     XLabel: [1x1 Text]
                       XLim: [1 10]
                   XLimMode: 'auto'
                 XMinorGrid: 'off'
                 XMinorTick: 'off'
                     XScale: 'linear'
                      XTick: [1 2 3 4 5 6 7 8 9 10]
                 XTickLabel: {10x1 cell}
             XTickLabelMode: 'auto'
         XTickLabelRotation: 0
                  XTickMode: 'auto'
              YAxisLocation: 'left'
                     YColor: [0.1500 0.1500 0.1500]
                 YColorMode: 'auto'
                       YDir: 'normal'
                      YGrid: 'off'
                     YLabel: [1x1 Text]
                       YLim: [-1 1]
                   YLimMode: 'auto'
                 YMinorGrid: 'off'
                 YMinorTick: 'off'
                     YScale: 'linear'
                      YTick: [1x11 double]
                 YTickLabel: {11x1 cell}
             YTickLabelMode: 'auto'
         YTickLabelRotation: 0
                  YTickMode: 'auto'
                     ZColor: [0.1500 0.1500 0.1500]
                 ZColorMode: 'auto'
                       ZDir: 'normal'
                      ZGrid: 'off'
                     ZLabel: [1x1 Text]
                       ZLim: [-1 1]
                   ZLimMode: 'auto'
                 ZMinorGrid: 'off'
                 ZMinorTick: 'off'
                     ZScale: 'linear'
                      ZTick: [-1 0 1]
                 ZTickLabel: ''
             ZTickLabelMode: 'auto'
         ZTickLabelRotation: 0
                  ZTickMode: 'auto'

----figure object-----
                 Alphamap: [1x64 double]
             BeingDeleted: 'off'
               BusyAction: 'queue'
            ButtonDownFcn: ''
                 Children: [1x1 Axes]
                 Clipping: 'on'
          CloseRequestFcn: 'closereq'
                    Color: [0.9400 0.9400 0.9400]
                 Colormap: [64x3 double]
                CreateFcn: ''
              CurrentAxes: [1x1 Axes]
         CurrentCharacter: ''
            CurrentObject: []
             CurrentPoint: [-192 -204]
                DeleteFcn: ''
             DockControls: 'on'
                 FileName: ''
        GraphicsSmoothing: 'on'
         HandleVisibility: 'on'
            IntegerHandle: 'on'
            Interruptible: 'on'
           InvertHardcopy: 'on'
              KeyPressFcn: ''
            KeyReleaseFcn: ''
                  MenuBar: 'figure'
                     Name: ''
                 NextPlot: 'add'
                   Number: 1
              NumberTitle: 'on'
         PaperOrientation: 'portrait'
            PaperPosition: [0.2500 2.5000 8 6]
        PaperPositionMode: 'manual'
                PaperSize: [8.5000 11]
                PaperType: 'usletter'
               PaperUnits: 'inches'
                   Parent: [1x1 Root]
                  Pointer: 'arrow'
        PointerShapeCData: [16x16 double]
      PointerShapeHotSpot: [1 1]
                 Position: [793 383 560 420]
                 Renderer: 'opengl'
             RendererMode: 'auto'
                   Resize: 'on'
            SelectionType: 'normal'
           SizeChangedFcn: ''
                      Tag: ''
                  ToolBar: 'auto'
                     Type: 'figure'
            UIContextMenu: []
                    Units: 'pixels'
                 UserData: []
                  Visible: 'on'
      WindowButtonDownFcn: ''
    WindowButtonMotionFcn: ''
        WindowButtonUpFcn: ''
        WindowKeyPressFcn: ''
      WindowKeyReleaseFcn: ''
     WindowScrollWheelFcn: ''
              WindowStyle: 'normal'

OUTPUT

使用 set 修改，

x4 = 1:0.3:10;
y4 = sin(x4);
h = plot(x4,y4);   圖像如下：

現在要修改它的屬性，將 x 限制在 [pi,3 pi] y 限制在 [-1.5 ,1.5 ]。

x4 = 1:0.3:10;
y4 = sin(x4);
h = plot(x4,y4); % line object 
set(gca,'xlim',[pi,3*pi]) % x in [pi,3pi] set(gca,'ylim',[-1.5,1.5]) % y in [-1.5,1.5] Just a test!

Setting Font and Tick of Axes

x4 = 1:0.3:10;
y4 = sin(x4);
h = plot(x4,y4); % line object 
set(gca,'xlim',[pi,3*pi]) % x in [pi,3pi]
set(gca,'ylim',[-1.5,1.5]) % y in [-1.5,1.5] Just a test! 
set(gca,'fontsize',15) % set tick' font size equal to 15px

x4 = 1:0.3:10;
y4 = sin(x4);
h = plot(x4,y4); % line object 
set(gca,'xlim',[pi,3*pi]) % x in [pi,3pi]
set(gca,'ylim',[-1.5,1.5]) % y in [-1.5,1.5] Just a test! 
set(gca,'fontsize',15) % set tick' font size equal to 15px
set(gca,'xtick',pi:pi/2:3*pi) % use pi on tick. set(gca,'xticklabel',180:90:3*180)

Setting Line Style

x5 = 1:0.3:10;
y5 = sin(x5);
h = plot(x5,y5); % line object 
set(h,'linestyle','-.',... 'linewidth',7.0,... 'color','g')

delete(h); 將會刪除圖像， delete： Delete file or graphics object.

修改Maker 的 Style:

x5 = 1:10;
y5 = sin(x5);
h = plot(x5,y5,'-d');

x5 = 1:10;
y5 = sin(x5);
h = plot(x5,y5,'-d'); % 設置Marker 的前提是，figure 首先要有 Marker，這里指定Marker 為d,即diamond.
set(h,'markeredgecolor','r',... 'markerfacecolor','g',... 'markersize',10) % set marker style  效果如下：

小練習：

完成圖1 到圖2 的轉換：

%% figure01 
%x = 1:0.02*pi:2;% can't get 2!
x = linspace(1,2,90); % can get 2! num is 90,
y1 = x.^2;
y2 = sin(2*pi*x);
hold on; % add not replace!  
l1 = plot(x,y1,'k');% color is black
l2 = plot(x,y2,'ro');% marker shape is circle,color is red
legend([l1,l2],'t^2','sin(2\pi t)','location','northwest')% legend 
title('Mini Assignment #1')
xlabel('Time (ms)')
ylabel('f(t)')
set(gca,'xtick',linspace(1,2,6));%[1,2] num is 6

%% figure02
%x = 1:0.02*pi:2;% can't get 2!
x = linspace(1,2,90); % can get 2! num is 90,
y1 = x.^2;
y2 = sin(2*pi*x);
hold on; % add not replace!  
l1 = plot(x,y1,'k');% color is black
l2 = plot(x,y2,'ro');% marker shape is circle,color is red
legend([l1,l2],'t^2','sin(2\pi t)','location','northwest')% legend 
title('Mini Assignment #1')
xlabel('Time (ms)')
ylabel('f(t)')
set(gca,'xtick',linspace(1,2,6));%[1,2] num is 6  
% to figure2 set(gca,'fontsize',16) % setting axes's fontsize set(l1,'linewidth',4) % setting l1's linewidth to 4 set(l2,'color',[160,32,240]/255,... % change color to purple 'marker','.',... % change marker to point 'markersize',16)

Multiple Figures

Note : 需要注意的是，此時會有兩個figure，如果使用gcf gca 僅僅指的是最后一個figure，前面的figure就不能訪問了，

%% multipe figures

x = 1:0.1:10;
y1 = x.^2;
y2 = exp(x);
%====
figure('position',[100,200,400,400]);% figure 1
plot(x,y1);
%====
figure('position',[600,600,400,400]);% figure 2
plot(x,y2);

Small figures in One Figure

%% small figures in one figure 
x = linspace(-1,1,50);
y = sqrt(1-x.^2);
x = [x,fliplr(x)];%fliplr Flip array in left/right direction.%reverse 
y = [y,-fliplr(y)];
subplot(221)
plot(x,y);
axis normal 
title('normal[default]')

subplot(222)
plot(x,y);
axis square  %makes the current axis box square in size. x axis's length == y axis's length 
title('square')

subplot(223)
plot(x,y);
axis equal  % x axis's unit == y axis's unit 
title('equal')

%% small figures in one figure 
x = linspace(-1,1,50);
y = sqrt(1-x.^2);
x = [x,fliplr(x)];%fliplr Flip array in left/right direction.%reverse 
y = [y,-fliplr(y)];
subplot(221)
plot(x,y);
axis normal 
title('normal[default]')
grid on box off 

subplot(222)
plot(x,y);
axis square  %makes the current axis box square in size. x axis's length == y axis's length 
title('square')
axis off     

subplot(223)
plot(x,y);
axis equal  % x axis's unit == y axis's unit 
title('equal')
grid on

我的使用習慣：

grid on
box off

保存圖片到文件

Vector 格式是矢量圖放大不失真！

%% save to file
x = linspace(1,10,100);
y = sin(x);
plot(x,y);
saveas(gcf,'test01','png') % bitmap format
saveas(gcf,'test01','svg') % vector format

更強大的保存圖片參考 print ！

lecture 06 高級繪圖

先上一張圖：

Advanced 2D plots

x = logspace(-1,1,100); % yi 10 wei di 
y = x.^2;
subplot(221);
plot(x,y);
title('plot');

subplot(222);
semilogx(x,y);
title('semilogx');

subplot(223);
semilogy(x,y);
title('semilogy');

subplot(224);
loglog(x,y);
title('loglog');

plotyy 函數

Graph with y-tick labels on the left and right side

%% plotyy()  left y label and right y label
x = 1:0.1:10;
y1 = sin(x);
y2 = 2*sin(2*x) + 1; 
[ax,h1,h2] = plotyy(x,y1,x,y2);% the handles of the two axes
                    % AX(1) is the left axes and AX(2) is the right axes.
set(get(ax(1),'ylabel'),'string','Left Y-axis')
set(get(ax(2),'ylabel'),'string','Right Y-axis')
title('Labeling plotyy');
set(h1,'linestyle','--'); 
set(h2,'linestyle',':');  % 如下圖

Histogram 直方圖

看分布是如何的，與概率相關，它是看整體的情況

%% histogram 
x  = randn(1,1000);  % Normally distributed ,x is 1x1000 martix 
subplot(211);
hist(x,10);  % bins is 10. 10 groups  
title('bins is 10')

subplot(212);
hist(x,50);  % bins is 50. 50 groups  
title('bins is 50')  % 效果如下圖

Bar Charts 柱狀圖條形圖

%% bar charts
x = 1:3; % 3 months
y = randi([2,10],1,3); % random 2-10. shape is 1x3 matrix.
subplot(221)
bar(x,y);
set(gca,'xticklabel',{'jan','feb','mar'}); % change xtick to text
title('bar single')
disp(x);
disp(y);
disp('===');

subplot(222)
y2 = randi([2,10],1,3);
y3 = randi([2,10],1,3);
bar(x,[y;y2;y3]); % The length of X must match the number of rows of Y.
title('bar multiple')
disp(x);
disp([y;y2;y3]);
disp('===');

subplot(223)
y = 1:3; % y axis
z = randi([2,10],1,3); % z axis
bar3(y,z); % y and z
title('bar3 single')
disp(y);
disp(z);
disp('===');

subplot(224)
y = 1:3; % y axis
z = randi([2,10],1,3); % z axis
z2 = randi([2,10],1,3); % z axis
z3 = randi([2,10],1,3); % z axis
bar3(y,[z;z2;z3]); % The length of Y must match the number of rows of Z.
title('bar3 multiple')
disp(y);
disp([z;z2;z3]);
disp('===');

OUTPUT
===
     1     2     3

     2     2     4

===
     1     2     3

     2     2     4
     2     5     5
     6    10     4

===
     1     2     3

     5     9     5

===
     1     2     3

     3     6     9
     5     8     7
     6     6    10

===

Stacked and Horizontal Bar Charts

%% stacked and horizontal bar charts 

x = 1:3; % 3 months 
y = randi([2,10],1,3); % random 2-10. shape is 1x3 matrix. 
y2 = randi([2,10],1,3); 
y3 = randi([2,10],1,3); 
subplot(131)
bar(x,[y;y2;y3],'stacked'); 
title('stacked bar')

subplot(132);
barh(x,[y;y2;y3]);
title('horizontal bar');

subplot(133);
barh(x,[y;y2;y3],'stacked');
title('stacked horizontal bar');

disp(x);
disp([y;y2;y3]);

OUTPUT
     1     2     3

     6     6     7
     8     8     4
     8     7     3

Pie Charts

%% Pie Charts 
x = [10 5 20 30];
subplot(221);
pie(x);


subplot(222);
pie(x,[0,0,0,1]);% 將最后一份分割出來 

subplot(223);
pie(x,[1,1,1,1]);% 每一份都分割開來 

subplot(224);
pie3(x,[0,0,0,1]);% 3d 繪圖

Polar Chart 極坐標圖

%% Polar Chart 
theta = 0:0.1*pi:pi/2;
r = 2*cos(theta);
subplot(121)
polar(theta,r);

theta2 = 0:0.1*pi:pi;
r2 = 2*cos(theta2);
subplot(122)
polar(theta2,r2);

更多圖形：

Stairs and Stem Charts 階梯圖和桿圖

%% Stairs and Stem Charts
x = linspace(0,4*pi,40);
y = sin(x);
subplot(121);
stairs(x,y);  % 階梯圖 
title('stairs');

subplot(122);
stem(x,y);  % 桿圖 
title('stem');

Boxplot and Error Bar 箱圖和誤差圖

箱型圖略.

%% Error Bar 
x = 0:pi/10:pi;
y = sin(x);
err = std(y)*ones(size(x));
errorbar(x,y,err);

fill() 函數

%% fill
theta = [1:2:15]*pi/8;
x = cos(theta);
y = sin(theta);
subplot(121)
plot(x,y);axis square;


subplot(122)
fill(x,y,'r');axis square;

Color space

set(lineObj,'color',[160,32,240]/255);

imagesc 函數， Scale data and display as image

%% imagesc 
d = [0 2 4 3; 8 10 12 14; 16 18 20 22];
imagesc(d)
colorbar  % 顯示 color bar

d = [0 2 4 3; 8 10 12 14; 16 18 20 22];
imagesc(d)
colormap(hot); % hot cool gray  
colorbar

內置colormaps

%% built-in colormaps 
x = randi([0,255],10,1);
y = randi([0,255],1,10);
d = mod(x*y,255);
% imagesc(d);  % default colormap is PARULA

%colormap('jet');
%imagesc(d);

colormap('hsv');
imagesc(d);

這個要自定義 ColorMap ,

%% custom colormap 
x = [1:10; 3:12; 5:14];
myColorMap = zeros(256,3);
myColorMap(:,2) = (0:255)/255;
colormap(myColorMap);

imagesc(x);
colorbar;

3D plots

2D vs. 3D

%% 2d vs 3d

x=0:0.1:2*pi;

subplot(121);
plot(x,sin(x));
grid on
xlabel('x-axis'); ylabel('y-axis');

subplot(122);
plot3(x,sin(x),zeros(size(x)));
grid on
xlabel('x-axis'); ylabel('y-axis'); zlabel('z-axis');

Note :3維空間坐標系之---左手坐標系和右手坐標系的區別：

右手坐標系的z軸指向內，左手指向外。

Matlab 使用的是右手坐標系。

3D surface 圖：

%% surface plots
x = linspace(-3,3,20);
y = linspace(-4,4,20);
[X,Y] = meshgrid(x,y);
Z = X.*exp(-X.^2-Y.^2);
subplot(121);
mesh(X,Y,Z);%3-D mesh surface. 3d 網格表面

subplot(122);
surf(X,Y,Z);%3-D colored surface.3d 顏色表面

contour() 函數

contour 是輪廓的意思，

Projection of equal heights of 3D plot onto a 2D plane

%% contour 
x = linspace(-3,3,20);
y = linspace(-4,4,20);
[X,Y] = meshgrid(x,y);
Z = X.*exp(-X.^2-Y.^2);
subplot(121);
surf(X,Y,Z);%3-D colored surface.
axis square;

subplot(122);
contour(X,Y,Z);
axis square;

More Contour Plots

%% more contour 
x = linspace(-3,3,20);
y = linspace(-4,4,20);
[X,Y] = meshgrid(x,y);
Z = X.*exp(-X.^2-Y.^2);
subplot(121);
surf(X,Y,Z);%3-D colored surface.
axis square;

subplot(122);
[C,h] = contour(X,Y,Z);
clabel(C,h); % 標上數值 
axis square

%% more contour

x = linspace(-3,3,20);

y = linspace(-4,4,20);

[X,Y] = meshgrid(x,y);

Z = X.*exp(-X.^2-Y.^2);

subplot(121);

surf(X,Y,Z);%3-D colored surface.

axis square;

subplot(122);

[C,h] = contourf(X,Y,Z);% 塗上顏色 f: fill

clabel(C,h);

axis square;

meshc() and surfc() 函數

Combination of surface/mesh and contours。

%% meshc() and surfc() 
x = linspace(-3,3,20);
y = linspace(-4,4,20);
[X,Y] = meshgrid(x,y);
Z = X.*exp(-X.^2-Y.^2);
subplot(121);
meshc(X,Y,Z);
axis square;

subplot(122);
surfc(X,Y,Z);
axis square;

View Angle: view()

觀察的角度，

Azimuth and Elevation:

Note: AZ = -37.5, EL = 30 is the default 3-D view.
AZ = 0, EL = 90 is directly overhead and the default 2-D view.

%% view angle :view()
% Note:     AZ = -37.5, EL = 30 is the default 3-D view.
%           AZ = 0, EL = 90 is directly overhead and the default 2-D view.
x = linspace(-3,3,20);
y = linspace(-4,4,20);
[X,Y] = meshgrid(x,y);
Z = X.*exp(-X.^2-Y.^2);
subplot(121);
surf(X,Y,Z);
axis square;


subplot(122);
surf(X,Y,Z);
axis square;
view(0,90);

light() 函數：

light Create light.

打光的角度

%% light
x = linspace(-3,3,20);
y = linspace(-4,4,20);
[X,Y] = meshgrid(x,y);
Z = X.*exp(-X.^2-Y.^2);
subplot(221);
surf(X,Y,Z);
axis square;
xlabel('x');
ylabel('y');
zlabel('z');

subplot(222);
surf(X,Y,Z);
axis square;
light('position',[0,0,0.5]);
hold on;  
plot3(0,0,0.5,'ro');

subplot(223);
surf(X,Y,Z);
axis square;
light('position',[5,5,0.5]);
hold on;  
plot3(5,5,0.5,'ro');

subplot(224);
surf(X,Y,Z);
axis square;
light('position',[-5,-5,0.5]);
hold on;  
plot3(-5,-5,0.5,'ro');

圖中紅點為打光的位置！

patch() 函數

patch 是補丁的意思

A graphical object containing polygons

%% patch()  
x = [0,1,1,0];
y = [0,0,1,1];

subplot(231);
patch(x,y,'red')

x2 = [2,5; 2,5; 8,8];
y2 = [4,0; 8,2; 4,0];
subplot(232);
patch(x2,y2,'green')


v = [0, 0; 1 ,0; 1, 1; 0, 1];
f = [1,2,3];  % 指定要上顏色的面(3 個點)
subplot(233);
patch('faces',f,'vertices',v,'facecolor','red')


v = [0 0 0; 1 0 0 ; 1 1 0; 0 1 0; 0.25 0.25 1; ...
0.75 0.25 1; 0.75 0.75 1; 0.25 0.75 1]; 
f =  [1 2 3 4; 5 6 7 8; 1 2 6 5; 2 3 7 6; 3 4 8 7; 4 1 5 8]; % 
subplot(234);
patch('faces',f,'vertices',v, 'facecolor', 'cyan');  
view(3); %3-d view 

v = [0 0 0; 1 0 0 ; 1 1 0; 0 1 0; 0.25 0.25 1; ...
0.75 0.25 1; 0.75 0.75 1; 0.25 0.75 1]; 
f =  [1 2 3 4; 5 6 7 8;  2 3 7 6; 3 4 8 7; 4 1 5 8]; % 沒有這個面 1 2 6 5;
subplot(235);
patch('faces',f,'vertices',v, 'facecolor', 'cyan');  
view(3); %3-d view

lecture 07 GUI 界面

如何開啟？

設置顯示組件的名稱：

有疑問自行谷歌，（會QT的人感覺到很舒服）

對齊工具：

修改文字內容：

.fig 對應的 .m 文件：

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
x = 1:0.1:10;
y = sin(x);
plot(x,y);





% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
x = 1:0.1:10;
y = cos(x);
plot(x,y);



% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
x = 1:0.1:10;
y = tan(x);
plot(x,y);

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
x = 1:0.1:10;
y = sin(x);
plot(handles.axes1,x,y); % 指定在 axes1 上作圖



% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
x = 1:0.1:10;
y = cos(x);
axes(handles.axes4); % 指定在 axes4 上作圖
plot(x,y);

% --- Executes on slider movement.
function slider1_Callback(hObject, eventdata, handles)
% hObject    handle to slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
val1 = get(hObject,'Value');
val2 = get(handles.slider2,'Value');
set(handles.res_text,'String',num2str(int16(val1+val2)));  % key -value 


% --- Executes on slider movement.
function slider2_Callback(hObject, eventdata, handles)
% hObject    handle to slider2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
val1 = get(hObject,'Value');
val2 = get(handles.slider1,'Value');
set(handles.res_text,'String',num2str(int16(val1+val2)));  % key -value

全局變量

編譯GUI程序為可執行文件：

需要使用 deploytool 指令，

選擇第一個，

Lecture 08 Image Processing 圖像處理 I：

Introduction to digital image

Types of Digital Image

三種具體的Image

Read and show images

Read and Show An Image

clear;
close all; % 關閉所有的figure

I = imread('pout.tif');% 圖片為matlab 內置的圖片
imshow(I); % show

我們可以通過改變這些數值來修改圖片，

例如如下程序：

clear;
close all; % 關閉所有的figure 

I = imread('pout.tif');% 圖片為matlab 內置的圖片 
subplot(121);
imshow(I); % show 

% 對圖像的數據進行操作
for i = 1:size(I,1)%size(I)[1]
    for j = 1:size(I,2)%size(I)[2]
        if mod(i,2) == 0 && mod(j,2)  == 0
            I(i,j) = 0; % 0 是黑色的  
        end
    end
end

subplot(122);
imshow(I); % show

imageinfo（）函數

它可以查看圖像的信息

%% imageinfo() 
%I = imread('pout.tif');% 圖片為matlab 內置的圖片 
imageinfo('pout.tif'); % 獲取文件的相關信息

從中我們也可以看出 pout.tif 是內置的圖像，可以通過路徑把它找出來，

imtool() 函數

%% imtool() 
imtool('pout.tif'); % isplay image in the Image Tool.

Image arithmetic（Image 的四則運算）

Image Multiplication: immultiply()

%% immultiply() image 乘法 
I = imread('rice.png');% 內置的圖片
subplot(121);
imshow(I);

J = immultiply(I,1.5); % 其實是將像素值 都乘以1.5    
subplot(122);
imshow(J);

Image Addition: imadd()

%% imadd() image 加法 
I = imread('rice.png');% 內置的圖片
J = imread('cameraman.tif');% 內置
K = imadd(I,J);

subplot(131);
imshow(I);

subplot(132);
imshow(J);

subplot(133);
imshow(K);

注意，相加的時候，這兩個圖片的 dim 要相同。

Image Histogram: imhist() （直方圖）

isplay histogram of image data.

%% imhist() isplay histogram of image data.
I = imread('pout.tif');% gray image
subplot(121);
imshow(I);
subplot(122);
imhist(I);

注：一般都是展示灰度圖的直方圖，

Histogram Equalization: histeq() 直方圖均衡化

Enhances the contrast of the image ，增強圖像的對比度

那么為什么要做圖像均衡化，

把集中在小區間的像素值變為大區間的！

Geometric Transformation 幾何變換

這和前面不一樣的是，這里改變的是 像素的位置， 前面改變的像素值。

一遍是通過矩陣進行改變：

左乘

更方便的是使用已經提供的內置函數，

例如做旋轉：

Image Rotation: imrotate()

%% Image rotate 
I = imread('rice.png'); 
subplot(121);
imshow(I); 
J = imrotate(I, 35, 'bilinear'); % 雙線性插值 參見： https://www.cnblogs.com/linkr/p/3630902.html
subplot(122); 
imshow(J);
size(I)
size(J)

注：看起來兩個圖片的大小一樣，其實是不一樣的，第二個的要比第一個大一些，

練習：不使用 imrotate() 實現圖片的旋轉，要求使用旋轉矩陣

jiaodu=45;                       %要旋轉的角度，旋轉方向為順時針
img=imread('rice.png');       %這里v為原圖像的高度，u為原圖像的寬度
subplot(121)
imshow(img);                    %這里y為變換后圖像的高度，x為變換后圖像的寬度
[h w]=size(img);

theta=jiaodu/180*pi;%轉換為弧度
rot=[cos(theta) -sin(theta) 0;
    sin(theta) cos(theta) 0;
    0 0 1];
pix1=[1 1 1]*rot;               %變換后圖像左上點的坐標
pix2=[1 w 1]*rot;               %變換后圖像右上點的坐標
pix3=[h 1 1]*rot;               %變換后圖像左下點的坐標
pix4=[h w 1]*rot;               %變換后圖像右下點的坐標

height=round(max([abs(pix1(1)-pix4(1))+0.5 abs(pix2(1)-pix3(1))+0.5]));     %變換后圖像的高度
width=round(max([abs(pix1(2)-pix4(2))+0.5 abs(pix2(2)-pix3(2))+0.5]));      %變換后圖像的寬度
imgn=zeros(height,width);

delta_y=abs(min([pix1(1) pix2(1) pix3(1) pix4(1)]));            %取得y方向的負軸超出的偏移量
delta_x=abs(min([pix1(2) pix2(2) pix3(2) pix4(2)]));            %取得x方向的負軸超出的偏移量

for i=1-delta_y:height-delta_y
    for j=1-delta_x:width-delta_x
        pix=[i j 1]/rot;                                %用變換后圖像的點的坐標去尋找原圖像點的坐標，
        %否則有些變換后的圖像的像素點無法完全填充
        float_Y=pix(1)-floor(pix(1));
        float_X=pix(2)-floor(pix(2));
        
        if pix(1)>=1 && pix(2)>=1 && pix(1) <= h && pix(2) <= w
            
            pix_up_left=[floor(pix(1)) floor(pix(2))];          %四個相鄰的點
            pix_up_right=[floor(pix(1)) ceil(pix(2))];
            pix_down_left=[ceil(pix(1)) floor(pix(2))];
            pix_down_right=[ceil(pix(1)) ceil(pix(2))];
            
            value_up_left=(1-float_X)*(1-float_Y);              %計算臨近四個點的權重
            value_up_right=float_X*(1-float_Y);
            value_down_left=(1-float_X)*float_Y;
            value_down_right=float_X*float_Y;
            
            imgn(i+delta_y,j+delta_x)=value_up_left*img(pix_up_left(1),pix_up_left(2))+ ...
                value_up_right*img(pix_up_right(1),pix_up_right(2))+ ...
                value_down_left*img(pix_down_left(1),pix_down_left(2))+ ...
                value_down_right*img(pix_down_right(1),pix_down_right(2));
        end
        
    end
end
subplot(122)
imshow(uint8(imgn))

參考代碼

Write Image: imwrite()

Lecture 09 Image Processing 圖像處理 II：

問題的提出：

rice.png 中有多少顆米，以及每粒米的大小？

Image thresholding （圖像閾值化）

A gray-level image can be turned into a binary image by using a threshold

graythresh() and im2bw()

graythresh() 用來計算最優的threshold 值

im2bw() Convert image to binary image by thresholding. b and w in im2bw() is b: black 0 w: white 1

%% 對一個灰度圖進行 threshold 
I = imread('rice.png');
subplot(131);
imshow(I);
% 首先找出一個最優的threshold level 值,這里使用 graythresh() 
threshold_level = graythresh(I);%LEVEL is a normalized intensity value that lies in the range [0, 1]

subplot(132);
imhist(I); % 直方圖
threshold_val = threshold_level*255; % 閾值  
ylim=get(gca,'Ylim'); % 獲取當前圖形的縱軸的范圍
hold on;
plot([threshold_val,threshold_val],ylim);   

subplot(133);
bw = im2bw(I,threshold_level);
imshow(bw);

產生的問題：

1，在 Backgroud 上存在一些白色的點（這是因為圖片上面比較亮，導致背景被誤判前景 ）

2，一些米不見了（這是因為圖片下面比較暗，導致前景被誤判背景 ）

如何解決這些問題呢？

1，要Estimate the background

2, 要Substract the background

Background image estimation （背景估計）

Estimation for the gray level of the background

%% background estimation
I = imread('rice.png');
subplot(131);
imshow(I);
BG = imopen(I,strel('disk',15));
subplot(132);
imshow(BG);

% BG in 3-d view
x = 1:8:size(BG,1);
y = 1:8:size(BG,2);
[X,Y] = meshgrid(x,y);
Z = BG(x,y);
ax3 = subplot(133);
surf(X,Y,Z); % 3d 圖
colormap(ax3,'jet');
set(gca,'xlim',[1,256])% image is 256x256
set(gca,'ylim',[1,256])
set(gca,'zlim',[0,255])%gray value is [0,255]
xlabel('x'); ylabel('y');zlabel('z');
axis square;

可以看出背景的上面值高，下面值低。

Background image subtraction（背景剔除）

%% background 剔除 
I = imread('rice.png');
subplot(131);
imshow(I);
BG = imopen(I,strel('disk',15)); 
subplot(132);
imshow(BG);
I2 = I - BG;  % 這時 整體的圖像就比較均勻了，就改善了 亮度 不統一的問題
subplot(133);
imshow(I2);

此時我們再次進行 threshold 就會好很多

%% thresholding on Background Removed Image 
I = imread('rice.png');
subplot(221);
imshow(I);

subplot(222);
threshold_level = graythresh(I);
bw = im2bw(I,threshold_level);
imshow(bw);

subplot(223);
imshow(I);

subplot(224);
% 去除背景
BG = imopen(I,strel('disk',15)); 
I2 = I - BG;

threshold_level2 = graythresh(I2);
bw2 = im2bw(I2,threshold_level2);
imshow(bw2);

Component-connected labeling （連接元素的標定）

現在怎么算呢？可以用 Component-connected labeling 算法，

Matlab 為我們提供了一個函數 bwlabel()

%% bwlabel()  
I=imread('rice.png');
BG=imopen(I, strel('disk', 15));
I2=imsubtract(I, BG); 
threshold_level=graythresh(I2);
BW=im2bw(I2, threshold_level);
[labeled, numObjects]=bwlabel(BW, 8);% 參數 8是8連通（上下左右加四角）。也可以是4，4連通（上下左右）
                                                            % labeled
                                                            % 是標定的矩陣，它的大小和BW一致
                                                            %  numObjects是總數  
disp(numObjects)

現在我們知道了米粒的總數為 99 ，

What is the size of the largest grain?

What is the mean size of the grains？

可以通過兩個循環來統計 labeled 這個矩陣中的個數。

%% What is the size of the largest grain?  What is the mean size of the grains?
I = imread('rice.png');
BG = imopen(I,strel('disk',15));
I2 = imsubtract(I,BG);
level = graythresh(I2);
BW = im2bw(I2,level);
[labeled,numObjects] = bwlabel(BW,8);
[heigh,width] = size(labeled);

largest = 0;%用於存儲最大米粒的像素個數
sum = 0;%用於計算平均數前的總數
for num = 1:1:numObjects %用於迭代米粒的個數標號
    calc = 0;  %初始化每粒米個數
    for i = 1:1:heigh
        for j = 1:1:width
            if labeled(i,j) == num
                calc = calc + 1;
            end
        end
    end
    if calc > largest
        largest = calc;
    end
    sum = calc + sum;
end
disp(largest);
average = sum/numObjects;
disp(average);

參考代碼

Note：參考代碼都是網上隨便找的，可能質量不是很好

label2rgb() 函數

可以用它來給米粒標上顏色

%% label2rgb()
I=imread('rice.png');
BG=imopen(I, strel('disk', 15));
I2=imsubtract(I, BG); 
threshold_level=graythresh(I2);
BW=im2bw(I2, threshold_level);
[labeled, numObjects]=bwlabel(BW, 8);
RGB_label=label2rgb(labeled); imshow(RGB_label);

小練習：

如何變為color 圖的紅米呢？具體做法是在 r層標紅

I=imread('rice.png');
BG=imopen(I, strel('disk', 15));
I2=imsubtract(I, BG); level=graythresh(I2);
BW=im2bw(I2, level);
[labeled, numObjects]=bwlabel(BW, 8);
No=zeros(1,99);
for i=1:256;
    for j=1:256;
        for n=1:numObjects;
            if labeled(i,j)==n;
                No(n)=1+No(n);
            else
            end
        end
    end
end
max(No)
mean(No)
perfectnum=0;
for k=1:99;
    if No(k)<80;
        perfectnum=perfectnum+0;
    else if No(k)>350;
            perfectnum=perfectnum+2;
        else perfectnum=perfectnum+1;
        end
    end
end
disp(perfectnum);
%以上為顯示優化計算后的米的數量，下面顯示柱狀圖與紅米；
subplot(121);histogram(No);
RGB_label = label2rgb(BW,'flag','k');
subplot(122);imshow(RGB_label);

參考代碼

Note：參考代碼都是網上隨便找的，可能質量不是很好

regionprops() 函數：

%% regionprops() 
I=imread('rice.png');
BG=imopen(I, strel('disk', 15));
I2=imsubtract(I, BG); 
threshold_level=graythresh(I2);
BW=im2bw(I2, threshold_level);
[labeled, numObjects]=bwlabel(BW, 8);
graindata = regionprops(labeled, 'basic');
disp(graindata);

OUTPUT：
99x1 struct array with fields:

    Area
    Centroid
    BoundingBox

Interactive Selection: bwselect()

Lets you select objects using the mouse

%% bwselect
I=imread('rice.png');
BG=imopen(I, strel('disk', 15));
I2=imsubtract(I, BG); 
threshold_level=graythresh(I2);
BW=im2bw(I2, threshold_level);
[labeled, numObjects]=bwlabel(BW, 8);

ObjI = bwselect(BW); 
imshow(ObjI);  % 用鼠標 左鍵選定后， 按 右鍵

Lecture 10 Differentiation and Integration（微分和積分）：

微分的定義：

微分：differentitation 是 △ 的意思。可以簡單理解為對其在某一點求導。

Polynomial differentiation（多項式的微分）

polyval() 函數

%% polyval()   計算x 在 f多項式下的值  
f = [9,-5,3,7]; % 多項式的系數  
x = -2:1:5; % x 值 
y = polyval(f,x);% y  值 。它按照 y = 9x^3 -5x^2+3x +7 計算y值   
plot(x,y);

polyder() 函數

%% polyder()  求一個多項式的導數（微分）
f = [9,-5,3,7]; % 多項式的系數  
f_d = polyder(f);% f' 多項式的系數，即函數的導數 
disp(f_d);

OUTPUT
27   -10     3

小練習：給定 𝑓(𝑥) =9x^3 -5x^2+3x +7，求𝑓′(7)?

%% f(x) = 9x^3 -5x^2+3x +7 求 f'(7)
f_d = polyder([9,-5,3,7]);  % f'(x) 的系數  
ret = polyval(f_d,[7,]); % 求f'(x) 這個多項式在 x = 7 的值 
disp(ret);

OUTPUT 
1256

小練習：

%% conv() 函數   Convolution and polynomial multiplication. 卷積與多項式乘法 （一維卷積）
% 所謂兩個向量卷積，說白了就是多項式乘法
f1 = [5,-7,5,10]; %  第一個多項式的系數
f2 = [4,12,-3]; %  第二個多項式的系數
f = conv(f1,f2); % 一維卷積 f1*f2 之后 多項式的系數 

x = -2:0.1:1; % x值 
y = polyval(f,x); % y值
plot(x,y,'b--','linewidth',2);
hold on;

f_d = polyder(f);% f'(x) 的系數 
y2 = polyval(f_d,x); % y' 值 
plot(x,y2,'r-','linewidth',2);

legend('f(x)','f ''(x)');  % 輸出單引號使用 轉義字符無效，使用兩個單引號即可

這里涉及到了卷積【擴展】，

參考文章： https://www.cnblogs.com/hyb221512/p/9276621.html 和 https://blog.csdn.net/leviopku/article/details/80327478

一維卷積 conv : 簡單記就是兩個多項式相乘

二維卷積 conv2 或多維卷積 convn : 分三種 full , same ,valid

以二維卷積展示三種模式：

full mode:

same mode:

valid mode:

%% conv2() 函數  二維卷積 
X = ones([3,3]); % 3x3 
K = ones([3,3]); % 3x3 

disp('======full=======');  
ret = conv2(X,K); % default is full mode 
disp(ret);

disp('======same=======');  
ret2 = conv2(X,K,'same'); % default is full mode 
disp(ret2);

disp('======valid=======');  
ret3 = conv2(X,K,'valid'); % default is full mode 
disp(ret3);

OUTPUT 
======full=======
     1     2     3     2     1
     2     4     6     4     2
     3     6     9     6     3
     2     4     6     4     2
     1     2     3     2     1

======same=======
     4     6     4
     6     9     6
     4     6     4

======valid=======
     9

關於卷積輸出的大小計算公式：

1.full: 滑動步長為1，圖片大小為N1xN1，卷積核大小為N2xN2，卷積后圖像大小：N1+N2-1 x N1+N2-1

如滑動步長為1，圖片大小為3x3，卷積核大小為3x3，卷積后圖像大小：5x5

2.same: 滑動步長為1，圖片大小為N1xN1，卷積核大小為N2xN2，卷積后圖像大小：N1xN1

3.valid:滑動步長為S，圖片大小為N1xN1，卷積核大小為N2xN2，卷積后圖像大小：(N1-N2)/S+1 x (N1-N2)/S+1

如滑動步長為1，圖片大小為3x3，卷積核大小為3x3，卷積后圖像大小：1x1

Polynomial Integration（多項式的積分）

polyint() 函數

%% polyint() 函數  
f = [5,0,-2,0,1]; % f(x) 多項式的系數 
f_int = polyint(f,3); % 常數項 constant 設置為3  
disp(f_int);

% 求 ∫ f(7)?
ret = polyval(f_int,7);
disp(ret);

Numerical differentiation（數值的微分）

diff() 函數：

%% diff() 函數 
x = [1,2,5,2,1];
ret = diff(x); %calculates the differences between adjacent elements of a vector。
                                             % 相鄰的  
disp(ret);

小練習：

%% Exercise: obtain the slope of a line between 2 points (1,5) and (2,7)
             % 獲取      斜率 
x = [1,2];
y = [5,7];
slope = diff(y)/diff(x); % 斜率
disp(slope);

%%  
x0 = pi/2;
h = 0.1;
x = [x0,x0+h];
y = sin(x);
res = diff(y)/diff(x);
disp(res);

How does ℎ affect accuracy?

%% 
hs = [0.1,0.01,0.001,0.0001,0.00001,0.000001,0.0000001];
% f = @(x) x+1;

x0 = pi/2;
f = @(h) diff(sin([x0,x0+h]))/diff([x0,x0+h]);
for h = hs
    ret = f(h) - cos(pi/2); % cos(pi/2) 是精確值 
    disp(ret);
end

OUTPUT （Error Value）
   -0.0500

   -0.0050

  -5.0000e-04

  -5.0000e-05

  -5.0000e-06

  -5.0004e-07

  -4.9960e-08

How to Find the 𝑓′ over An Interval [0,2𝜋]?

%%  前面是一個點的導數,這里是一個區間的導數 
h = 0.5;
x  = 0:h:2*pi;
y = sin(x);
ret = diff(y)./diff(x); % 多個值 除就要 加 . 了 ！  
disp(ret);

hold on;
plot(x,y);
plot(x(1:size(x,2)-1),ret,'ro-'); % x 要少一個 

legend('sin(x)','sin''(x)','location','best'); % 會自動選擇一個最佳的位置

Second and Third Derivatives

%% 二次 和 三次 微分(導數)  
h = 0.5;
x = 0:h:2*pi;
%y = x.^3;
y = sin(x);
y_1 = diff(y)./diff(x);
y_2 = diff(y_1)./diff(x(1:end-1));
y_3 = diff(y_2)./diff(x(1:end-2));
hold on;
plot(x,y,'r.-');
plot(x(1:end-1),y_1,'g*-');
plot(x(1:end-2),y_2,'b+-');
plot(x(1:end-3),y_3,'cx-');

legend('sin(x)','sin''(x)','sin''''(x)','sin''''''(x)','location','best');

Numerical and integration （數值的積分）

第一個是用矩形來做預測，第二個是用梯形來做預測。

Midpoint Rule Using sum()

%%  ∫ 4x^3 dx [0,2]
h = 0.05; 
x = 0:h:2;
midpoint = ( x(1:end-1) + x(2:end) ) ./2;
s = 4*midpoint.^3 * h;
A =  sum(s);
disp(A);  %    15.9950

Trapezoid Rule Using sum()

%%  ∫ 4x^3 dx [0,2]  trapezoid  
f = @(x) 4*x.^3; % x is a vector 
h = 0.05; 
x = 0:h:2;
x1 = x(1:end-1);
x2 = x(2:end);

s = ( f(x1) + f(x2) ) /2 * h;  
A = sum(s);
disp(A); %   16.0100

也可以使用trapz() 函數，

h = 0.05; 
x = 0:h:2;
trapezoid = trapz(4*x.^3);% 它接收的參數是 每個點的y值  
                          % 返回的是 sum ( f(x0)+f(x1) ) /2
A = trapezoid*h;  
disp(A); %    16.0100

Second-order Rule: 1/3 Simpson’s

第三種積分方法，它一次積兩個的面積。所以它要求等分的數目一定是偶數個！

%% 1/3 simpson's rule 
h = 0.05; 
x = 0:h:2; %  必須是 EVEN 段 ,也就是 EVEN + 1 個點。而且要取到 0 和 2  
y = 4*x.^3;
A = h/3*(y(1) +  4*sum(y(2:2:end)) + 2*sum(y(3:2:end-2))  + y(end));  
disp(A); % 16 精准！

Numerical Integration: integral() 函數

%% integral() 函數 
y = @(x) sin(x);  % 匿名函數  其實y 是個函數指針 
res = integral(y,0,pi);
disp(res); % 2.0000 

y2 = @(x) x.^2;  % 當然也可以對 多項式來積分  
res2 = integral(y2,0,1);
disp(res2); %     0.3333

Double and Triple Integrals （二重三重積分）

%% 二重/ 三重積分  
f1 = @(x,y) y.*sin(x) + x.*cos(x);
res1 = integral2(f1,pi,2*pi,0,pi);
disp(res1); % -3.5864


f2 = @(x,y,z) y.*sin(x) + z.*cos(y);
res2 = integral3(f2,0,pi,0,1,-1,1);
disp(res2); % 2.0000

Lecture11 Root Finding （方程式求根）

Problem Statement

Symbolic approach （符號逼近）（解析解）

Symbolic Root Finding Approach

%%  1 使用 syms 或者 sym()  

syms x; % 此時 x 就是一個符號了  
ret = x + x + x;
disp(ret);

ret2 = ret/4;
disp(ret2);
%OUTPUT:
%3*x
%(3*x)/4

sym('x'); % 此時 x 就是一個符號了  
ret = x + x + x;
disp(ret);

ret2 = ret/4;
disp(ret2);
%OUTPUT:
%3*x
%(3*x)/4

%% to define y = x^2 - 2x -8 
x = sym('x');
y = x^2 - 2*x -8; % y 也自動是個 sym 了  
disp(y);

Symbolic Root Finding: solve()

%% solve()  y = x*sin(x) -x = 0

syms x; % x 是個符號 
res = solve('x*sin(x) - x',x);
disp(res); % res 也是個符號 

% 也可以如下 
syms x;
y = x*sin(x) -x; % y 也是個符號
res2 = solve(y,x); 
disp(res2);

%OUTPUT
%    0
% pi/2
% 
%    0
% pi/2

Solving Multiple Equations

%% solve multiple equations 
syms x y;
eq1 = x - 2*y - 5;
eq2 = x +  y  - 6;
res = solve(eq1,eq2,x,y);
disp(res); % res 是個 結構體 
disp(res.x); 
disp(res.y); 

%OUTPUT:
%    x: [1x1 sym]
%    y: [1x1 sym]
%
%    x: [1x1 sym]
%    y: [1x1 sym]
%
%17/3
% 
%1/3

Solving Equations Expressed in Symbols

%% 
syms a b x;
y = a*x^2 - b;
res1 = solve(y); % 默認以 x 為 未知數  
disp(res1);


res2 = solve(y,a); % 以 a 為 未知數  
disp(res2);


res3 = solve(y,b); % 以 b 為 未知數  
disp(res3);

%OUTPUT
%  b^(1/2)/a^(1/2)
% -b^(1/2)/a^(1/2)
% 
%b/x^2
% 
%a*x^2

%% Exercise   
% 1 
syms a b x y r;
eq = (x-a)^2 + (y-b)^2 - r^2;
res = solve(eq,x);
disp(res);

%OUTPUT 
% a + (b + r - y)^(1/2)*(r - b + y)^(1/2)
% a - (b + r - y)^(1/2)*(r - b + y)^(1/2)


% 2 
syms a b c d;
m = [a b;c d];
res2 = inv(m); % inv() Matrix inverse
disp(res2);
disp(m*res2);

%OUTPUT
%[  d/(a*d - b*c), -b/(a*d - b*c)]
%[ -c/(a*d - b*c),  a/(a*d - b*c)]
%  
%[ (a*d)/(a*d - b*c) - (b*c)/(a*d - b*c),                                     0] % 沒有化簡 
%[                                     0, (a*d)/(a*d - b*c) - (b*c)/(a*d - b*c)]

Symbolic Differentiation and Integration （符號微分和積分）

%% symbolic differentitation  diff()  
syms x;
y = 4*x^5;
res = diff(y); 
disp(res);  

%OUTPUT:
%20*x^4

%% Exercise
%1
syms x y;
f = exp(x^2) /(x^3 - x + 3);
res = diff(f,x);
disp(res);

%2
g = (x^2 + x*y - 1 )/(y^3 + x+ 3);
res2 = diff(g,x);
disp(res2);

%OUTPUT:
%(2*x*exp(x^2))/(x^3 - x + 3) - (exp(x^2)*(3*x^2 - 1))/(x^3 - x + 3)^2
%
%(2*x + y)/(y^3 + x + 3) - (x^2 + y*x - 1)/(y^3 + x + 3)^2

%% symbolic integration  int()  
syms x;
y = x^2*exp(x);
z = int(y);
disp(z);
z = z-subs(z, x, 0); % subs 是 替代的意思，在z 中，x用0取代
disp(z);

%OUTPUT 
%exp(x)*(x^2 - 2*x + 2)
% 
%exp(x)*(x^2 - 2*x + 2) - 2

%% Exercise 
syms x;
f = (x^2 - x+1) /(x+3);
F = int(f);
disp(F);
res = subs(F,x,10) - subs(F,x,0);
disp(res);

%OUTPUT:
%13*log(x + 3) - 4*x + x^2/2
%13*log(13) - 13*log(3) + 10

Numeric root solvers （數值根求解）

Symbolic vs. Numeric

fsolve() 函數

%% fsolve() 函數  
f = @(x) 1.2*x + 0.3 + x*sin(x);
ret = fsolve(f,0); % 初始根  （猜測一個數）
disp(ret);

%OUTPUT
%Equation solved.
%
%fsolve completed because the vector of function values is near zero
%as measured by the default value of the function tolerance, and
%the problem appears regular as measured by the gradient.
%
%<stopping criteria details>
%
%   -0.3500

%% Exercise 
f = @(x) [2*x(1) - x(2) - exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];
res = fsolve(f,[-5,-5]);
disp(res);

%OUTPUT
%Equation solved.
%
%fsolve completed because the vector of function values is near zero
%as measured by the default value of the function tolerance, and
%the problem appears regular as measured by the gradient.
%
%<stopping criteria details>
%
%    0.5671    0.5671

fzero() 函數

%% fzero() 函數  
f = @(x) x.^2;
res1 = fzero(f,0.2); % guess value is 0.2 
disp(res1); 

res2 = fsolve(f,0.2); 
disp(res2);

%OUTPUT:
%    NaN 因為沒有 cross the x axis 
%
%    0.0063

fsolve() fzero() 的options

%% fsolve() fzero() 的options 
f = @(x) x.^2;
                                  %tolerance 偏差 
options = optimset('MaxIter',1e3,'TolFun',1e-10);
res = fsolve(f,0.1,options);
disp(res);
res2 = fzero(f,0.1,options);
disp(res2);

%OUTPUT:
%Equation solved.
%
%fsolve completed because the vector of function values is near zero
%as measured by the selected value of the function tolerance, and
%the problem appears regular as measured by the gradient.
%
%<stopping criteria details>
%
%   1.9532e-04
%
%Exiting fzero: aborting search for an interval containing a sign change
%    because NaN or Inf function value encountered during search.
%(Function value at -1.37296e+154 is Inf.)
%Check function or try again with a different starting value.
%   NaN

%% 
f = @(x) x.^2 - 2; % cross the x axis 

res3 = fzero(f,1.4); % guess value is 1.4 
disp(res3); 

res4 = fsolve(f,1.4);
disp(res4); 

res3 = fzero(f,-1.4); % guess value is -1.4 
disp(res3); 

res4 = fsolve(f,-1.4);
disp(res4); 

%OUTPUT:
%
%    1.4142
%
%
%Equation solved.
%
%fsolve completed because the vector of function values is near zero
%as measured by the default value of the function tolerance, and
%the problem appears regular as measured by the gradient.
%
%<stopping criteria details>
%
%    1.4142
%
%   -1.4142
%
%
%Equation solved.
%
%fsolve completed because the vector of function values is near zero
%as measured by the default value of the function tolerance, and
%the problem appears regular as measured by the gradient.
%
%<stopping criteria details>
%
%   -1.4142

Finding Roots of Polynomials: roots()

%% roots() 函數  only works for polynomials 

res = roots([1,-3.5,2.75,2.125,-3.875,1.25]);
disp(res);

%OUTPUT:
%   2.0000 + 0.0000i
%  -1.0000 + 0.0000i
%   1.0000 + 0.5000i
%   1.0000 - 0.5000i
%   0.5000 + 0.0000i

How Do These Solvers Find the Roots?（上面的fsolve 和 fzero 都是如何找到根的呢？）

第一個方法： Bisection Method（Bracketing）（二分法）

l : lower u:upper

第二個方法： Newton-Raphson Method (Open)

不收斂的情況：

The tangent lines of x^3 − 2x + 2 at 0 and 1 intersect the x-axis at 1 and 0 respectively, illustrating why Newton's method oscillates between these values for some starting points.

它一直在 1 和 0 之間震盪，這種時候應該采用二分法來解決。

Lecture 12 Linear equations and Linear system(線性方程和線性系統)：

Linear Equation（線性方程）

Why Matrix Form?

更方便計算，也更簡潔

如何解線性方程組？

常用的兩種方法：

1. Successive elimination (through factorization) 逐次消除（通過因子分解）

又可以分為 Gaussian Elimination 和 LU Factorization

a Gaussian Elimination

MATLAB 中使用 Gaussian Elimination – rref()函數

%% Gaussian Elimination rref()  Reduced row echelon form.
A = [1,2,1;2,6,1;1,1,4];
b = [2;7;3];
res = rref([A,b]);
disp(res);

%OUTPUT:
%     1     0     0    -3
%     0     1     0     2
%     0     0     1     1

b LU Factorization lu() 函數

參考里面的LU分解部分： https://www.cnblogs.com/zach0812/p/13368110.html

%% LU Factorization
A = [1,1,1;2,3,5;4,6,8];
[L,U,P] = lu(A); % permutation(交換) matrix P so that P*A = L*U
                        % 交換矩陣 實現行交換  (因為主元位置出現0的情況需要進行 行交換 )
                        
disp(L);
disp(U);
disp(P);
disp('====');
disp(P*A);
disp(L*U);

%OUTPUT：
%    1.0000         0         0
%    0.2500    1.0000         0
%    0.5000         0    1.0000
%
%    4.0000    6.0000    8.0000
%         0   -0.5000   -1.0000
%         0         0    1.0000
%
%     0     0     1
%     1     0     0
%     0     1     0
%
%====
%     4     6     8
%     1     1     1
%     2     3     5
%
%     4     6     8
%     1     1     1
%     2     3     5

Matrix Left Division: \ or mldivide()

Ax = b => x = A \ b;

Matlab 的左除 \ 原理是使用 Successive elimination 來計算，不是使用 Cramer’s method （克拉默法則）

2. Cramer’s method （克拉默法則）Ax = b , x =^-1*b

矩陣C是“代數余子式矩陣”（cofactor matrix）

C 的轉置矩陣 C^T。即伴隨矩陣（adjoint matrix）。參考： https://www.cnblogs.com/zach0812/p/13368110.html

克萊姆法則 Cramer’s Rule x = A^(-1)b

現在我們已經知道了 A的求逆公式，那么對於Ax = b ,那么

Example:

𝑥 + 2𝑦 + 𝑧 = 2

2𝑥 + 6𝑦 + 𝑧 = 7

𝑥 + 𝑦 + 4𝑧 = 3

%% Cramer's method  inv() 函數  
A = [1,2,1;2,6,1;1,1,4];
b = [2;7;3];
x = inv(A)*b; 
disp(x);

%OUTPUT:
%   -3.0000
%    2.0000
%    1.0000

Cramer’s Method 的問題：

A的inverse 不一定存在，此時 A矩陣就是一個 singular （奇異）矩陣。

下面哪個矩陣更健康？

肯定是B，

%% cond() 函數，來計算 一個矩陣 k值  越小越好 
A = [ 1 2 3; 2 4.0001 6; 9 8 7]; 
cond(A)
B = [ 1 2 3; 2 5 6; 9 8 7]; 
cond(B)

%OUTPUT
%ans =
%   4.3483e+05
%ans =
%   45.5623

Linear System （線性系統）

Eigenvalues and Eigenvectors 特征值和特征向量

作用一：簡化Ab 的計算，將矩陣向量的乘法轉化為向量的加法

作用二：闡述系統A的放大，縮小作用

A 是一個系統，，

當 λ > 1 的時候就是系統 A就是放大的作用

當 λ < 1 的時候就是系統 A就是縮小的作用

eig() 函數 Find the eigenvalues and eigenvectors

%% eig() 函數 
A = [2,-12;1,-5];
[vec,diag] = eig(A);

disp(vec);
disp(diag);

%OUTPUT:
%    0.9701    0.9487  % 兩個特征向量 [0.9701;0.2425]
%    0.2425    0.3162               [0.9487;0.3162]
%
%    -1     0          % 兩個特征值    -1  
%     0    -2                        -2

Lecture 13 Statistics and Data_Analysis（統計與數據分析）

Main Statistical Methodologies

Descriptive Statistics（描述統計學）

Mean, Median, Mode, and Quartile

%% mean() 函數 求平均數  
d = [1,3,5,5,5,5,7,9,9,9,10,13,14];
res = mean(d);
disp(res);

%OUTPUT:
%    7.3077



%% median() 函數 求中位數  
d = [1,3,5,5,5,5,7,9,9,9,10,13,14];
res = median(d);
disp(res);
%OUTPUT:
%    7



%% mode() 函數 Most frequent values in array  
d = [1,3,5,5,5,5,7,9,9,9,10,13,14];
res = mode(d);
disp(res);
%OUTPUT:
%    5

關於四分位數：

四分位數（英語：Quartile）是統計學中分位數的一種，即把所有數值由小到大排列並分成四等份，處於三個分割點位置的數值就是四分位數。

第一四分位數（ $Q_{1}$
第二四分位數（ $Q_{2}$
第三四分位數（ $Q_{3}$

第三四分位數與第一四分位數的差距又稱四分位距（InterQuartile Range, IQR）。

計算方法：

離散分布
對於離散分布，在選擇四分位數值方面沒有普遍的共識。

方法1

使用中位數將有序數據集分為兩半。
1. 　　如果原始排序數據集中的數據點數奇數，請不要在任何一半中包括中位數（排序表中的中心值）。
2. 　　如果原始排序的數據集中有偶數個數據點，則將此數據集精確地分成兩半。
下四分位值是數據下半部分的中位數。高四分位值是數據上半部分的中位數。

方法2

使用中位數將有序數據集分為兩半。
1. 　　如果原始排序數據集中的數據點數量奇數，請在兩個半數中均包括中位數（排序表中的中心值）。
2. 　　如果原始排序的數據集中有偶數個數據點，則將該數據集精確地分成兩半。
下四分位值是數據下半部分的中位數。高四分位值是數據上半部分的中位數。

方法3 （Matlab 中使用）

如果數據點的數量為偶數，則方法3與以上任一方法相同。
如果有（4 n +1）個數據點，則較低的四分位數為第n個數據值的25％加上第（n +1）個數據值的75％；高四分位數是第（3 n +1）個數據點的75％加第（3 n +2）數據點的25％。
如果有（4 n +3）個數據點，則較低的四分位數是第（n +1）個數據值的75％加第（n +2）個數據值的25％；高四分位數是第（3 n +2）個數據點的25％加第（3 n +3）數據點的75％。

方法4

%% prctile() 函數 Percentiles of data set. 百分數
d = [6, 7, 15, 36, 39, 40, 41, 42, 43, 47, 49];
res = prctile(d,50); % 50  
disp(res);

res = prctile(d,25); % 25  
disp(res);

res = prctile(d,75); % 75  
disp(res);

%OUTPUT
%    40
%
%   20.2500
%
%   42.7500

Range and Interquartile Range （范圍和四分位間距）

%% range = max - min
d = [6, 7, 15, 36, 39, 40, 41, 42, 43, 47, 49];
range = max(d) - min(d);
disp(range);


%% interquartile range = q3 - q1
q1 = prctile(d,25); % 25  
q3 = prctile(d,75); % 75  
interquartile_range = q3 - q1;
disp(interquartile_range);

Variance and Standard Deviation （方差和標准差）

%% 方差 
d = [6, 7, 15, 36, 39, 40, 41, 42, 43, 47, 49];
S = var(d);
disp(S);
%OUTPUT:
%  251.9636


%% 標准差  
d = [6, 7, 15, 36, 39, 40, 41, 42, 43, 47, 49];
s = std(d);
disp(s);
%OUTPUT:
%   15.8734

Figures Are Always More Powerful

x = 1:14;
freqy = [1 0 1 0 4 0 1 0 3 1 0 0 1 1];
subplot(1,3,1); bar(x,freqy); xlim([0 15]);
subplot(1,3,2); area(x,freqy); xlim([0 15]);
subplot(1,3,3); stem(x,freqy); xlim([0 15]);

Boxplot 箱形圖

%% 箱形圖 Boxplot 
marks = [80 81 81 84 88 92 92 94 96 97];
boxplot(marks)
res = prctile(marks, [25 50 75]);
disp(res);

%OUTPUT:
%    81    90    94

Skewness（偏度）

%% 偏度 skewness()  
d = randn([100,3]); % 100x3 正態分布
idx= d(:,1) < 0;
d(idx,1) = 0;  % 第一列 小於0 都賦值為0  眾數和最小值都為0  
d(d(:,3)>0,3) = 0;  % 第三列大於0 的都賦值為0  眾數和最大值都為0  
boxplot(d, {'Right-skewed', 'Symmetric', 'Left-skewed'});
res = skewness(d);
disp(res);

%OUTPUT:
%    1.6558    0.2162   -1.9950

注：偏度與平均值和中位數之間的關系不直接相關：偏度為負（左偏）的分布的均值可以大於或小於中位數，正偏度（右偏）也相同。

許多教科書教給了一個經驗法則，指出均值在右偏斜下位於中間值的右側，而在左偏斜下位於中間值的左側。

該規則以令人驚訝的次數失敗。

它可能在多峰分布或尾巴較長而另一尾巴較重的分布中失敗。

例如，在美國家庭中成年居民的分布中，它是個右偏分布。但是，由於大多數情況都小於或等於眾數，均值位於較重的左尾。結果，在右偏的情況下，平均值在中位數的右邊的經驗法則失敗了。

Kurtosis（峰度）

%% kurtosis 峰度
d = randn([100,3])*3; % 100x3 正態分布
k1 = kurtosis(d(:,1));
k2 = kurtosis(d(:,2));
k3 = kurtosis(d(:,3));
disp(k1);
disp(k2);
disp(k3);

%OUTPUT:
% 2.2985
%
% 3.2461
%
% 2.9984

Inferential Statistics（推斷統計學）

置信區間和

假設檢驗

省略這部分！

Lecture 14 Curve Fitting and Interpolation(曲線擬合和插值)

Simple Linear Regression

這里的目的就是要找到兩個 β0 和 β1 使得擬合的效果最好！

Linear Regression Formulation (線性回歸公式)

需要對它做微分！（求導）

注意：這里只有 β0 和 β1 兩個是未知的，

繼續解：

MATLAB中具體可以使用polyfit() 來解，如下：

Polynomial Curve Fitting: polyfit()

%% polyfit() 函數 
x =[-1.2 -0.5 0.3 0.9 1.8 2.6 3.0 3.5];
y =[-15.6 -8.5 2.2 4.5 6.6 8.2 8.9 10.0];
fit = polyfit(x,y,1); %第三個參數是1次的擬合  返回擬合好的兩個參數 ， 
disp(fit); % 將擬合好的結果 畫出來 
xfit = [x(1):0.1:x(end)]; 
yfit = fit(1)*xfit + fit(2); 
plot(x,y,'ro',xfit,yfit); %(x,y) 畫點 (xfit,yfit)畫線 
set(gca,'FontSize',14); 
legend('data points','best-fit');

scatter() 和 corrcoef() 函數

%% scatter() and corrcoef()
x =[-1.2 -0.5 0.3 0.9 1.8 2.6 3.0 3.5];
y =[-15.6 -8.5 2.2 4.5 6.6 8.2 8.9 10.0];
scatter(x,y); % 畫散點圖
box on;
axis square;
r = corrcoef(x,y); % x 和 y 的相關系數
disp(r);

%OUTPUT:
% x y
% x 1.0000(x 和 x) 0.9202 (x 和 y)
% y 0.9202(y 和 x) 1.0000 (y 和 y)

%% corrcoef() 擴展
x = randn(6,1);
y = randn(6,1);
z = randn(6,1);
A = [x,y,z]; % x,y,z 分別是一個列向量
r = corrcoef(A); % x 和 y 和 z 的相關系數
disp(A);
disp(r);

%OUTPUT:
% -0.0450 0.6989 0.6261
% -0.1024 -0.0259 1.1436
% -0.3022 0.3281 -1.8217
% -0.8536 -0.7011 0.2969
% -0.7193 -2.1663 0.1585
% -0.9170 0.0498 0.0461
% x y z
% x 1.0000(xx) 0.5515(xy) 0.1103(xz)
% y 0.5515(yx) 1.0000(yy) -0.1167(yz)
% z 0.1103(zx) -0.1167(zy) 1.0000(zz)

不是越高越好，可能會過擬合，

Multiple Linear regression

上面是一個變量的回歸，下面是多個變量

What If There Exists More Variables?

對兩個變量的線性回歸直觀的表示就要在三維的空間中，

%% multiple linear regression 
load carsmall;
y = MPG;
x1 = Weight; x2 = Horsepower;
X = [ones(length(x1),1) x1 x2];
b = regress(y,X);
x1fit = min(x1):100:max(x1);
x2fit = min(x2):10:max(x2);
[X1FIT,X2FIT]=meshgrid(x1fit,x2fit);
YFIT=b(1)+b(2)*X1FIT+b(3)*X2FIT;
scatter3(x1,x2,y,'filled'); hold on;
mesh(X1FIT,X2FIT,YFIT); hold off;
xlabel('Weight');
ylabel('Horsepower');
zlabel('MPG'); view(50,10);

% [b,bint,r,rint,stats]=regress(y,X);

adding column of ones to make it suitable for matrix multiplication.

小練習：用上面同樣的數據，用下面的公式來擬合

矩陣形式為：

%% multiple linear regression 
load carsmall;
y = MPG;
x1 = Weight; x2 = Horsepower;
X = [ones(length(x1),1) x1 x2 x1.^2 x2.^2 x1.*x2];
b = regress(y,X);
x1fit = min(x1):100:max(x1);
x2fit = min(x2):10:max(x2);
[X1FIT,X2FIT]=meshgrid(x1fit,x2fit);
YFIT=b(1)+b(2)*X1FIT+b(3)*X2FIT + b(4)*X1FIT.^2 + b(5)*X2FIT.^2 + b(6)*X1FIT.*X2FIT;
scatter3(x1,x2,y,'filled'); hold on;
mesh(X1FIT,X2FIT,YFIT); hold off;
xlabel('Weight');
ylabel('Horsepower');
zlabel('MPG'); view(50,10);

% [b,bint,r,rint,stats]=regress(y,X);

上面所涉及的都是線性回歸，如果不是線性的呢？

可以使用cftool()

Interpolation（內插，插值）

回歸的目的是擬合整體的數據

插值的目的是數據不夠，需要估計出新的數據，即插值

In the mathematical field of numerical analysis, interpolation is a type of estimation, a method of constructing new data points within the range of a discrete set of known data points.

In engineering and science, one often has a number of data points, obtained by sampling or experimentation, which represent the values of a function for a limited number of values of the independent variable.

It is often required to interpolate, i.e., estimate the value of that function for an intermediate value of the independent variable.

Common Interpolation Approaches (常見的插值方法)

Linear Interpolation: interp1() 函數

%% interp1()   1-D interpolation 
X = 0:0.5:10; 
V = sin(X); 
hold on;
plot(X,V,'ro-');
Xq = 0:0.3:3; % new Xq 
Vq = interp1(X,V,Xq); 
plot(Xq,Vq,'g*')

Spline Interpolation（樣條插值）: spline()

%% spline()   Cubic spline data interpolation.
X = 0:0.5:10; 
V = sin(X); 
hold on;
plot(X,V,'ro-');
Xq = 0:0.3:3; % new Xq 
Vq = spline(X,V,Xq); 
plot(Xq,Vq,'g*')

What Are Splines?

2D Interpolation: interp2() 函數

%% interp2() 2D Interpolation
x = -2:.5:2; y = -2:.5:3;
[X,Y] = meshgrid(x,y);
Z = X.*exp(-X.^2-Y.^2);
hold on;
surf(X,Y,Z);
plot3(X,Y,Z,'ro');

x_i = 0:0.3:2; % new x
y_i = 0:0.3:3; % new y
[X_i,Y_i] = meshgrid(x_i,y_i);
Z_i = interp2(X,Y,Z,X_i,Y_i); % new Z_i

plot3(X_i,Y_i,Z_i,'g*');

view(3); % 3-d view

2D Interpolation Using Spline

%% 2D Interpolation Using Spline
x = -2:.5:2; y = -2:.5:3;
[X,Y] = meshgrid(x,y);
Z = X.*exp(-X.^2-Y.^2);

x_i = 0:0.3:2; % new x
y_i = 0:0.3:3; % new y
[X_i,Y_i] = meshgrid(x_i,y_i);
Z_i = interp2(X,Y,Z,X_i,Y_i); % default is linear
Z_i_2 = interp2(X,Y,Z,X_i,Y_i,'spline'); % new Z_i

subplot(121);
hold on;
surf(X,Y,Z);
plot3(X,Y,Z,'ro');
plot3(X_i,Y_i,Z_i,'g*-');
view(3);

subplot(122);
hold on;
surf(X,Y,Z);
plot3(X,Y,Z,'ro');
plot3(X_i,Y_i,Z_i_2,'g*-');
view(3);

結束

最后編輯於 2021-01-18 12:40 歷時16天於家

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 台大郭彥甫MATLAB教學-個人筆記（一）《MATLAB從入門到放棄》打通 “矩陣” 障礙放棄MATLAB！簡述winpython為什么比MATLAB更方便 Matlab入門 MATLAB 入門總結 | MATLAB 速成 MATLAB入門教程 MATLAB入門筆記 matlab mex入門簡介 MATLAB入門（一）數組一、MATLAB之基礎入門代碼