请各位大侠帮我看看这段程序
%Example program routine to generate FFT plots and determine the dynamic performance of%a high-speed dataconverter from the data records taken with a HP16500C Logic Analyzer
%System. Data was extracted through the HPIB interface and read into the following MATLAB
%program routine. The same data can be extracted from the controller interface of the LA
%and simply be copied to a floppy disk—a rather time-consuming way, but possible.
%Start MAX1448 Dynamic Performance Test Routine
disp('HP16500C State Card');
filename=input('Type a:\filename or Press RETURN for HPIB Data Transfer: ');
if isempty(filename), filename = 'data.m'; end
fid=fopen(filename,'r');
numpt=input('Data Record Size (Number of Points)? ');
fclk=input('Sampling Frequency (MHz)? ');
%MAX1448 - 10-bit data converter
numbit=10;
%Discard first 13 lines from the data file, which do not contain data
for i=1:13, fgetl(fid); end
=fscanf(fid,'%f',); fclose(fid);
v1=v1'; code=v1(:,2);
%Display a warning, when the input generates a code greater than full-scale
if (max(code)==2^numbit-1) | (min(code)==0), disp('Warning: ADC may be clipping!!!'); end
%Plot results in the time domain
figure;
plot(,code); title('TIME DOMAIN'); xlabel('SAMPLES'); ylabel('DIGITAL OUTPUT CODE');
%Recenter the digital sine wave
Dout=code-(2^numbit-1)/2;
%If no window function is used, the input tone must be chosen to be unique and with
%regard to the sampling frequency. To achieve this prime numbers are introduced and the
%input tone is determined by fIN = fSAMPLE * (Prime Number / Data Record Size).
%To relax this requirement, window functions such as HANNING and HAMING (see below) can
%be introduced, however the fundamental in the resulting FFT spectrum appears 'sharper'
%without the use of window functions.
Doutw=Dout;
%Doutw=Dout.*hanning(numpt);
%Doutw=Dout.*hamming(numpt);
%Performing the Fast Fourier Transform
Dout_spect=fft(Doutw);
%Recalculate to dB
Dout_dB=20*log10(abs(Dout_spect));
%Display the results in the frequency domain with an FFT plot
figure; maxdB=max(Dout_dB(1:numpt/2));
%For TTIMD, use the following short routine, normalized to —6.5dB full-scale.
%plot(.*fclk/numpt,Dout_dB(1:numpt/2)-maxdB-6.5);
plot(.*fclk/numpt,Dout_dB(1:numpt/2)-maxdB);
grid on; a1=axis; axis();
title('FFT PLOT'); xlabel('ANALOG INPUT FREQUENCY (MHz)'); ylabel('AMPLITUDE (dB)');
%Calculate SNR, SINAD, THD and SFDR values
%Find the signal bin number, DC = bin 1
fin=find(Dout_dB(1:numpt/2)==maxdB);
%Span of the input frequency on each side
span=max(round(numpt/200),5);
%Approximate search span for harmonics on each side
spanh=2;
%Determine power spectrum
spectP=(abs(Dout_spect)).*(abs(Dout_spect));
%Find DC offset power
Pdc=sum(spectP(1:span));
%Extract overall signal power
Ps=sum(spectP(fin-span:fin+span));
%Vector/matrix to store both frequency and power of signal and harmonics
Fh=[];
%The 1st element in the vector/matrix represents the signal, the next element represents
%the 2nd harmonic, etc.
Ph=[];
%Find harmonic frequencies and power components in the FFT spectrum
for har_num=1:10
%Input tones greater than fSAMPLE are aliased back into the spectrum
tone=rem((har_num*(fin-1)+1)/numpt,1);
if tone>0.5
%Input tones greater than 0.5*fSAMPLE (after aliasing) are reflected
tone=1-tone;
end
Fh=;
%For this procedure to work, ensure the folded back high order harmonics do not overlap
%with DC or signal or lower order harmonics
har_peak=max(spectP(round(tone*numpt)-spanh:round(tone*numpt)+spanh));
har_bin=find(spectP(round(tone*numpt)-spanh:round(tone*numpt)+spanh)==har_peak);
har_bin=har_bin+round(tone*numpt)-spanh-1;
Ph=;
end
%Determine the total distortion power
Pd=sum(Ph(2:5));
%Determine the noise power
Pn=sum(spectP(1:numpt/2))-Pdc-Ps-Pd;
format;
A=(max(code)-min(code))/2^numbit
AdB=20*log10(A)
SINAD=10*log10(Ps/(Pn+Pd))
SNR=10*log10(Ps/Pn)
disp('THD is calculated from 2nd through 5th order harmonics');
THD=10*log10(Pd/Ph(1))
SFDR=10*log10(Ph(1)/max(Ph(2:10)))
disp('Signal & Harmonic Power Components:');
HD=10*log10(Ph(1:10)/Ph(1))
%Distinguish all harmonics locations within the FFT plot
hold on;
plot(Fh(2)*fclk,0,'mo',Fh(3)*fclk,0,'cx',Fh(4)*fclk,0,'r+',Fh(5)*fclk,0,'g*',Fh(6)*fclk,0,'bs',Fh(7)*fclk,0,'bd',Fh(8)*fclk,0,'kv',Fh(9)*fclk,0,'y^');
legend('1st','2nd','3rd','4th','5th','6th','7th','8th','9th');
hold off;
各位大侠，我对此程序有点不太懂，不知那位大侠帮我标注一下，不胜感激！！！
[ 本帖最后由 ChaChing 于 2009-5-31 10:59 编辑 ]
回复 楼主 maxuewei2004 的帖子
那一点不太懂? 标出来!?建议自己喂资料一行一行执行, 不是就知道吗? 我是看到标题想来学习下的
刚好之前也看到了一段类似的程序
其中也有下面的语句
span=max(round(numpt/200),5);
一直不甚明白
应该是计算谐波的有效频率范围吧
numpt/200，200是怎么来的？
5又代表什么
[ 本帖最后由 ChaChing 于 2009-5-14 13:19 编辑 ] max(round(numpt/200),5); 取两数round(numpt/200)及5的最大值 谢谢chaching的回复。
我的理解，这句话应该是解释计算信号功率时，除了计算幅值最大的频率点外，左右还要加上几个点。
要加的点的个数就是max(round(numpt/200),5);
但是我不知道两个数中除5之外，另一个数是的含义是什么？为什么要除以200？
请高手帮忙解惑下，谢谢！
[ 本帖最后由 ChaChing 于 2009-5-31 11:05 编辑 ]
回复 5楼 swifthy 的帖子
这已经是专业问题了!程序又长, 时间及个人水平专业有限
待高人路过 还是谢谢chaching
等待高手路过，我自己也钻研下！
页:
[1]