sfft.m

function [f, S, P] = sfft(x, Fs, draw)
% 计算单边幅值谱
% 输入: x 信号向量, Fs 采样频率
% 输出: f 频率向量(Hz), S 单边幅值谱

if nargin < 3
    draw = false;
end

N   = length(x);
X   = fft(x);

halfN = floor(N/2);
idx = 1:(halfN + 1);

P = angle(X(idx));
S = abs(X(idx)) / N;

% 单边谱幅值修正：DC 不翻倍；偶数长度时 Nyquist 点也不翻倍
if N > 1
    if mod(N, 2) == 0
        if numel(S) > 2
            S(2:end-1) = 2 * S(2:end-1);
        end
    else
        S(2:end) = 2 * S(2:end);
    end
end

f = (0:halfN) * (Fs / N);

if draw
    figure;
    subplot(2,1,1);
    eval([draw '(f, S);']);
    xlabel('Frequency (Hz)'); xlim([0 Fs/2]);
    ylabel('Amplitude');
    title('Single-Sided Amplitude Spectrum');
    subplot(2,1,2);
    eval([draw '(f, rad2deg(P));']);
    xlabel('Frequency (Hz)'); xlim([0 Fs/2]);
    ylabel('Phase (degrees)');
end