# Lecture 7 – 一維數據擬合 (1D Data Fitting)；Exp7 聲波拍頻 (Sound Wave Beating)

## 課程目標
- 學會利用 MATLAB 對實驗數據進行擬合 (data fitting)
- 練習自訂擬合模型，並理解數據與物理模型的關聯
- 觀察與分析聲波拍頻現象，並找出其頻率成分

---

## 一、數據擬合 (Data Fitting)

### 1. 基本概念
- 實驗數據因為有限取樣率，通常是 **離散點**。
- 擬合 (fitting)：找出一條「最佳函數曲線」來描述數據點。
- 可用統計方法中的 **迴歸分析 (regression analysis)** 完成。

> 延伸閱讀：[Curve fitting](https://en.wikipedia.org/wiki/Curve_fitting)、[Regression analysis](https://en.wikipedia.org/wiki/Regression_analysis)

---

### 2. 範例：自由落體求 g
已知自由落體公式：
$$
x = \tfrac{1}{2} g t^2
$$

若量測高度 \(x\) 與時間 \(t\) 數據如下：

| 高度 x (m) | 1.38 | 1.28 | 1.18 | 1.08 | 0.98 | 0.88 | 0.78 | 0.68 | 0.58 |
|------------|------|------|------|------|------|------|------|------|------|
| 時間 t (s) | 0.53 | 0.51 | 0.49 | 0.47 | 0.45 | 0.43 | 0.40 | 0.37 | 0.34 |

在 MATLAB 中可用 **Basic Fitting Tool** 或 **程式碼**進行擬合。

---

#### Method 1：Basic Fitting Toolbox
- `Figure → Tools → Basic Fitting` → 選擇 `linear` 擬合
- 顯示方程式，斜率 \( s = 0.2045 \)
- 由公式 \( g = 2/s \)，得 \( g = 9.78 \, m/s^2 \)

---

#### Method 2：程式碼擬合

x = [1.38, 1.28, 1.18, 1.08, 0.98, 0.88, 0.78, 0.68, 0.58];
t = [0.53, 0.51, 0.49, 0.47, 0.45, 0.43, 0.40, 0.37, 0.34];

t_s = t.^2;

% 內建 linear 擬合
f_poly1 = fit(x', t_s', 'poly1');

% 自訂強迫通過 (0,0)
ft1 = fittype({'x'});
f_ft1 = fit(x', t_s', ft1);

% 繪圖
X = linspace(0, 1.1*max(x), 200);
plot(x, t_s, 'ro'); hold on;
plot(X, f_poly1(X), 'b--', X, f_ft1(X), 'k')
xlabel('Height (m)')
ylabel('t^2 (s^2)')
legend({'Raw Data', 'poly1', 'Through (0,0)'})

👉 兩種方法差異不大，但強迫通過原點更符合物理模型。
### 3. polyfit 與 polyval

x = 1:20;
y = 3*x.^2 + 0.25*x - 19;

p = polyfit(x,y,2);  % 2階多項式擬合
f = polyval(p,x);

plot(x,y,'ro', x,f,'b-')


- polyfit：回傳係數陣列 [$p_2$,$p_1$,$p_0$]
- polyval：計算擬合函數在指定 x 的值

---
### 4. 進階範例
**Quadratic fit – 上拋運動**

y0 = 50; vy0 = 5; g0 = 9.8;
Nexp = 15;
Tfly = vy0/g0 + sqrt(vy0^2+2*g0*y0)/g0;
t_exp = linspace(0, Tfly, Nexp);

y_ideal = y0 + vy0.*t_exp - 0.5*g0.*(t_exp.^2);
y_exp = y_ideal + 2*(rand(1,Nexp)-0.5); % 加入雜訊

f_poly2 = fit(t_exp', y_exp', 'poly2');
plot(f_poly2, t_exp, y_exp, 'bo')


---
**Sinusoidal fit – 簡諧振盪**

A = 2; T = 3; w = 2*pi/T; phi = 0.32*pi;
Nexp = 60;
t_exp = linspace(0,2*T,Nexp);
x_ideal = A*cos(w*t_exp + phi);
x_exp = x_ideal + 0.3*A*(rand(1,Nexp)-0.5);

f_sin = fit(t_exp', x_exp', 'sin1');
plot(f_sin, t_exp, x_exp, 'ro')


👉 可由擬合結果得到 振幅 A、角頻率 ω、相位 φ。

**自定義 fittype – 阻尼振子**

A=3; beta=0.3; T=10; w0=2*pi/T; phi=0.2*pi;
Nexp=30;
t_exp = linspace(0,5*T,Nexp);
x_ideal = A*exp(-beta*w0*t_exp).*sin(sqrt(1-beta^2)*w0*t_exp+phi);
x_exp = x_ideal + 0.02*A*(rand(1,Nexp)-0.5);

s = fitoptions('Method','NonlinearLeastSquares','Startpoint',[]);
damp_osc = fittype('a1*exp(-a2*a3*x).*sin(sqrt(1-a2^2)*a3*x+a4)',...
    'independent','x','dependent','Y','options',s);

f_dampOsc = fit(t_exp', x_exp', damp_osc);
plot(f_dampOsc, t_exp, x_exp, 'bo')

---
## 二、Exp 7：聲波拍頻 (Sound Wave Beating)

### 實驗目的
- 觀察 **拍頻 (beat)** 現象
- 使用 **Phyphox app** 錄製聲音訊號並輸出數據
- 在 MATLAB 中分析聲音組成頻率

---

### 實驗原理
當兩個頻率相近的聲波 \( f_1, f_2 \) 疊加時，會形成 **拍頻** 現象。

數學式：
$$
y = \sin(2\pi f_1 t) + \sin(2\pi f_2 t)
$$

可化簡為：
$$
y = 2 \cos\!\left(2\pi \tfrac{f_1-f_2}{2} t\right) \cdot \sin\!\left(2\pi \tfrac{f_1+f_2}{2} t\right)
$$

- 外部包絡 (低頻) → 拍頻頻率
$$
f_{\text{beat}} = |f_1 - f_2|
$$

- 內部振盪 (高頻) → 平均頻率
$$
f_{\text{avg}} = \frac{f_1 + f_2}{2}
$$

👉 聽起來像「忽大忽小」的聲音強弱變化。

---

### 實驗步驟
1. 打開 [Tone Generator](https://onlinetonegenerator.com/multiple-tone-generator.html)，設定兩個接近的頻率 (例如 440 Hz 與 445 Hz)。
2. 使用 **Phyphox app → 聲學工具 → 時變相對音強** 來錄製聲音。
3. 將錄音數據匯出為 `.csv` 或 Excel 檔案。
4. 在 MATLAB 匯入數據並繪圖
5. 使用 fit 或傅立葉轉換 (fft) 分析決定出合成訊號的頻率、振幅以及相位差資訊，如圖所示。
6. 驗證觀測到的拍頻是否符合理論值 |$f_1-f_2$|
---
### 延伸思考
- 為什麼樂器調音時會聽到拍頻？
- 當兩頻率差變大時，拍頻的感受會如何改變？
- 若兩聲波振幅不同，拍頻是否仍清晰？
---
### 參考資料
- https://en.wikipedia.org/wiki/Beat_(acoustics)