之前找到的别人的代码,感觉多余无用的地方太多了,重写了一版很简单明了的,附了详细的注释解析
注意不是用于直接解析数字码的,而是解析数字码经过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=()/
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 = % dB \nENOB = % 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 - ((Dout) + (Dout)) / 2
numpt2 = int((FFTN / 2))
Window = (FFTN)
Doutw = Dout * Window
DoutFFT = (Doutw, FFTN)
PSD = ((DoutFFT)) ** 2
FinBin = (PSD[:numpt2]) # 获取最大值索引
span = 3
DcPower = (PSD[0:span])
if FinBin > span:
SignalPower = (PSD[FinBin - span:FinBin + span])
else:
SignalPower = 0
SignalHarmonyNumbers = int((FFTN/2/FinBin))
SignalHarmonyFrequency = (SignalHarmonyNumbers)
SignalHarmonyPower = (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] = (PSD[HarmonyBin - span:HarmonyBin + span])
SignalHarmonyFrequency[i] = (HarmonyBin) * Fs / FFTN
AllPower = (PSD[0:numpt2])
NoisePower = AllPower - DcPower - SignalPower
SNDR = 10 * (SignalPower / NoisePower)
ENOBloc = (SNDR - ) /
HarmonyPower = (SignalHarmonyPower[1:SignalHarmonyNumbers]) - SignalHarmonyPower[0]
THD = 10 * (((SignalHarmonyPower[1:]**2)) / SignalHarmonyPower[0])
print(THD)
SFDR = 10 * (SignalHarmonyPower[0] / (SignalHarmonyPower[1:SignalHarmonyNumbers]))
print(SFDR)
(1)
(Dout[:FFTN])
('Sample Point')
('Dout')
('Dout Signal')
()
maxPower = (PSD)
PSD_dB = 10 * (PSD / maxPower)
(2)
((numpt2) * Fs / FFTN, PSD_dB[:numpt2], '-')
(Fs / 3, -20, f'SNDR = {SNDR:.1f} dB \nENOB = {ENOBloc:.2f} bits',
bbox=dict(facecolor='red', edgecolor='red', linewidth=3, pad=10))
('Normalized Power Spectrum')
('Frequency (Hz)')
('Normalized Power (dB)')
()
()
Copyright 2015-2022 华中质量网 版权所有 备案号:京ICP备12018864号-26 联系邮箱:2 913 236 @qq.com