用matlab做adc的动态参数仿真

kongc

本人最近在研究ADC，现在是用matlab做fft分析，用网上找的代码加上自己的数据，运行后频谱图结果是出来了，但是那些动态参数一个都没有计算出来，而且都提示“下标索引必须为正整数类型或逻辑类型。”，求大佬指点迷津。

kongc

代码如下：

1. %***********************************************************************%
   
2. % The following program code plots the FFT spectrum of a desired test tone.
   
3. % Test tone based on coherent sampling criteria, and computes SNR, SNDR, THD and SFDR.
   
4. % This program is believed to be accurate and reliable.
   
5. % This program may get altered without prior notification.;
   
6. % Company: xxx analog Mixed-signal Group
   
7. % Author:
   
8. % Date: 2016-01-09
   
9. %***********************************************************************%
   
10. 
    
11. 
    
12. 
    
13. 
    
14. %***********************************************************************%
    
15. % 采样时钟240MHz，输入信号频率86.25MHz，FFT点数64，ADC分辨率10bit，仿真时间400ns。
    
16. % 从cadence输出时间隔是4.16666667ns，共计97个点，从第20个点开始取64点做FFT。
    
17. % cadence输出的是归一化的结果范围是0-1.6V，需除以1.6再乘以1024，转换为数字码。
    
18. %***********************************************************************%
    
19. 
    
20. 
    
21. 
    
22. 
    
23. clear all;
    
24. clc;
    
25. datafile='ADC_result8.csv'
    
26. spectP_file='spec.csv';
    
27. 
    
28. 
    
29. 
    
30. 
    
31. %***********************************************************************%
    
32. % 输入采样时钟、样本点数、分辨率等变量
    
33. %***********************************************************************%
    
34. fs=240e6;        %采样时钟
    
35. Data_Num=64;   %样本点数
    
36. numbit=10;       %ADC分辨率
    
37. data_start=20; %取点起始位置
    
38. fclk=fs/1e6;     %x坐标轴数值显示
    
39. numpt=Data_Num;
    
40. fres=fclk/numpt; %Desired frequency resolution of FFT[MHz], fres=fclk/2^N
    
41. 
    
42. 
    
43. %***********************************************************************%
    
44. % 读取数据
    
45. %***********************************************************************%
    
46. d_in=csvread('ADC_result8.csv',1,1);
    
47. d_in=d_in/1.6*1024;
    
48. code=zeros(1,numpt);
    
49. code(1:numpt)=d_in(data_start:data_start+numpt-1);
    
50. 
    
51. 
    
52. %***********************************************************************%
    
53. % Plot output code
    
54. %***********************************************************************%
    
55. figure;
    
56. plot(code);
    
57. title(sprintf('ADC Digital Output'));
    
58. 
    
59. 
    
60. %***********************************************************************%
    
61. % Recenter the digital sine wave, for 2's complement code
    
62. %***********************************************************************%
    
63. m_ean=mean(code);
    
64. for hk=1:length(code)
    
65.     code(hk)=code(hk)-m_ean;
    
66. end
    
67. 
    
68. 
    
69. %***********************************************************************%
    
70. % Display a warning, when the input generates a code greater than
    
71. % full-scale, for 2's complement code
    
72. %***********************************************************************%
    
73. max_code=max(code)
    
74. min_code=min(code)
    
75. if (max(code)>2^(numbit-1)) | (min(code)<(0-2^(numbit-1)))
    
76. %if (max(code)==2numbit -1) | (min(code)==0)
    
77.     disp('Warning: ADC may be clipping!');
    
78. end
    
79. 
    
80. 
    
81. %***********************************************************************%
    
82. % Normalize input signal relative to full-scale
    
83. %***********************************************************************%
    
84. fin_dB=20*log10((max_code-min_code)/(2^numbit));
    
85. 
    
86. 
    
87. %***********************************************************************%
    
88. % 对数据样本加窗函数处理
    
89. %***********************************************************************%
    
90. % If no window function is used, the input tone must be chosen to be unique and with
    
91. % regard to the sampling frequency. To achieve this prime numbers are introduced and the
    
92. % input tone is determined by Fin = Fsample * (Prime Number / Data Record Size).
    
93. % To relax this requirement, window functions such as HANNING and HAMING (see below) can
    
94. % be introduced, however the fundamental in the resulting FFT spectrum appears 'sharper'
    
95. %without the use of window functions.
    
96. Dout=code';
    
97. Doutw=Dout;
    
98. % Doutw=Dout.*hanning(numpt);
    
99. % Doutw=Dout.*hamming(numpt);
    
100. % Doutw=Dout.*blackman(numpt);
     
101. 
     
102. 
     
103. %***********************************************************************%
     
104. % Performing the Fast Fourier Transform [FFT]
     
105. %***********************************************************************%
     
106. span=0; %Span of the input frequency on each side; span=max(round(numpt/200),5);
     
107. spanh=0; %Approximate search span for harmonics on each side
     
108. spandc=0; %Approximate search span for DC on right side
     
109. Dout_spect=fft(Doutw);
     
110. Dout_dB=20*log10(abs(Dout_spect)); %Recalculate to dB abs(Dout_spect)
     
111. spectP=(abs(Dout_spect)).*(abs(Dout_spect)); %Determine power spectrum
     
112. maxdB=max(Dout_dB(1+spandc:numpt/2));
     
113. fin=find(Dout_dB(1:numpt/2)==maxdB); %Find the signal bin number, DC=bin1
     
114. 
     
115. 
     
116. %***********************************************************************%
     
117. % Calculate SNR, SNDR, THD and SFDR values.
     
118. %***********************************************************************%
     
119. fw=fopen(spectP_file,'w'); %write the power soectrum to file
     
120. fprintf(fw,'%12.9e\n',spectP);
     
121. fclose('all');
     
122. %Find DC offset power
     
123. Pdc=sum(spectP(1:span));
     
124. %Extract overall signal power
     
125. idx1=fin-span;
     
126. idx2=fin+span;
     
127. if(idx1<=0)
     
128.     idx1 = 1;
     
129. end
     
130. Ps=sum(spectP(idx1:idx2));
     
131. %Vector/matrix to store both frequency and power of signal and harmonics
     
132. Fh=[];
     
133. %The 1st element in the vector/matrix represents the signal, the next element represents the 2nd harmonic
     
134. Ph=[];
     
135. %Vector/matrix to store the sampling points responding to the harmonics
     
136. Nh=[];
     
137. %Ah represents signal and harmonic amplitude
     
138. Ah=[];
     
139. %Find harmonic frequencies and power components in the FFT spectrum
     
140. %For this procedure to work, ensure the folded back high order harmonics do not overlap
     
141. %with DC or signal or lower order harmonics , so it should be modified according to the actual condition
     
142. for har_num=1:5
     
143.     tone=rem((har_num*(fin-1)+1)/numpt,1); %Input tones greater than fSAMPLE are aliased back into the spectrum
     
144.     if tone>0.5
     
145.         tone=1-tone; %Input tones greater than 0.5*Fsample (after aliasing) are reflected
     
146.     end
     
147.     Fh=[Fh tone];
     
148.     %Check Nh to see the bin of the harmonics
     
149.     Nh=[Nh round(tone*numpt)];
     
150.     %For this procedure to work, ensure the folded back high order harmonics do not overlap
     
151.     %with DC or signal or lower order harmonics
     
152.     har_peak=max(spectP(round(tone*numpt)-spanh:round(tone*numpt)+spanh));
     
153.     har_bin=find(spectP(round(tone*numpt)-spanh:round(tone*numpt)+spanh)==har_peak);
     
154.     har_bin=har_bin+round(tone*numpt)-spanh-1;
     
155.     Ph=[Ph sum(spectP(har_bin-spanh:har_bin+spanh))];
     
156.     Ah=[Ah Dout_dB(har_bin)];
     
157. end
     
158. %Determine the total distortion power, it should be modified according to the actual condition.
     
159. Pd=sum(Ph(2:5));
     
160. %Determine the noise power
     
161. Pn=sum(spectP(1:numpt/2))-Pdc-Ps-Pd;
     
162. %Determine the next largest component
     
163. spur_max=max(max(spectP(spandc+1:fin-span-1)),max(spectP(fin+span+1:numpt/2)));
     
164. spur_bin=find(spectP(1:numpt/2)==spur_max)
     
165. %**********************计算动态特性结果**********************%
     
166. format; %设置输出格式
     
167. SFDR = 10*log10(max(spectP(1:numpt/2))/spur_max); %-fin_dB
     
168. THD = 10*log10(Pd/Ps); %+fin_dB
     
169. SNR = 10*log10(Ps/Pn); %-fin_dB
     
170. SNDR = 10*log10(Ps/(Pn+Pd)); %-fin_dB
     
171. ENOB = (SNDR-1.76)/6.02;
     
172. %disp('Note: THD is calculated from 2nd through 10th order harmonics.');
     
173. %*********************标示信号和谐波位置*********************%
     
174. %hold on;
     
175. %plot((Nh(2:10)-1).*fres,Ah(2:10)-maxdB+fin_dB,'rs');
     
176. % 标示信号
     
177. bins=(Nh(1)-1)*fres;
     
178. Ahs=Ah(1)-maxdB+fin_dB;
     
179. % 标示2次谐波
     
180. bin2=(Nh(2)-1)*fres;
     
181. Ah2=Ah(2)-maxdB+fin_dB;
     
182. % 标示3次谐波
     
183. bin3=(Nh(3)-1)*fres;
     
184. Ah3=Ah(3)-maxdB+fin_dB;
     
185. % 标示4次谐波
     
186. bin4=(Nh(4)-1)*fres;
     
187. Ah4=Ah(4)-maxdB+fin_dB;
     
188. % 标示5次谐波
     
189. bin5=(Nh(5)-1)*fres;
     
190. Ah5=Ah(5)-maxdB+fin_dB;
     
191. % 在FFT频谱图中追加标示
     
192. figure;
     
193. plot(bins,Ahs,'rs',bin2,Ah2,'rd',bin3,Ah3,'r^',bin4,Ah4,'r*',bin5,Ah5,'rx');
     
194. legend('SINGAL','HD2','HD3','HD4','HD5');
     
195. %**********************图表显示**********************%
     
196. ylabel('Full-Scale Normalized Magnitude[dB]')
     
197. xlabel('Frequency [MHz]')
     
198. title(sprintf('ADC FFT Spectrum (%g points)\nFs = %g MSps, Fin = %g MHz (%1.2gdBFS)', Data_Num,fs/1e6,(fin-1)*fres,fin_dB));
     
199. grid on;
     
200. box on;
     
201. ylim([-110 10]);
     
202. set(gca,'xgrid', 'off');
     
203. set(gca, 'GridLineStyle' ,'-');
     
204. set(gca,'yTick',[-110:10:10]);
     
205. %****************************************************%
     
206. %Display the results in the frequency domAIn with an FFT plot.
     
207. for i=0:1:numpt/2-1)
     
208.     hold on;
     
209.     line([i*fres,i*fres],[-110,Dout_dB(i+1)-maxdB+fin_dB],'LineWidth',2);
     
210.     hold off;
     
211. end
     
212. %***********************在图中打印结果***********************%
     
213. hold on;
     
214. s1=sprintf('SFDR = %4.1fdB\n',SFDR);
     
215. s2=sprintf('THD = %4.1fdB\n',THD);
     
216. s3=sprintf('SNR   = %4.1fdB\n',SNR);
     
217. s4=sprintf('SNDR = %4.1fdB\n',SNDR);
     
218. s5=sprintf('ENOB = %4.2fbit\n',ENOB);
     
219. text(25,-10,s1);
     
220. text(25,-20,s2);
     
221. text(25,-30,s3);
     
222. text(25,-40,s4);
     
223. text(25,-50,s5);
     
224. hold off;
     

复制代码

fightshan

没有必要用网上找的matlab代码，直接用波形窗口的频谱分析工具不香吗？

kongc

那个频谱分析工具可以算ADC的动态参数吗？我之前也研究过，但好像就是因为只有频谱图没有参数的计算才考虑代码的

kongc

原因找到了，下列代码中的fin值过小。fft后的dout_dB的值前两个都很大，导致fin值很小，后面的矩阵下标都是负值了，所以改了spandc的值把前两个dout_dB忽略掉就正常了。但是这样我有两个问题：一是这样忽略会不会有问题？二是为什么fft后前两个点的幅度怎么这么大，明显比后面的值大很多。

1. Dout_spect=fft(Doutw);
   
2. Dout_dB=20*log10(abs(Dout_spect)); %Recalculate to dB abs(Dout_spect)
   
3. spectP=(abs(Dout_spect)).*(abs(Dout_spect)); %Determine power spectrum
   
4. maxdB=max(Dout_dB(1+spandc:numpt/2));
   
5. fin=find(Dout_dB(1:numpt/2)==maxdB);
   

复制代码

GSHANG

顶一个，遇到过同样的问题，增大输入信号频率之后解决了

jade_ysz

楼主你好赞。现在刚刚做到这里

jade_ysz

楼主，你好，请问一下代码第26行的文件是什么文件啊，我只有一个数据文件

17315768922

jade_ysz 发表于 2021-11-24 20:26
楼主，你好，请问一下代码第26行的文件是什么文件啊，我只有一个数据文件 ...

你好，请问一下，你问的这个文件是什么文件

linda313417

楼主你好，请问spectP_file='spec.csv';这行代码中的文件是从哪来的？是什么文件呢？