# 實驗一 人機介面
## 工作日誌
### 遭遇問題
* 不熟悉Matlab與ASA的溝通方式

### 解決方法
* 詢問助教和同學

## 程式碼
### DataAgent輔助開發狀態空間計算系統
#### ASA

#include "c4mlib.h"
# define K 100
int main()
{
	C4M_DEVICE_set();
	
	/*double A[2][2];
	double B[2][1];
	double C[1][2];
	double D[1][1];
	double x_0[2][1];
	double u[1][K];
	double y[1][K];*/
	
	double A[2][2] = {{1.35 ,0.55}, {-0.45, 0.35}};
	double B[2][1] = {0.5, 0.5};
	double C[1][2] = {3, 1};
	double D[1][1] = {1};
	double x_0[2][1] = {0, 0};
	double u[1][K];
	double y[1][K];
	
	//設定狀態空間方程式的A,B,C,D,x_0和數列u(k)
	/*printf("Please give matrix A\n");
	HMI_snget_matrix(8, 2, 2, &A);
	printf("Please give matrix B\n");
	HMI_snget_matrix(8, 2, 1, &B);
	printf("Please give matrix C\n");
	HMI_snget_matrix(8, 1, 2, &C);
	printf("Please give matrix D\n");
	HMI_snget_matrix(8, 1, 1, &D);
	printf("Please give matrix x_0\n");
	HMI_snget_matrix(8, 2, 1, &x_0);*/
	printf("Please give matrix u\n");
	HMI_snget_matrix(8, 1, K, &u);
	
	//設定x1,x2達到數列x(k)作用
	double x1[2][1];
	double x2[2][1];
	x1[0][0]=x_0[0][0];
	x2[0][0]=x_0[0][0];
	x1[1][0]=x_0[1][0];
	x2[1][0]=x_0[1][0];
	
	//狀態更新方程式x(k+1)=Ax(k)+Bu(k)、y(k)=Cx(k)+Du(k)
	for(int i = 0; i < K; i++){
		x2[0][0]=A[0][0]*x1[0][0]+A[0][1]*x1[1][0]+B[0][0]*u[0][i];
		x2[1][0]=A[1][0]*x1[0][0]+A[1][1]*x1[1][0]+B[1][0]*u[0][i];
		y[0][i]=C[0][0]*x1[0][0]+C[0][1]*x1[1][0]+D[0][0]*u[0][i];
		x1[0][0]=x2[0][0];
		x1[1][0]=x2[1][0];
		printf("%f\n",y[0][i]);
	}
	HMI_snput_matrix(8, 1, K, &y);
	printf("end\n");
}}

#### MATLAB generate_u

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all  % 
clc     %   
close all


f = 10;%頻率
omega = 2*pi*f;
K = 100;%100個點

for i = 0:1:K-1
     u(i + 1) = sin(omega*i/K);
end

u = cast(u,'single');
save datau u;

#### MATLAB plot

clear   %清除記憶體內儲存變數
clc       %清除command window訊息
close all   %清除已開啟的作圖視窗

load datau;
load datay.mat;

k = 0.01:0.01:1;
plot(k , u, 'B-o', k, data0, 'R-x');

xlabel('t');  
ylabel('u(k) , y(k)');  
title('u(k) and y(k)');  
legend('u(k)=sin(wk)','y(k)');

### Matlab Remo put get 巨集輔助開發狀態空間計算系統
#### ASA

#include "c4mlib.h"

int main()
{
	C4M_DEVICE_set();
	float A[2][2] = {{0.0, 0.0}, {0.0, 0.0}};
	float B[2][1] = {{0.0}, {0.0}};
	float C[1][2] = {0.0, 0.0};
	float D[1][1] = {0.0};
	float x_0[2][1] = {{0.0}, {0.0}};
	float u[1][100];
	float y[1][100];
	
	HMI_snput_matrix(8, 2, 2, &A);
	HMI_snput_matrix(8, 2, 1, &B);
	HMI_snput_matrix(8, 1, 2, &C);
	HMI_snput_matrix(8, 1, 1, &D);
	HMI_snput_matrix(8, 2, 1, &x_0);

	HMI_snget_matrix(8, 2, 2, &A);
	HMI_snget_matrix(8, 2, 1, &B);
	HMI_snget_matrix(8, 1, 2, &C);
	HMI_snget_matrix(8, 1, 1, &D);
	HMI_snget_matrix(8, 2, 1, &x_0);
	HMI_snget_matrix(8, 1, 100, &u);
	
	float x1[2][1];
	float x2[2][1];
	x1[0][0]=x_0[0][0];
	x2[0][0]=x_0[0][0];
	x1[1][0]=x_0[1][0];
	x2[1][0]=x_0[1][0];
	
	for(int k = 0; k < 100; k++){
		x2[0][0]=A[0][0]*x1[0][0]+A[0][1]*x1[1][0]+B[0][0]*u[0][k];
		x2[1][0]=A[1][0]*x1[0][0]+A[1][1]*x1[1][0]+B[1][0]*u[0][k];
		y[0][k]=C[0][0]*x1[0][0]+C[0][1]*x1[1][0]+D[0][0]*u[0][k];
		x1[0][0]=x2[0][0];
		x1[1][0]=x2[1][0];
	}
	HMI_snput_matrix(8, 1, 100, &y);

}

#### MATLAB

clear       %清除記憶體內儲存變數
clc         %清除command window訊息
close all   %清除已開啟的作圖視窗

load datau;
port = remo_open(9);

A = remo_snget_matrix(port);
B = remo_snget_matrix(port);
C = remo_snget_matrix(port);
D = remo_snget_matrix(port);
x_0 = remo_snget_matrix(port);

A = [1.35 0.55;
    -0.45 0.35];
B = [0.5;
     0.5];
C = [3 1];
D = [1];
x_0 = [0;
       0];
   
A = cast(A,'single');
B = cast(B,'single');
C = cast(C,'single');
D = cast(D,'single');
x_0 = cast(x_0,'single');
   
remo_snput_matrix(port, A);
remo_snput_matrix(port, B);
remo_snput_matrix(port, C);
remo_snput_matrix(port, D);
remo_snput_matrix(port, x_0);
remo_snput_matrix(port, u);

y = remo_snget_matrix(port);
remo_close(port);

k = 0.01:0.01:1;
plot(k, u, 'B-o', k, y, 'R-x');

xlabel('t');  
ylabel('u(k) and y(k)')  
title('u(k) and y(k)');  
legend('u(k)=sin(wk)','y(k)');

## 流程圖
### DataAgent輔助開發狀態空間計算系統

### Matlab Remo put get 巨集輔助開發狀態空間計算系統

## 實驗數據
### DataAgent輔助開發狀態空間計算系統
#### 第一組數據

f = 5HZ
w = 2 * pi * f
u=sin(wt)


#### 第二組數據
f = 10HZ
w = 2 * pi * f
u=sin(wt)


### Matlab Remo put get 巨集輔助開發狀態空間計算系統
#### 第一組數據
f = 5HZ
w = 2 * pi * f
u=sin(wt)


#### 第二組數據
f = 10HZ
w = 2 * pi * f
u=sin(wt)


## 驗收題目
:::info
> [name=cy023]

- 題目:
1. 請先了解 [Moving Average](https://en.wikipedia.org/wiki/Moving_average) 演算法。
2. 由 Matlab 傳送 500 筆 `uint8_t` 型態變數至 M128 端。
3. 在 M128 端實作 **Exponential Moving Average** 演算法，將上述資料做濾波後送回 Matlab。
4. Matlab 將原始數據與 M128 計算結果疊圖輸出，比較濾波效果。

- 要求:
1. 數據不限於亂數產生、蒐集股價資訊 ... 等
2. 至少需 3 種不同 **Exponential Moving Average** 參數，並討論其差異。
3. 需與以額外方式驗證演算法正確性，不限於使用 Matlab built-in function、網路上現成計算機...等。

- 延伸閱讀:
[An Easy-to-Use Digital Filter](https://blog.stratifylabs.co/device/2013-10-04-An-Easy-to-Use-Digital-Filter/)
[股票投資新手攻略：解構股票用語　技術分析篇](https://www.moneyhero.com.hk/blog/zh/%E8%82%A1%E7%A5%A8%E6%8A%95%E8%B3%87%E6%96%B0%E6%89%8B%E6%94%BB%E7%95%A5-%E8%A7%A3%E6%A7%8B%E8%82%A1%E7%A5%A8%E7%94%A8%E8%AA%9E-%E6%8A%80%E8%A1%93%E5%88%86%E6%9E%90%E7%AF%87)
:::

###
股票選擇:1101台泥
日期:2021/12/1~2022/9/23 共200天
#### 第一組數據
EMA = 5
前5日EMA設為0，從第11日開始畫

藍線為收盤價，棕線為EMA
## Matlab Plot

## Tradingview 驗證

程式碼

#include "c4mlib.h"
# define K 200
# define n 5


int main()
{
	C4M_DEVICE_set();
	float x[1][K];
	float y[1][K];
	float temp[1][K];
	
	float alpha;//平滑constant
	
	alpha = 2.0 / (n + 1);
	//printf("alpha = %f\n", alpha);
	
	printf("Please give matrix x\n");
	HMI_snget_matrix(8, 1, K, &x);
	
	
	for(int i = 0; i < K; i++){
		if(i < n){
			temp[0][i] = 0;
		}
		else{
			//printf("%f\n",x[0][i]);
			temp[0][i] = temp[0][i - 1]*(1 - alpha) + x[0][i]* alpha;
		}
		if(i < 2*n){
			y[0][i] = 0;
		}
		else{
			y[0][i] = temp[0][i];
		}
		//printf("%f\n",y[0][i]);
	}
	HMI_snput_matrix(8, 1, K, &y);
	printf("end\n");

}

#### 第二組數據
EMA = 20
前20日EMA設為0，從第21日開始畫

藍線為收盤價，棕線為EMA
## Matlab Plot

## Tradingview 驗證

程式碼

#include "c4mlib.h"
# define K 200
# define n 20


int main()
{
	C4M_DEVICE_set();
	float x[1][K];
	float y[1][K];
	float temp[1][K];
	
	float alpha;//平滑constant
	
	alpha = 2.0 / (n + 1);
	//printf("alpha = %f\n", alpha);
	
	printf("Please give matrix x\n");
	HMI_snget_matrix(8, 1, K, &x);
	
	
	for(int i = 0; i < K; i++){
		if(i < n){
			temp[0][i] = 0;
		}
		else{
			//printf("%f\n",x[0][i]);
			temp[0][i] = temp[0][i - 1]*(1 - alpha) + x[0][i]* alpha;
		}
		if(i < 2*n){
			y[0][i] = 0;
		}
		else{
			y[0][i] = temp[0][i];
		}
		//printf("%f\n",y[0][i]);
	}
	HMI_snput_matrix(8, 1, K, &y);
	printf("end\n");

}

#### 第三組數據
EMA = 50
前50日EMA設為0，從第51日開始畫

藍線為收盤價，棕線為EMA
## Matlab Plot

## Tradingview 驗證

程式碼

#include "c4mlib.h"
# define K 200
# define n 50


int main()
{
	C4M_DEVICE_set();
	float x[1][K];
	float y[1][K];
	float temp[1][K];
	
	float alpha;//平滑constant
	
	alpha = 2.0 / (n + 1);
	//printf("alpha = %f\n", alpha);
	
	printf("Please give matrix x\n");
	HMI_snget_matrix(8, 1, K, &x);
	
	
	for(int i = 0; i < K; i++){
		if(i < n){
			temp[0][i] = 0;
		}
		else{
			//printf("%f\n",x[0][i]);
			temp[0][i] = temp[0][i - 1]*(1 - alpha) + x[0][i]* alpha;
		}
		if(i < 2*n){
			y[0][i] = 0;
		}
		else{
			y[0][i] = temp[0][i];
		}
		//printf("%f\n",y[0][i]);
	}
	HMI_snput_matrix(8, 1, K, &y);
	printf("end\n");

}

#### MATLAB

clear   %清除記憶體內儲存變數
clc       %清除command window訊息
close all   %清除已開啟的作圖視窗

load datax;
load datay.mat;


k = 1:1:200;

plot(k , x, 'B-', k, data0, '.');



xlabel('Day');  
ylabel('收盤價 , EMA5');  
title('台泥1101');  
legend('收盤價','EMA5');

# 討論
EMA是取加權平均，算出來的值會比SMA更靈活。所選天數越小越貼近股價。由於N天的EMA前N天沒辦法算平均，因此我把前N天設為0，而由於前N天用0去算，EMA會有誤差，所以取2N天才畫EMA線。