cal snr

lhlbluesky

下面是求SNR的代码，但是我实在是没看懂，望高人能指点一下，具体各个部分都是干什么的
function SNR=cal_snr(v1)
v1=v1';
code=v1(:,2);
numbit=14;
numpt=length(code);
fclk=50e6
%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([1:numpt],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([0:numpt/2-1].*fclk/numpt,Dout_dB(1:numpt/2)-maxdB-6.5);

%plot([0:numpt/2-1].*fclk/numpt,Dout_dB(1:numpt/2)-maxdB);
%grid on;
%title('FFT PLOT');
%xlabel('Analog INPUT FREQUENCY (MHz)');
%ylabel('AMPLITUDE (dB)');
%a1=axis; axis([a1(1) a1(2) -120 a1(4)]);

%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(2);
%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=[Fh tone];
%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=[Ph sum(spectP(har_bin-1:har_bin+1))];
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;

%Dynamic-Range Specifications, TTIMD

%Two-tone IMD can be a tricky measurement, because the additional equipment required (a power combiner to combine two input frequencies) can contribute unwanted intermodulation products that falsify the ADC's intermodulation distortion. You must observe the following conditions to optimize IMD performance, although they make the selection of proper input frequencies a tedious task.

%First, the input tones must fall into the passband of the input filter. If these tones are close together (several tens or hundreds of kilohertz for a megahertz bandwidth), an appropriate window function must be chosen as well. Placing them too close together, however, may allow the power combiner to falsify the overall IMD readings by contributing unwanted 2nd- and 3rd-order IMD products (depending on the input tones' location within the passband). Spacing the input tones too far apart may call for a different window type that has less frequency resolution.

%The setup also requires a minimum of three phase-locked signal generators. This requirement seldom poses a problem for test labs, but generators have different capabilities for matching frequency and amplitude. Compensating such mismatches to achieve (for example) a -0.5dB FS two-tone envelope and signal amplitudes of -6.5dB FS will increase your effort and test time (see the following program-code extraction).

%For TTIMD, use the following short routine, normalized to -6.5dB full-scale.
%plot([0:numpt/2-1].*fclk/numpt,Dout_dB(1:numpt/2)-maxdB-6.5);

%plot([0:numpt/2-1].*fclk/numpt,Dout_dB(1:numpt/2)-maxdB);
%grid on;
%title('FFT PLOT');
%xlabel('ANALOG INPUT FREQUENCY (MHz)');
%ylabel('AMPLITUDE (dB)');
%a1=axis; axis([a1(1) a1(2) -120 a1(4)]);

lhlbluesky

有没有 牛人 给指点 一二  

谢谢了先

mydreamhouse

这是MAXIM公司做AD测试的程序，
程序主要先将信号做DFT变换，再将信号功率和谐波功率以及噪声功率分别算出，其中难点就在谱泄漏的处理。MAXIM采用了一个近似方法做估计，即取基波和谐波周围SPAN范围内功率。

luyan

学习了啊

foolish

我先复制了，回头分析去。

emerlai

Yes, it calculate the spectrum leakage by hanning window first, then average it by get a span of spectrum.

raylee

在学习SNR中，看不懂！希望将来能彻底明白！

ruyi_278

谢谢共享，现在正需要 ！！

hustlsc

回复 1# lhlbluesky

我最近也在用这个代码计算snr等参量。但是遇到点问题，想请教下。就是我用这个代码计算出来的Pd（谐波能量）很大，大于我的信号能量，这样我在计算Pn的时候，用信号能量减去Ps-Pd的时候就会出现负数，不知道这是怎么回事，求解答。谢谢。

hustlsc

回复 3# mydreamhouse


   我最近也在用这个代码计算snr等参量。但是遇到点问题，想请教下。就是我用这个代码计算出来的Pd（谐波能量）很大，大于我的信号能量，这样我在计算Pn的时候，用信号能量减去Ps-Pd的时候就会出现负数，不知道这是怎么回事，求解答。谢谢。