---

1. 信號預處理 - 去直流分量

data = data - mean(data);

數學原理：
去除信號的直流分量（DC offset），即：
$$x_{centered}[n]=x[n]?\mu$$

其中：

• $\mu =\frac{1}{N}{\sum }_{n=0}^{N?1}x[n]$ 是信號均值
• $N$ 是信號長度

作用： 去除0Hz分量，避免在頻譜圖中出現很大的直流峰值，影響其他頻率成分的觀察。

2. 快速傅里葉變換（FFT）

Y = fft(data, nfft);
f = (0:nfft-1) * fs / nfft;

數學原理：
離散傅里葉變換（DFT）公式：
$$X[k]=\sum _{n=0}^{N?1}x[n]?e^{?j2\pi kn/N}$$

其中：

• $k=0,1,...,N?1$ 是頻率索引
• $N=nfft=2^{12}=4096$ 是FFT點數

頻率分辨率：
$$\Delta f=\frac{f_s}{nfft}$$

頻率向量計算：
$$f[k]=k?\frac{f_s}{nfft},k=0,1,...,nfft?1$$

3. 功率譜密度（PSD）計算

PSD = abs(Y).^2 / (fs * nfft);
PSD = PSD(1:nfft/2+1);
f = f(1:nfft/2+1);
PSD(2:end-1) = 2 * PSD(2:end-1);

數學原理：

Step 1: 計算功率譜
$$PSD[k]=\frac{|X[k]{|}^2}{f_s?N}$$

這里除以 $f_s?N$ 是為了得到功率譜密度（單位：功率/Hz）。

Step 2: 單邊譜轉換
由于實信號的FFT結果具有共軛對稱性：
$$X[k]=X^{?}[N?k]$$

所以只需要保留前一半頻率（0到$f_s/2$），但要將功率翻倍（除了0Hz和Nyquist頻率）：

PSD(2:end-1) = 2 * PSD(2:end-1);

數學表達：
$$PS{D}_{single}[k]=\{\begin{array}{ll}PSD[k],& k=0\text{?或?}k=N/2\\ 2?PSD[k],& 1\le k<N/2\end{array}$$

4. 主頻率檢測

[~, peak_idx] = max(PSD(2:end));
main_freq = f(peak_idx + 1);

數學原理：
找到功率譜密度的最大值對應的頻率：
$$f_{main}=\arg \underset{f>0}{\max }PSD(f)$$

注意：從索引2開始是為了排除直流分量（0Hz）。

5. 譜質心計算

centroid = sum(f_col .* PSD_col) / sum(PSD_col);

數學原理：
譜質心（Spectral Centroid）是頻譜的"重心"：
$$f_{centroid}=\frac{{\sum }_{k=0}^{N/2}f[k]?PSD[k]}{{\sum }_{k=0}^{N/2}PSD[k]}$$

這是一個加權平均，權重是各頻率的功率。

物理意義：

• 譜質心反映了信號能量在頻域的分布中心
• 值越大，說明高頻成分越多
• 常用于音頻信號分析中表征音色"明亮度"

6. 對數譜顯示

semilogx(f, 10*log10(PSD));

數學原理：
將功率譜密度轉換為分貝（dB）單位：
$$PS{D}_{dB}=10{\log }_{10}(PSD)$$

使用對數坐標的原因：

1. 人耳對頻率的感知是對數的
2. 可以同時顯示大范圍的幅值差異
3. 更容易觀察低幅值的頻率成分

完整的信號處理流程

原始信號 x[n]↓ (去直流)
中心化信號 x_c[n] = x[n] - μ↓ (FFT)
頻域信號 X[k]↓ (功率計算)
雙邊功率譜 |X[k]|2/(fs·N)↓ (單邊譜轉換)
單邊功率譜密度 PSD[k]↓ (特征提取)
主頻率 + 譜質心

實際應用示例

假設采樣率 fs = 1000 Hz，信號包含50Hz和120Hz兩個正弦波：

t = 0:1/fs:1;
x = sin(2*pi*50*t) + 0.5*sin(2*pi*120*t);

經過上述處理后：

• 主頻率將是 50 Hz（因為其幅度更大）
• 譜質心將在 50-120 Hz 之間，偏向50Hz