adc fft代码求助

wstudouoooooo

filename = '128point.csv';
data0 = readmatrix(filename);

data1 = data0(:, 2);
ns = 128; % 需要的样本数量
np = 1; % 相关采样的参数
if length(data1) >= ns
    data1 = data1(end - ns + 1:end); % 提取最后128个数据
else
    error('数据长度不足128个样本，请检查数据文件。');
end
% 对时间序列数据减去直流分量
data = data1 - mean(data1);
fs = 1e6/3; % 采样频率
%% 5. 供 FFT 使用的关键量
point = ns;                     % 样点数
fin   = np * fs / point;        % 信号中心频 (Hz)
fprintf('Fin ≈ %.2f Hz\n', fin);
% ------------------  修正后的 FFT 段 ------------------
%% FFT
Fs = fs;
Fin = np * Fs / ns;
OSR = 1;
band = Fs / (2 * OSR);
data = data .* hann(length(data), 'periodic')'; % 加汉宁窗
y = fft(data);
point = length(data);
pyy = y .* conj(y) / point^2 * 16; % 功率谱计算
% 关键修正1：将fbin转为整数（四舍五入），确保索引合法
fbin = find(pyy == max(pyy), 1);  % 直接定位功率谱中最大值对应的索引（最准确）
% 额外校验：确保fbin在有效范围内（1到数组长度）
fbin = max(1, min(fbin, length(pyy)));
% 计算信号功率和噪声功率
pf = sum(pyy(max(1, fbin-1) : min(length(pyy), fbin+1))); % 取附近3个点（避免越界）
ND = sum(pyy([3 : max(1, fbin-2), min(length(pyy), fbin+2) : end])); % 排除信号主峰范围
band_bin = ceil(band / Fs * length(data));
ND_inband = sum(pyy([3 : max(1, fbin-2), min(band_bin, fbin+2) : band_bin]));
% 计算性能指标
SNDR = 10 * log10(pf / ND);
SNDR_inband = 10 * log10(pf / ND_inband);
hd_top = floor(band / Fin);
hd = zeros(1, hd_top);
for j = 2 : hd_top
    harmonic_bin = round(j * fbin); % 谐波频点转为整数
    harmonic_bin = max(1, min(harmonic_bin, length(pyy))); % 边界保护
    hd(j) = sum(pyy(max(1, harmonic_bin-1) : min(length(pyy), harmonic_bin+1)));
end
thd = -10 * log10(pf / sum(hd));
snr = -10 * log10((ND_inband - sum(hd)) / pf);
% 绘图部分
figure(2)
fplot = (1 : point/2) / length(data) * Fs; % 频率轴刻度
% 标注参数文本框
str = {
    sprintf('Nsample = %d', ns),
    sprintf('M = %d', np),
    sprintf('Fsignal = %.2f Hz', fin),
    sprintf('BW = %.2e Hz', band),
    sprintf('RBW = %.2e Hz', fplot(2)-fplot(1)),
    sprintf('SNDR = %.2f dB', SNDR_inband),
    sprintf('SNR = %.2f dB', snr),
    sprintf('THD = %.2f dB', thd)
};
dim = [0.6 0.1 0.3 0.3];
annotation('textbox', dim, 'String', str, 'FitBoxToText', 'on');
% 功率谱转换为dB
pyydb = 10 * log10([pyy(2 : point/2), pyy(point/2)/2]);
% 绘制整体频谱
semilogx(fplot, pyydb, 'b', 'LineWidth', 2);
hold on;
% 关键修正2：信号主峰绘制（添加边界判断，避免索引越界）
start_idx = max(1, fbin - 2);
end_idx = min(length(pyydb), fbin + 2);
semilogx(fplot(start_idx : end_idx), pyydb(start_idx : end_idx), 'k', 'LineWidth', 2);
% 绘制谐波（同样添加边界保护）
for j = 2 : hd_top
    harmonic_bin = round(j * fbin);
    h_start = max(1, harmonic_bin - 2);
    h_end = min(length(pyydb), harmonic_bin + 2);
    semilogx(fplot(h_start : h_end), pyydb(h_start : h_end), 'r', 'LineWidth', 2);
end
% 绘制带宽边界
plot([band, band], [-300, 100], 'g', 'LineWidth', 2);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('频率谱');
grid on;   
disp('DONE');

好心人帮忙看看这个adc fft代码行不行呀

                               
登录/注册后可看大图

wstudouoooooo

为什么我跑出来的结果是这样的