用于ADC信号测试的PSD分析代码
之前找到的别人的代码，感觉多余无用的地方太多了，重写了一版，附了详细的注释解析
注意不是用于直接解析数字码的，而是解析数字码经过DAC后转成的阶梯波的。
如图，蓝色部分的Dout
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT calculate SNDR,ENOB %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 以12bit为例
Dout=reshape(Dout,'',1)';%转置成行向量
Dout=Dout-(max(Dout(1:end))+min(Dout(1:end)))/2;%去除直+交中的直流能量。
Dout=Dout(1:FFTN);%取出前“采样点数个”信号，如开始未稳定，可改为（x:FFTN+x）
Window=hann(FFTN)';%汉宁窗口，注意转置
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;%第一个Bin是直流，信号Bin在1+31，所以要减去1，得到实际的Bin
if FinBin > span % 信号不在直流点内
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);%信号带宽内的谐波数量
SignalHarmonyPower(1)=SignalPower;% 信号能量
SignalHarmonyFrequency(1)=Fs*(FinBin)/FFTN; %信号频率
for i = 1:SignalHarmonyNumbers%遍历每一个谐波（注意循环内只计算谐波，不计算信号）
HarmonyBin=FinBin*(i+1);%从第二个谐波点开始，计算下一个谐波位于的Bin
if (HarmonyBin <= FFTN/2)%如果Bin在BW's Bin内，执行下述操作
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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 归一化，仅画图时用到 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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'); % 将当前坐标轴改为对数轴，去除该行则线性x轴坐标。
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): # 遍历每一个谐波（注意循环内只计算谐波，不计算信号）
HarmonyBin = FinBin * (i + 1) # 从第二个谐波点开始，计算下一个谐波位于的Bin
if HarmonyBin <= FFTN / 2: # 如果Bin在BW's Bin内，执行下述操作
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()
