f_sndr.m

在通讯系统中

% function [din_p din_f] = f_sndr(din, fs, winname, fbw, fin, mslope, nodc);
%
% Description : Get the frequency domain information of signal
% 
% Inputs    : 
%   din     : Input data time domain
%   fs      : Data sample rate, default 2, that fs/2=1
%   winname : Window name, default none, else hanning ....
%   fbw     : Bband width of interest
%   fin     : Input data frequency guide setting
%   mslope  : Frequency signal main slope, default 10
%   nodc    : Get rid of DC, default 0, keep DC
%
% Outputs   : 
%   din_p   : Power spectrum of frequency
%   din_f   : Frequency

function [din_p din_f] = f_sndr(din, fs, winname, fbw, fin, mslope, nodc);

len_din = length(din);
din = din(:);

if (nargin <7) | isempty(nodc)
    nodc = 0;
end

if (nargin < 6) | isempty(mslope)
    mslope = 5;
end

if (nargin < 5) | isempty(fin)
    fin = 0; % Input signal frequency
end

if (nargin < 4) | isempty(fbw)
    if nargin < 2
        fbw = 1;
    else
        fbw = fs/2; % Signal band width around Input signal frequency
    end
end

if (nargin < 3) | isempty(winname)
    dwin = ones(1, len_din)';
else
     dwin = window(winname, len_din);
    end

if (nargin < 2) | isempty(fs)
    fs = 2; % Signal band width around Input signal frequency
end

len_din_half = round(len_din/2);
din = din.*dwin;
din_p_temp = (abs(fft(din/len_din_half))).^2; % normalize power, 0db corresponding to full scale [-1 +1] sin wave
din_f_temp = fs/len_din*(0:len_din-1);

if (fin - fbw <0)
    din_bw_low = 1;
else
    din_bw_low = round(len_din/fs * (fin-fbw));
end

if (fin + fbw > fs/2)
    din_bw_high = len_din_half;
else
    din_bw_high = round(len_din/fs * (fin+fbw));
end

din_p = din_p_temp(din_bw_low:din_bw_high);
if (nodc)
    din_p(1:5) = 0; % get rid of dc
end
din_f = din_f_temp(din_bw_low:din_bw_high);
len_bw = length(din_p);

fin_ind = round(len_din/fs*fin);

% SNDR
if (fin~=0)
    [y i] = max(din_p(fin_ind-1:fin_ind+1));
else
    [y i] = max(din_p);
end

if (i-mslope<=0)
    low_ind = 1;
else
    low_ind = i-mslope;
end
if (i+mslope>=length(din_p))
    high_ind = length(din_p);
else
    high_ind = i+mslope;
end
s = sum(din_p(low_ind:high_ind));

nd = sum(din_p) - s;
sndr = 10*log10(s/nd);
disp(strcat('Max Signal frequency is: ', num2str(din_f(i)), 'Hz')); 
disp(strcat('SNDR is: ', num2str(sndr), 'db'));
disp(strcat('Frequency resolution is: ', num2str(fs/len_din), 'Hz'));

% THD and SNR
thd_p = 0;
for k=i*2:i:din_bw_high-10
    if (k-mslope<=0)
        low_ind = 1;
    else
        low_ind = k-mslope;
    end
    if (k+mslope>=length(din_p))
        high_ind = length(din_p);
    else
        high_ind = k+mslope;
    end        
    thd_p = thd_p + sum(din_p(low_ind:high_ind));
end

n = nd - thd_p;
if( thd_p ~=0)
    thd = 10*log10(thd_p/s);
else
    thd = -200;
end

snr = 10*log10(s/n);
disp(strcat('THD is: ', num2str(thd), 'db'));
disp(strcat('SNR is: ', num2str(snr), 'db'));

pdb = zeros(1, len_bw);
for k=1:len_bw
    if (din_p(k) < 1e-20)
        pdb(k) = -200;
    else
        pdb(k) = 10*log10(din_p(k));
    end
end
figure;
plot(din_f, pdb);
xlabel('Frequency [Hz]');
ylabel('Amplitude [Db]');