Matlab 频谱分析 (Spectral Analysis)

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]\cdot 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\cdot \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\cdot N}$$

这里除以 $f_s\cdot N$ 是为了得到功率谱密度（单位：功率/Hz）。

Step 2: 单边谱转换
由于实信号的FFT结果具有共轭对称性：
$$X[k]=X^{\ast }[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\cdot 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]\cdot 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]|²/(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