adc的fft仿真

kongc

最近在做一个saradc的fft仿真，在网上找的matlab代码，加上自己的数据，发现得到的频谱图很离谱，有没有大佬有兴趣拿着代码和数据试一下，帮我看看哪里出错了？，谢谢！

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. close all;
    
25. %spectP_file='spec.csv';
    
26. 
    
27. 
    
28. 
    
29. 
    
30. %***********************************************************************%
    
31. % 输入采样时钟、样本点数、分辨率等变量
    
32. %***********************************************************************%
    
33. fs=100e6;        %采样时钟
    
34. Data_Num=32768;   %样本点数
    
35. numbit=12;       %ADC分辨率
    
36. data_start=1; %取点起始位置
    
37. fclk=fs/1e6;     %x坐标轴数值显示
    
38. numpt=Data_Num;
    
39. fres=fclk/numpt; %Desired frequency resolution of FFT[MHz], fres=fclk/2^N
    
40. 
    
41. 
    
42. %***********************************************************************%
    
43. % 读取数据
    
44. %***********************************************************************%
    
45. %d_in=dlmread('E:\kongcaho\ADC\MATLAB\adc_data_0715.csv','%s')';
    
46. %d_in=d_in*2^numbit;
    
47. 
    
48. d_in=csvread('E:\kongcaho\ADC\MATLAB\adc_data.csv');%以字符形式打开文件
    
49. d11=1;d10=2;d9=3;d8=4;d7=5;d6=6;d5=7;d4=8;d3=9;d2=10;d1=11;d0=12;%标示样本矩阵的列
    
50. vref=1;
    
51. %adc_data=d_in;
    
52. adc_data=round(d_in);
    
53. 
    
54. for i=1:numpt
    
55.     dout(i)=adc_data(i,d11)/2+adc_data(i,d10)/4+adc_data(i,d9)/8 ...
    
56.     +adc_data(i,d8)/16+adc_data(i,d7)/32+adc_data(i,d6)/64 ...
    
57.     +adc_data(i,d5)/128+adc_data(i,d4)/254+adc_data(i,d3)/512 ...
    
58.     +adc_data(i,d2)/1024+adc_data(i,d1)/2048+adc_data(i,d0)/4096;%通过转换器出来的结果，恢复原来的波形
    
59.     doute(i)=dout(i)-0.5;%将正弦波的共模电平偏置为0
    
60.     %doute(i)=dout(i)-mean(dout(i));
    
61. end
    
62. 
    
63. 
    
64. code=zeros(1,numpt);  %创建一个全为0元素的数组
    
65. %for i=1:1:numpt
    
66. %code(1:numpt)=d_in(data_start:data_start+numpt-1);
    
67. %end
    
68. for i=1:numpt
    
69.     code(i)=adc_data(i,d11)*2048+adc_data(i,d10)*1024+adc_data(i,d9)*512 ...
    
70.     +adc_data(i,d8)*256+adc_data(i,d7)*128+adc_data(i,d6)*64+adc_data(i,d5)*32 ...
    
71.     +adc_data(i,d4)*16+adc_data(i,d3)*8+adc_data(i,d2)*4+adc_data(i,d1)*2+adc_data(i,d0)*1;%通过转换器出来的结果，恢复原来的波形
    
72.    
    
73. end
    
74. 
    
75. dcadc(code,numpt);
    
76. 
    
77. %***********************************************************************%
    
78. % Plot output code
    
79. %***********************************************************************%
    
80. figure;
    
81. plot(code);
    
82. title(sprintf('ADC Digital Output'));
    
83. 
    
84. 
    
85. figure;
    
86. plot(dout);
    
87. title(sprintf('ADC ANALOG Output'));
    
88. 
    
89. 
    
90. %***********************************************************************%
    
91. % Recenter the digital sine wave, for 2's complement code
    
92. %***********************************************************************%
    
93. %m_ean=mean(code);   %求code数组的平均值
    
94. %for hk=1:length(code)
    
95. %    code(hk)=code(hk)-m_ean;
    
96. %end
    
97. 
    
98. 
    
99. %***********************************************************************%
    
100. % Display a warning, when the input generates a code greater than
     
101. % full-scale, for 2's complement code
     
102. %***********************************************************************%
     
103. max_code=max(code)
     
104. min_code=min(code)
     
105. if (max(code)==2^numbit-1) | (min(code)==0)
     
106. %if (max(code)==2^numbit -1) | (min(code)==0)
     
107.     disp('Warning: ADC may be clipping!');
     
108. end
     
109. 
     
110. 
     
111. %***********************************************************************%
     
112. % Normalize input signal relative to full-scale
     
113. %***********************************************************************%
     
114. fin_dB=20*log10((max_code-min_code)/(2^numbit));
     
115. 
     
116. 
     
117. %***********************************************************************%
     
118. % 对数据样本加窗函数处理
     
119. %***********************************************************************%
     
120. % If no window function is used, the input tone must be chosen to be unique and with
     
121. % regard to the sampling frequency. To achieve this prime numbers are introduced and the
     
122. % input tone is determined by Fin = Fsample * (Prime Number / Data Record Size).
     
123. % To relax this requirement, window functions such as HANNING and HAMING (see below) can
     
124. % be introduced, however the fundamental in the resulting FFT spectrum appears 'sharper'
     
125. %without the use of window functions.
     
126. Dout=doute';
     
127. %Doutw=Dout;
     
128. Doutw=Dout.*hanning(numpt);
     
129. % Doutw=Dout.*hamming(numpt);
     
130. % Doutw=Dout.*blackman(numpt);
     
131. 
     
132. 
     
133. %***********************************************************************%
     
134. % Performing the Fast Fourier Transform [FFT]
     
135. %***********************************************************************%
     
136. %span=5;     %Span of the input frequency on each side;
     
137. span=max(round(numpt/200),5);
     
138. spanh=2;    %Approximate search span for harmonics on each side
     
139. spandc=12;   %Approximate search span for DC on right side
     
140. Dout_spect=fft(Doutw);
     
141. Dout_dB=20*log10(abs(Dout_spect)); %Recalculate to dB abs(Dout_spect)
     
142. spectP=(abs(Dout_spect)).*(abs(Dout_spect)); %Determine power spectrum
     
143. maxdB=max(Dout_dB(1+spandc:numpt/2));
     
144. fin=find(Dout_dB(1:numpt/2)==maxdB); %Find the signal bin number, DC=bin1
     
145. %fin=200;
     
146. 
     
147. %***********************************************************************%
     
148. % Calculate SNR, SNDR, THD and SFDR values.
     
149. %***********************************************************************%
     
150. %fw=fopen(spectP_file,'w'); %write the power soectrum to file
     
151. %fprintf(fw,'%12.9e\n',spectP);
     
152. %fclose('all');
     
153. %Find DC offset power
     
154. Pdc=sum(spectP(1:span));
     
155. %Extract overall signal power
     
156. idx1=fin-span;
     
157. idx2=fin+span;
     
158. if(idx1<=0)
     
159.     idx1 = 1;
     
160. end
     
161. Ps=sum(spectP(idx1:idx2));    %信号功率
     
162. %Vector/matrix to store both frequency and power of signal and harmonics
     
163. Fh=[];
     
164. %The 1st element in the vector/matrix represents the signal, the next element represents the 2nd harmonic
     
165. Ph=[];
     
166. %Vector/matrix to store the sampling points responding to the harmonics
     
167. Nh=[];
     
168. %Ah represents signal and harmonic amplitude
     
169. Ah=[];
     
170. %Find harmonic frequencies and power components in the FFT spectrum
     
171. %For this procedure to work, ensure the folded back high order harmonics do not overlap
     
172. %with DC or signal or lower order harmonics , so it should be modified according to the actual condition
     
173. for har_num=1:5
     
174.     tone=rem((har_num*(fin-1)+1)/numpt,1); %Input tones greater than fSAMPLE are aliased back into the spectrum
     
175.     if tone>0.5
     
176.         tone=1-tone; %Input tones greater than 0.5*Fsample (aft er aliasing) are reflected
     
177.     end
     
178.     Fh=[Fh tone];
     
179.     %Check Nh to see the bin of the harmonics
     
180.     Nh=[Nh round(tone*numpt)];
     
181.     %For this procedure to work, ensure the folded back high order harmonics do not overlap
     
182.     %with DC or signal or lower order harmonics
     
183.     har_peak=max(spectP(round(tone*numpt)-spanh:round(tone*numpt)+spanh));
     
184.     har_bin=find(spectP(round(tone*numpt)-spanh:round(tone*numpt)+spanh)==har_peak);
     
185.     har_bin=har_bin+round(tone*numpt)-spanh-1;
     
186.     Ph=[Ph sum(spectP(har_bin-spanh:har_bin+spanh))];
     
187.     Ah=[Ah Dout_dB(har_bin)];
     
188. end
     
189. %Determine the total distortion power, it should be modified according to the actual condition.
     
190. Pd=sum(Ph(2:5));                      %总谐波功率
     
191. %Determine the noise power
     
192. Pn=sum(spectP(1:numpt/2))-Pdc-Ps-Pd;  %噪声功率
     
193. %Determine the next largest component
     
194. spur_max=max(max(spectP(spandc+1:fin-span-1)),max(spectP(fin+span+1:numpt/2)));
     
195. spur_bin=find(spectP(1:numpt/2)==spur_max);
     
196. %**********************计算动态特性结果**********************%
     
197. format; %设置输出格式
     
198. SFDR = 10*log10(max(spectP(1:numpt/2))/spur_max); %-fin_dB
     
199. THD = 10*log10(Pd/Ps); %+fin_dB
     
200. SNR = 10*log10(Ps/Pn); %-fin_dB
     
201. SNDR = 10*log10(Ps/(Pn+Pd)); %-fin_dB
     
202. ENOB = (SNDR-1.76)/6.02;
     
203. %disp('Note: THD is calculated from 2nd through 10th order harmonics.');
     
204. %*********************标示信号和谐波位置*********************%
     
205. %hold on;
     
206. %plot((Nh(2:10)-1).*fres,Ah(2:10)-maxdB+fin_dB,'rs');
     
207. % 标示信号
     
208. bins=(Nh(1)-1)*fres;
     
209. Ahs=Ah(1)-maxdB+fin_dB;
     
210. % 标示2次谐波
     
211. bin2=(Nh(2)-1)*fres;
     
212. Ah2=Ah(2)-maxdB+fin_dB;
     
213. % 标示3次谐波
     
214. bin3=(Nh(3)-1)*fres;
     
215. Ah3=Ah(3)-maxdB+fin_dB;
     
216. % 标示4次谐波
     
217. bin4=(Nh(4)-1)*fres;
     
218. Ah4=Ah(4)-maxdB+fin_dB;   
     
219. % 标示5次谐波
     
220. bin5=(Nh(5)-1)*fres;
     
221. Ah5=Ah(5)-maxdB+fin_dB;
     
222. % 在FFT频谱图中追加标示
     
223. figure;
     
224. plot(bins,Ahs,'rs',bin2,Ah2,'rd',bin3,Ah3,'r^',bin4,Ah4,'r*',bin5,Ah5,'rx');
     
225. legend('SINGAL','HD2','HD3','HD4','HD5');
     
226. %**********************图表显示**********************%
     
227. ylabel('Full-Scale Normalized Magnitude[dB]')
     
228. xlabel('Frequency [MHz]')
     
229. title(sprintf('ADC FFT Spectrum (%g points)\nFs = %g MSps, Fin = %g MHz (%1.2gdBFS)', Data_Num,fs/1e6,(fin-1)*fres,fin_dB));
     
230. grid on;
     
231. box on;
     
232. ylim([-110 100]);
     
233. set(gca,'xgrid', 'off');
     
234. set(gca, 'GridLineStyle' ,'-');
     
235. set(gca,'yTick',[-110:10:100]);
     
236. %****************************************************%
     
237. %Display the results in the frequency domAIn with an FFT plot.
     
238. for i=0:1:(numpt/2-1)
     
239.     hold on;
     
240.     line([i*fres,i*fres],[-110,Dout_dB(i+1)-maxdB+fin_dB],'LineWidth',2);
     
241.     hold off;
     
242. end
     
243. %***********************在图中打印结果***********************%
     
244. hold on;
     
245. s1=sprintf('SFDR = %4.1fdB\n',SFDR);
     
246. s2=sprintf('THD = %4.1fdB\n',THD);
     
247. s3=sprintf('SNR   = %4.1fdB\n',SNR);
     
248. s4=sprintf('SNDR = %4.1fdB\n',SNDR);
     
249. s5=sprintf('ENOB = %4.2fbit\n',ENOB);
     
250. text(25,-10,s1);
     
251. text(25,-20,s2);
     
252. text(25,-30,s3);
     
253. text(25,-40,s4);
     
254. text(25,-50,s5);
     
255. hold off;
     

复制代码

1. function[DNL,INL,vo]=dcdac(vi,num);
   
2. %vi_in=dlmread('E:\kongcaho\ADC\MATLAB\adc_data_0715.csv','%s')';
   
3. %num=4096;
   
4. vi_f=vi;
   
5. 
   
6. format long;
   
7. bit=12;%ADC精度
   
8. n=2^bit;%ADC输出码数目
   
9. N_negetive=0;%高于0的采样总数
   
10. N_positive=0;%低于0的采样总数
    
11. N_tot=0;%总采样数目
    
12. A=0.5;%正弦波幅度
    
13. temp=vi_f ; %读入数
    
14. %-----------------------------------------------%
    
15. %计算每个码的统计次数
    
16. %-----------------------------------------------%
    
17. hh=zeros(n,1);
    
18. for i=1:2^bit
    
19.     count=0;%附值
    
20.     for j=1:num
    
21.         if temp(j)==i-1
    
22.             count=count+1;
    
23.         else count=count;
    
24.         end
    
25.     end
    
26.     hh(i)=count;
    
27. end
    
28. h=[hh(1:2^bit)];
    
29. LSB=A*2.0/(2^bit);%理想LSB
    
30. %-----------------------------------------------%
    
31. %计算offset：vo
    
32. %-----------------------------------------------%
    
33. N_negetive=sum(hh(1:n/2));
    
34. N_positive=sum(hh(n/2+1:n));
    
35. N_tot=sum(hh(1:n))
    
36. vo=A*pi/2*sin((N_positive-N_negetive)/N_tot)
    
37. %-----------------------------------------------%
    
38. %计算INL、DNL
    
39. %-----------------------------------------------%
    
40. ch=zeros(n,1);
    
41. v=zeros(n,1);%possibility of the ith code occurs
    
42. in=zeros(n,1);%每级台阶的INL
    
43. dn=zeros(n,1);%每级台阶的DNL
    
44. for i=1:n
    
45.     for j=1:i
    
46.         ch(i)=ch(i)+h(j);
    
47.     end
    
48.     v(i)=-cos(pi*ch(i)/N_tot);
    
49. end
    
50. LSB=(v(n-1)-v(2))/(n-3);%实际LSB
    
51. for i=1:n-1
    
52.     in(i)=v(i)/LSB-i+1024;
    
53.   %  in(i)=(v(i)-vo)/LSB-i+2048;
    
54.     dn(i)=(v(i+1)-v(i))/LSB-1;
    
55. end
    
56. INL=max(abs(in))
    
57. DNL=max(abs(dn))
    
58. %-----------------------------------------------%
    
59. %初始化转移曲线坐标
    
60. %-----------------------------------------------%
    
61. vv(1)=-A;
    
62. vv(2)=v(1);
    
63. xx(1)=0;
    
64. xx(2)=0;
    
65. vv(2*(n-1)+1)=v(n-1);
    
66. xx(2*(n-1)+1)=n-1;
    
67. for i=1:n-2
    
68.     vv(2*i+1)=v(i);
    
69.     vv(2*i+2)=v(i+1);
    
70.     xx(2*i+1)=i;
    
71.     xx(2*i+2)=i;
    
72. end
    
73. %-----------------------------------------------%
    
74. figure;
    
75. subplot(2,1,1);%画INL
    
76. plot(2:n-2,in(2:n-2),'-.rd');
    
77. xlabel('DIGITAL OUTPUT CODE');
    
78. ylabel('INL(LSB)');
    
79. title('INL vs.DIGITAL OUTPUT CODE');
    
80. grid on;
    
81. subplot(2,1,2);%画DNL
    
82. plot(1:n-3,dn(1:n-3));
    
83. xlabel('DIGITAL OUTPUT CODE');
    
84. ylabel('DNL(LSB)');
    
85. title('DNL vs.DIGITAL OUTPUT CODE');
    
86. grid on;
    
87. figure;
    
88. subplot(2,1,1);%画码密度柱状图
    
89. bar(1:n,hh,'g');
    
90. xlabel('DIGITAL OUTPUT CODE');
    
91. ylabel('Count of code');
    
92. title('Histogram Graph of Code Density');
    
93. axis([-100 2100 0 1000]);
    
94. grid on;
    
95. subplot(2,1,2);%画转移曲线
    
96. plot(vv(1:2*(n-1)),xx(1:2*(n-1)));
    
97. %plot(vv(1:n),xx(1:n));
    
98. xlabel('Input Voltage(V)');
    
99. ylabel('DIGITAL OUTPUT CODE');
    
100. title('Transfer Function');
     
101. axis([-1,1,0,2^bit-1]);
     
102. grid on;
     
103. 
     

复制代码

buxiangwanla

兄弟，你这什么原因找到了吗？