当前位置: 首页 > 时讯

用于ADC信号测试的PSD/FFT测试代码

发布时间:2023-08-06 17:28:47 来源:哔哩哔哩

之前找到的别人的代码,感觉多余无用的地方太多了,重写了一版很简单明了的,附了详细的注释解析

注意不是用于直接解析数字码的,而是解析数字码经过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