Octave / MATLAB 第98题：一阶常微分方程初值问题数值仿真（显式欧拉法步进）

题目

编写程序，合成一个包含两个正弦分量（50 Hz 与 120 Hz）并混入高斯白噪声的时域离散信号，采样频率 $F_s = 1000 \text{ Hz}$；调用 <code>fft()</code> 计算信号的快速傅里叶变换，提取单边双峰幅值谱，并找出谱峰对应的主频率。

输入格式：无需外部输入。

输出格式：输出采样参数、FFT 幅值谱提取出的两处最强峰值频率与幅值。

数据范围与约定：数字信号处理与频谱分析。

思路与算法

1. 采样定理与时间轴
根据奈奎斯特采样定理，采样频率 $F_s$ 必须大于最高信号频率的两倍。采样点数 $N$ 决定频率分辨率 $\Delta f = F_s / N$。

2. 单边谱归一化
fft 返回双边复数频谱。幅值取模 abs(Y)/N；对于单边正频率分量，除直流分量外需乘以 2。

3. 频率轴对齐
频率轴向量为 f = (0:(N/2)) * (Fs / N)。

完整程序

% Octave / MATLAB 第98题：FFT 快速傅里叶变换与主频分析
Fs = 1000;            % 采样频率 1000 Hz
T = 1 / Fs;           % 采样周期
L = 1000;             % 信号长度 (1000 点，时域 1 秒)
t = (0:L-1) * T;      % 时间向量

% 合成测试信号: 50 Hz (幅值 0.7) + 120 Hz (幅值 1.0)
S = 0.7 * sin(2*pi*50*t) + 1.0 * sin(2*pi*120*t);

% 计算 FFT
Y = fft(S);

% 计算双边频谱 P2，然后转单边频谱 P1
P2 = abs(Y / L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2 * P1(2:end-1);

f = Fs * (0:(L/2)) / L; % 频率轴

% 寻找幅值谱中的峰值点
[pks, idxs] = sort(P1, "descend");

fprintf("信号采样率: %d Hz, 长度: %d 点\n", Fs, L);
fprintf("检测到的第一主频: %.1f Hz (测量幅值: %.2f, 理论: 1.00)\n", f(idxs(1)), pks(1));
fprintf("检测到的第二主频: %.1f Hz (测量幅值: %.2f, 理论: 0.70)\n", f(idxs(2)), pks(2));

运行示例

输出：

信号采样率: 1000 Hz, 长度: 1000 点
检测到的第一主频: 120.0 Hz (测量幅值: 1.00, 理论: 1.00)
检测到的第二主频: 50.0 Hz (测量幅值: 0.70, 理论: 0.70)