snr2.m

自己写的信噪比计算工具snr2.m snr2(s2ana,fs,f0,fl,fh,bl,bh) s2ana: 待分析信号 fs:采样频率 f0:信号频率 [fl,fh]:信号能量分散带

%自己写的信噪比计算工具snr2.m 
%snr2(s2ana,fs,f0,fl,fh,bl,bh)
%s2ana: 待分析信号
%fs:采样频率
%f0:信号频率
%[fl,fh]:信号能量分散带


function snr = snr2(s2ana,fs,f0,fl,fh,bl,bh);
[sizea,sizeb] = size(s2ana);
if (sizea==1)
    s = s2ana';
else
    s = s2ana;
end
len = length(s);

s_window = s.*hann(len)/max(s2ana);
fft_s = fft(s_window,len)/len;
psd_s = fft_s .*conj(fft_s);
len1 = fix(len/2);
psd1_s = psd_s(1:len1)*65/2;
figure;
plot((0:len1-1)/len*fs,10*log10(psd1_s));grid;

%compute the fft line corresponding to fl,fh,bl,bh;
fp_fl = fix(fl*len/fs)+1;
if (fp_fl==0)
    fp_fl = 1;
end

fp_fh = fix(fh*len/fs)+1;
if (fp_fh>len1)
    fp_fh = len1;
end

fp_bl = fix(bl*len/fs)+1;
if (fp_bl==0),
    fp_bl = 1;
end

fp_bh = fix(bh*len/fs)+1;
if (fp_bh>len1),
    fp_bh=len1;
end

energy_sig = 0;
energy_noise = 0;
energy_sig = sum(psd1_s(fp_fl:fp_fh));
energy_noise = sum(psd1_s(fp_bl:fp_bh)) - energy_sig ...
    + (fp_fh-fp_fl+1)*(psd1_s(fp_fl)+psd1_s(fp_fh))/2;

snr = 10*log10(energy_sig /  energy_noise);