之前找到的別人的代碼，感覺多余無用的地方太多了，重寫了一版，附了詳細(xì)的注釋解析

注意不是用于直接解析數(shù)字碼的，而是解析數(shù)字碼經(jīng)過DAC后轉(zhuǎn)成的階梯波的。

如圖，藍(lán)色部分的Dout

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%? ? ? ? ? ? ? ? ? ? FFT calculate SNDR,ENOB? ? ? ? ? ? ? ? ? ? ? ?%?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 以12bit為例

Dout=reshape(Dout,'',1)';%轉(zhuǎn)置成行向量

Dout=Dout-(max(Dout(1:end))+min(Dout(1:end)))/2;%去除直+交中的直流能量。

Dout=Dout(1:FFTN);%取出前“采樣點(diǎn)數(shù)個(gè)”信號，如開始未穩(wěn)定，可改為（x:FFTN+x）

Window=hann(FFTN)';%漢寧窗口，注意轉(zhuǎn)置

DoutW=Dout.*Window;%1*4096

DoutFFT=fft(DoutW,FFTN);%1*4096 complex

PSD=(abs(DoutFFT)).*(abs(DoutFFT));%功率譜密度

span = 3;

DcPower = sum(PSD(1:span))

[maxPower,FinBin] = max(PSD);%maxPower沒用到，只需要FinBin

FinBin = FinBin - 1;%第一個(gè)Bin是直流，信號Bin在1+31，所以要減去1，得到實(shí)際的Bin

if FinBin > span % 信號不在直流點(diǎn)內(nèi)

? ? SignalPower = sum(PSD(FinBin - span:FinBin + span));

else

? ? SignalPower = 0;

end

SignalHarmonyFrequency = zeros(1,fix(FFTN/2/FinBin));%存儲信號與諧波的能量

SignalHarmonyPower = zeros(1,fix(FFTN/2/FinBin));%存儲信號與諧波的頻率，未用到

SignalHarmonyNumbers = fix(FFTN/2/FinBin);%信號帶寬內(nèi)的諧波數(shù)量

SignalHarmonyPower(1)=SignalPower;% 信號能量

SignalHarmonyFrequency(1)=Fs*(FinBin)/FFTN; %信號頻率

for i = 1:SignalHarmonyNumbers%遍歷每一個(gè)諧波（注意循環(huán)內(nèi)只計(jì)算諧波，不計(jì)算信號）

? ??

? ? HarmonyBin=FinBin*(i+1);%從第二個(gè)諧波點(diǎn)開始，計(jì)算下一個(gè)諧波位于的Bin??

? ??

? ? if (HarmonyBin <= FFTN/2)%如果Bin在BW's Bin內(nèi)，執(zhí)行下述操作

? ? ? ? SignalHarmonyPower(i+1)=sum(PSD(HarmonyBin-span:HarmonyBin+span));

? ? ? ? SignalHarmonyFrequency(i+1)=(HarmonyBin)*Fs/FFTN;

? ? end

? ??

end

NoisePower=sum(PSD(1:FFTN/2))-DcPower-SignalPower;

SNDR=10*log10(SignalPower/NoisePower)

ENOBloc=(SNDR-1.76)/6.02

HarmonyPower = sum(SignalHarmonyPower(1:SignalHarmonyNumbers))- SignalHarmonyPower(1);

THD=10*log10(HarmonyPower/SignalHarmonyPower(1))

SFDR=10*log10(SignalHarmonyPower(1)/max(SignalHarmonyPower(2:SignalHarmonyNumbers)))

? ??

figure;

hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%? ? ? 歸一化，僅畫圖時(shí)用到? ? ?%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PSD_dB = 10*log10(PSD/max(PSD));

PSD_dB = PSD_dB(1:FFTN/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

semilogx([0:FFTN/2-1] * Fs / FFTN, PSD_dB(1:FFTN/2), '-');

text_handle = text(Fin * 5, -20, sprintf('SNDR = %4.1f dB \nENOB = %2.2f bits ', SNDR, ENOBloc), ...

? ? 'EdgeColor', 'red', 'LineWidth', 3, 'BackgroundColor', [1 1 1], 'Margin', 10);

title('FREQ DOMAIN');

xlabel('ADC FFT SPECTRUM');

ylabel('AMPLITUDE(dB)');

set(gca, 'XScale', 'log');? % 將當(dāng)前坐標(biāo)軸改為對數(shù)軸，去除該行則線性x軸坐標(biāo)。

hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Python版

######### FFT calculate SNR & ENOB ###########

# 以12bit為例

Dout = Dout[:FFTN]

Dout = Dout - (np.max(Dout) + np.min(Dout)) / 2

numpt2 = int(np.fix(FFTN / 2))

Window = np.hanning(FFTN)

Doutw = Dout * Window

DoutFFT = np.fft.fft(Doutw, FFTN)

PSD = (np.abs(DoutFFT)) ** 2

FinBin = np.argmax(PSD[:numpt2])? # 獲取最大值索引

span = 3

DcPower = np.sum(PSD[0:span])

if FinBin > span:

? ? SignalPower = np.sum(PSD[FinBin - span:FinBin + span])

else:

? ? SignalPower = 0

SignalHarmonyNumbers = int(np.fix(FFTN/2/FinBin))

SignalHarmonyFrequency = np.zeros(SignalHarmonyNumbers)

SignalHarmonyPower = np.zeros(SignalHarmonyNumbers)

SignalHarmonyPower[0] = SignalPower? # 信號能量

SignalHarmonyFrequency[0] = Fs * (FinBin) / FFTN? # 信號頻率

for i in range(1, SignalHarmonyNumbers):? # 遍歷每一個(gè)諧波（注意循環(huán)內(nèi)只計(jì)算諧波，不計(jì)算信號）

? ? HarmonyBin = FinBin * (i + 1)? # 從第二個(gè)諧波點(diǎn)開始，計(jì)算下一個(gè)諧波位于的Bin

? ? if HarmonyBin <= FFTN / 2:? # 如果Bin在BW's Bin內(nèi)，執(zhí)行下述操作

? ? ? ? SignalHarmonyPower[i] = np.sum(PSD[HarmonyBin - span:HarmonyBin + span])

? ? ? ? SignalHarmonyFrequency[i] = (HarmonyBin) * Fs / FFTN

? ? ? ??

AllPower = np.sum(PSD[0:numpt2])

NoisePower = AllPower - DcPower - SignalPower

SNDR = 10 * np.log10(SignalPower / NoisePower)

ENOBloc = (SNDR - 1.76) / 6.02

HarmonyPower = np.sum(SignalHarmonyPower[1:SignalHarmonyNumbers]) - SignalHarmonyPower[0]

THD = 10 * np.log10(np.sqrt(np.sum(SignalHarmonyPower[1:]**2)) / SignalHarmonyPower[0])

print(THD)

SFDR = 10 * np.log10(SignalHarmonyPower[0] / np.max(SignalHarmonyPower[1:SignalHarmonyNumbers]))

print(SFDR)

plt.figure(1)

plt.plot(Dout[:FFTN])

plt.xlabel('Sample Point')

plt.ylabel('Dout')

plt.title('Dout Signal')

plt.grid()

maxPower = np.max(PSD)

PSD_dB = 10 * np.log10(PSD / maxPower)??

plt.figure(2)

plt.semilogx(np.arange(numpt2) * Fs / FFTN, PSD_dB[:numpt2], '-')

plt.text(Fs / 3, -20, f'SNDR = {SNDR:.1f} dB \nENOB = {ENOBloc:.2f} bits',

? ? ? ? ?bbox=dict(facecolor='red', edgecolor='red', linewidth=3, pad=10))

plt.title('Normalized Power Spectrum')

plt.xlabel('Frequency (Hz)')

plt.ylabel('Normalized Power (dB)')

plt.grid()

plt.show()