Matlab 频谱分析 (Spectral Analysis)

研究笔记更新于 2026-08-16编辑本文

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}, \quad 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^*[N-k]$

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

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

数学表达：

PSD_{single}[k] = \begin{cases} PSD[k], & k = 0 \text{ 或 } k = N/2 \\ 2 \cdot PSD[k], & 1 \leq k < N/2 \end{cases}

4. 主频率检测

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

数学原理： 找到功率谱密度的最大值对应的频率： $f_{main} = \arg\max_{f > 0} 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）单位： $PSD_{dB} = 10 \log_{10}(PSD)$

使用对数坐标的原因：

人耳对频率的感知是对数的
可以同时显示大范围的幅值差异
更容易观察低幅值的频率成分

完整的信号处理流程

原始信号 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

来源：CSDN 原文，首次发布于 2025-07-06。