TI ADS7042 EVM ADC数据Python SNR计算结果异常排查

我从TI ADS7042 EVM套件获取了12位ADC输出数据，输出码值范围为0-4095，尝试计算SNR、THD+N、SINAD（单位为dB）等参数。数据文件为ADC-output2.txt。

我使用以下Python代码计算SNR和THD+N：

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import scipy.io

def preprocess_adc_data(adc_data, v_ref, adc_resolution):
    adc_data_np = np.array(adc_data)
    adc_data_np = adc_data_np - np.mean(adc_data_np)  

    # 7 term Blackman-Harris window
    window = signal.windows.blackmanharris(len(adc_data_np))
    adc_data_np = adc_data_np * window

    return (adc_data_np / (2 ** adc_resolution - 1)) * v_ref


def estimate_input_frequency_from_spectrum(adc_data_voltage, fs):
    # Spectrum Calculate
    f, Pxx = signal.periodogram(adc_data_voltage, fs)

    # Find the frequency corresponding to the maximum value in the power spectrum.
    max_index = np.argmax(Pxx)
    frequency = f[max_index]

    return frequency

def calculate_snr(adc_data_voltage, fs):
    f, Pxx = signal.periodogram(adc_data_voltage, fs)
    snr = 10 * np.log10(np.max(Pxx) / (np.sum(Pxx) - np.max(Pxx)))
    return snr

def plot_signal(adc_data_voltage, fs):
    t = np.arange(len(adc_data_voltage)) / fs
    plt.figure()
    plt.plot(t, adc_data_voltage)
    plt.xlabel('Time (s)')
    plt.ylabel('Voltage (V)')
    plt.title('Restored Signal Waveform')
    plt.show()

def plot_spectrum(adc_data_voltage, fs, f_input):
    nfft = int(fs)
    f, Pxx = signal.periodogram(adc_data_voltage, fs, nfft=nfft)
    Pxx_dBc = 10 * np.log10(Pxx / np.max(Pxx))
    # Maximum frequency
    f_limit = 1e5
    idx = np.where(f <= f_limit)
    f_limited = f[idx]
    Pxx_dBc_limited = Pxx_dBc[idx]

    fig, ax = plt.subplots()
    ax.plot(f_limited, Pxx_dBc_limited)
    ax.set_xlabel("Frequency (Hz)")
    ax.set_ylabel("Amplitude (dBC)")
    ax.set_title("Spectrum")

    # Limit X-axis and Y-axis
    ax.set_xlim([0, f_limit])
    ax.set_ylim([-160, 20])

    # Compute the fundamental frequency and harmonics.
    fundamental_freq = int(f_input)
    first_harmonic_freq = 2 * fundamental_freq
    second_harmonic_freq = 3 * fundamental_freq
    third_harmonic_freq = 4 * fundamental_freq

    # Locate the positions of the fundamental frequency and harmonics in the spectrum.
    fundamental_idx = np.where(f == fundamental_freq)
    first_harmonic_idx = np.where(f == first_harmonic_freq)
    second_harmonic_idx = np.where(f == second_harmonic_freq)
    third_harmonic_idx = np.where(f == third_harmonic_freq)

    # Extract the dBc values of the fundamental frequency and harmonics.
    fundamental_dBc = float(Pxx_dBc[fundamental_idx])
    first_harmonic_dBc = float(Pxx_dBc[first_harmonic_idx])
    second_harmonic_dBc = float(Pxx_dBc[second_harmonic_idx])
    third_harmonic_dBc = float(Pxx_dBc[third_harmonic_idx])

    print("H1：{:.0f} Hz, {:.2f} dBc".format(fundamental_freq, fundamental_dBc))
    print("H2：{:.0f} Hz, {:.2f} dBc".format(first_harmonic_freq, first_harmonic_dBc))
    print("H3：{:.0f} Hz, {:.2f} dBc".format(second_harmonic_freq, second_harmonic_dBc))
    print("H4：{:.0f} Hz, {:.2f} dBc".format(third_harmonic_freq, third_harmonic_dBc))

    plt.show()


def read_adc_data_from_file(file_path):
    with open(file_path, 'r') as file:
        lines = file.readlines()
    adc_data = [int(line.strip()) for line in lines]
    return adc_data


def calculate_thdn(adc_data_voltage, fs, fundamental_frequency, num_harmonics=7):
    f, Pxx = signal.periodogram(adc_data_voltage, fs)

    fundamental_power = None
    harmonics_power = 0

    for harmonic in range(1, num_harmonics + 1):
        harmonic_frequency = fundamental_frequency * harmonic
        index = np.argmin(np.abs(f - harmonic_frequency))
        power = Pxx[index]

        if harmonic == 1:
            fundamental_power = power
        else:
            harmonics_power += power

    thdn = np.sqrt(harmonics_power / fundamental_power)
    thdn_db = 20 * np.log10(thdn)
    return thdn_db


file_path = './ADC-output2.txt'  # File
adc_data = read_adc_data_from_file(file_path)


v_ref = 3.3
adc_resolution = 12
adc_data_voltage = preprocess_adc_data(adc_data, v_ref, adc_resolution)
sample_rate = 1e6  # 1 MHz

f_input = estimate_input_frequency_from_spectrum(adc_data_voltage, sample_rate)
print("Estimated input frequency: {:.2f} Hz".format(f_input))

# Plot the reconstructed signal waveform and the spectral analysis waveform.
plot_signal(adc_data_voltage, sample_rate)
plot_spectrum(adc_data_voltage, sample_rate, f_input)
# SNR
snr = calculate_snr(adc_data_voltage, sample_rate)
print("SNR: {:.2f} dB".format(snr))
# THD
thdn = calculate_thdn(adc_data_voltage, sample_rate, f_input)
print("THD+N: {:.2f} dB".format(thdn))

TI EVM套件的官方软件计算结果与我的代码结果存在显著差异，我的代码计算结果如下：

• H1：10009 Hz, -0.00 dBc
• H2：20018 Hz, -80.43 dBc
• H3：30027 Hz, -86.12 dBc
• H4：40036 Hz, -91.65 dBc
• SNR: -0.03 dB
• THD+N: -77.97 dB

Python计算的频谱图显示上述谐波分布。现恳请协助排查代码中导致SNR计算结果异常的问题。

内容的提问来源于stack exchange，提问作者Karsten