% 清除工作区变量、关闭所有图形窗口、清空命令行
clear; close all; clc;
%% 参数设置
Fs = 3200; % 采样频率 3200 Hz
N = 10000; % 数据点总数 10000
t = (0:N-1)/Fs; % 时间向量（用于生成示例信号）
%% 生成示例信号（用户应替换为自己的数据）
% 合成信号：包含10Hz和100Hz的正弦波 + 高斯噪声
f1 = 10; % 第一个频率分量 10Hz
f2 = 100; % 第二个频率分量 100Hz
data = 0.5*sin(2*pi*f1*t) + sin(2*pi*f2*t) + 0.2*randn(size(t));
%% 预处理（可选）
data = data - mean(data); % 去除直流分量（中心化处理）
%% 快速傅里叶变换（FFT）
Y = fft(data); % 执行FFT，结果Y是复数数组
P2 = abs(Y/N); % 计算双边频谱幅度，并归一化
%% 转换为单边频谱
P1 = P2(1:N/2+1); % 取前半部分频谱（0Hz ~ Fs/2）
P1(2:end-1) = 2*P1(2:end-1); % 幅度修正（除直流和Nyquist频率外均×2）
%% 构建频率轴
f = Fs*(0:(N/2))/N; % 创建频率向量（0Hz ~ Fs/2）
%% 绘制时域信号
figure;
subplot(2,1,1)
plot(t, data)
title('时域信号')
xlabel('时间 (s)')
ylabel('幅度')
xlim([0 0.3]) % 显示前0.3秒数据
grid on;
%% 绘制频域信号
subplot(2,1,2)
plot(f, P1)
title('单边幅度频谱')
xlabel('频率 (Hz)')
ylabel('幅度')
grid on;
xlim([0 Fs/2]) % 显示完整频率范围
%% 峰值标注示例（可选）
[~, locs] = findpeaks(P1, 'SortStr','descend', 'NPeaks',2); % 找前两个峰值
hold on;
plot(f(locs), P1(locs), 'ro'); % 标出峰值点
text(f(locs(1))+5, P1(locs(1)), [num2str(f(locs(1)), '%.1f'), ' Hz']);
text(f(locs(2))+5, P1(locs(2)), [num2str(f(locs(2)), '%.1f'), ' Hz']);
hold off;
