做ADC的动态性能时仿真遇到的问题，求助求助！

amazion

目前正在做ADC的动态特性仿真，计算SNR、THD等指标。从网上找到一个程序，如下所示。


file_input=('4.txt');
format long;
numpt=4096;
fclk=100;
x=load(file_input);
numpt2=numpt/2;
numbit=12;
FS=(2^numbit-1)/2;
length=numpt;
bias=0;
dc=(max(x)+min(x))/2+bias;
x=x-dc;
np=numpt;
np2=np/2;
Dout=x;
Doutw=Dout.*blackman(numpt);%°嬀/size]
Dout_spect=fft(Doutw);%±嬀/size]
spectP=(abs(Dout_spect)).*(abs(Dout_spect));%óù嬀/size]
Dout_dB=10*log10(spectP);%×
[maxdB,fin]=max(Dout_dB(1:numpt2));%×°·§×ó×ó嬀/size]
span=max(round(numpt/200),5);
spanh=2;
Pdc=sum(spectP(1:span));
if fin>span
    Ps=sum(spectP(fin-span:fin+span));
else
    Ps=sum(spectP(1:fin+span));
end
Fh=zeros(20);%
Ph=zeros(20);%
Ph(1)=sum(spectP(fin-1:fin+1));%计算基波
Fh(1)=fclk*(fin-1)/numpt;%计算信号频率
har_step=fin-1;%设定寻找高次谐波的间隔
har=fin;
for har_num=1:10    %取10次谐波
    har_ind=har_step*har_num+fin;    %计算出每个高阶谐波所在位置

if(har_ind<numpt2)
        [P_har har]=max(spectP(har_ind-spanh:har_ind+spanh));
        har=har_ind+har-spanh-1;

if har_num==1,
            bas_point=har;

end
        Ph(har_num+1)=sum(spectP(har-1:har+1));

Fh(har_num+1)=(har-1)*fclk/numpt;

end
end
f=Fh;
fsig=Fh(1)
har_pow_num=min(5,floor(fclk/(2*fsig)))
Pd=sum(Ph(1:har_pow_num))-Ph(1);
Pn=sum(spectP(1:numpt2))-Pdc-Ps-Pd;
har_range=min(10,floor(fclk/(2*fsig)));
A=(max(x)-min(x))/(2*FS);
AdB=20*log10(A)
SINAD=10*log10(Ps/(Pn+Pd))
SNR=10*log10(Ps/Pn)
THD=10*log10(Pd/Ph(1))
SFDR=10*log10(Ph(1)/max(Ph(2:har_range)))
ENOB=(SINAD-1.76)/6.02
figure;
subplot(211),plot([1:length],x);
a1=axis;
title('TIME DOMAIN')
xlabel('SAMPLES');
ylabel('DIGITAL OUTPUT CODE');
y_up=100*(floor(max(x)/100)+1);
y_down=0-y_up;
axis([a1(1) a1(2) y_down y_up]);
subplot(212),plot([0:numpt2-1]*fclk/numpt,Dout_dB(1:numpt2)-maxdB);
s2=sprintf('SNR   = %4.1fdB\n',SNR);
s3=sprintf('SINAD = %4.1fdB\n',SINAD);
s4=sprintf('ENOB = %4.1fBit\n',ENOB);
s5=sprintf('THD = %4.1fdB\n',THD);
s1=sprintf('SFDR = %4.1fdB\n',SFDR);
text(37,-20,s2);
text(37,-35,s3);
text(37,-50,s4);
text(37,-65,s5);
text(37,-80,s1);
title('FREQ.DOMAIN');
xlabel('ADC FFT SPECTRUM');
ylabel('AMPLITUDE(dB)');






在进行低频信号计算时，还算可以。但是计算奈奎斯特频率下，THD会出现溢出问题。
如输入信号频率为49.99999MHz时，则har_pow_num=1，找不到高阶谐波。则Pd=Ph(1)-Ph(1)=0.


求助，在奈奎斯特频率下，如何计算THD？
多谢多谢！

csmclee

哈哈，这个程序怎么看着这么眼熟啊，感觉和我以前写的程序很接近啊。有几个地方说一下：
1.如果不是过采样ADC，采4096点完全没必要，普通的NYQUIST AD/DA，如果遵循相干采样原则，64点就足够准确了，64点和4096点计算出来的SNR出入也就1,2个dB,点数多不过频谱看起来漂亮一点而已（注意前提是采用相干采样，同时电路初始要留足够的dummy time，我个人习惯是采80个点，用最后的64个点做FFT）；
2.同上，仿真如果采用相干采样，不会有频谱泄漏问题，无需加窗，不过测试里加blackman窗效果远比教科书推荐的hanning好；
3.按定义，一般是在你关心的频率范围内（FS/2），将最大的前9个谐波相加用于计算THD，这就出现你说的问题，输入信号较低时，THD计算没问题，但输入信号频率接近FS/2时，基频的谐波分量已经处于NYQUIST频带之外，由此计算出来的THD就会无穷小，最简单的办法就是忽略它不计，如果你实在想计算折射回FS/2以内的谐波能量，可以参见下图： ，不过我个人偷懒，也没有将这种谐波写进matlab程序，简单地说说方法吧：假设相干系数为64/3，意味着你用64个点采样3个正弦波，FFT之后，你的基波频率就处于第三个点（NYQUIST频率或者说FS/2频率就是第32个点），由此你关心的基频的谐波分量分别位于第6，第9，第12……到第30个点；
假设相关系数为31/64，则代表用64个点采样31个正弦波，输入信号频率已经接近NYQUIST频率，频谱上信号频率处于31点位置，而它相应谐波能量如上图会折射回FS/2频带内，你只需把前30个点中最大的九个加起来就是你希望得到的THD了。
希望我的回答可以给你一点帮助，呵呵。

amazion

回复 2# csmclee


   给俺讲的相当详细，多谢了。后来我并没有继续用MATLAB进行计算，而是用了国半的一个软件WAVEVISON直接进行了计算，结果类似。多谢了！

yongxie

s2=sprintf('SNR   = %4.1fdB\n',SNR);
s3=sprintf('SINAD = %4.1fdB\n',SINAD);
s4=sprintf('ENOB = %4.1fBit\n',ENOB);
s5=sprintf('THD = %4.1fdB\n',THD);
s1=sprintf('SFDR = %4.1fdB\n',SFDR);
以上几行能解释下吗？%4.1fdB是个什么输出格式？
先谢谢各位大牛了

amazion

回复 4# yongxie


    这是调整输出格式的语句

ah_yy

你好，我最近设计了一个10bit 的SSADC，但是不清楚如何仿真INL DNL SNR参数，查了一些资料，多数都是指导用matlAB，但是小弟都不清楚数据文件“4.txt”都包括什么，如果可以的话 希望你能指点一下~

will199321

学习了

chafenfa

学习了