www.gusucode.com > signal 工具箱matlab源码程序 > signal/sinad.m
function [r, totDistPow] = sinad(varargin) %SINAD Signal to Noise and Distortion ratio % R = SINAD(X) computes the signal to noise and distortion ratio (SINAD), % in dBc, of the real sinusoidal input signal X. The computation is % performed over a periodogram of the same length as the input using a % Kaiser window. % % R = SINAD(X, Fs) specifies the sampling rate, Fs. The % default value of Fs is 1. % % R = SINAD(Pxx, F, 'psd') specifies the input as a one-sided % PSD estimate, Pxx, of a real signal. F is a vector of frequencies % that correspond to the PSD estimates. % % R = SINAD(Sxx, F, RBW, 'power') specifies the input as a one-sided % power spectrum, Sxx, of a real signal. RBW is the resolution bandwidth % over which each power estimate is integrated. % % [R, TOTDISTPOW] = SINAD(...) also returns the total noise and harmonic % distortion power of the signal. % % SINAD(...) with no output arguments plots the spectrum of the signal % and annotates the fundamental signal and noise. The DC component is % removed before computing SINAD. % % % Example 1: % % Plot the SINAD ratio of a 2.5 kHz distorted sinusoid sampled % % at 48 kHz % load('sineex.mat','x','Fs'); % sinad(x,Fs) % % % Example 2: % % Generate the periodogram of a 2.5 kHz distorted sinusoid sampled % % at 48 kHz and compute SINAD % % % create the sinusoid and compute the power spectrum % load('sineex.mat','x','Fs'); % w = kaiser(numel(x),38); % [Sxx, F] = periodogram(x,w,numel(x),Fs,'power'); % % % compute SINAD % rbw = enbw(w,Fs); % [sineSINAD, distPower] = sinad(Sxx, F, rbw, 'power') % % % plot and annotate the spectrum % sinad(Sxx, F, rbw, 'power') % % See also THD SFDR SNR TOI. % Copyright 2013 The MathWorks, Inc. narginchk(1,4); % look for psd or power window compensation flags [esttype, varargin] = psdesttype({'psd','power','time'},'time',varargin); % check for unrecognized strings chknostropts(varargin{:}); % plot if no arguments are specified plotType = distplottype(nargout, esttype); switch esttype case 'psd' [r, totDistPow] = psdSINAD(plotType, varargin{:}); case 'power' [r, totDistPow] = powerSINAD(plotType, varargin{:}); case 'time' [r, totDistPow] = timeSINAD(plotType, varargin{:}); end %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [r, totDistPow] = timeSINAD(plotType, x, fs) % force column vector before checking attributes if max(size(x)) == numel(x) x = x(:); end validateattributes(x,{'numeric'},{'real','finite','vector'}, ... 'sinad','x',1); if nargin > 2 validateattributes(fs, {'numeric'},{'real','finite','scalar','positive'}, ... 'sinad','Fs',2); else fs = 1; end % remove DC component x = x - mean(x); n = length(x); % use Kaiser window to reduce effects of leakage w = kaiser(n,38); rbw = enbw(w,fs); [Pxx, F] = periodogram(x,w,n,fs,'psd'); [r, totDistPow] = computeSINAD(plotType, Pxx, F, rbw); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [r, totDistPow] = powerSINAD(plotType, Sxx, F, rbw) if F(1)~=0 error(message('signal:sinad:MustBeOneSidedSxx')); end % use the average bin width df = mean(diff(F)); validateattributes(rbw,{'numeric'},{'real','finite','positive','scalar','>=',df}, ... 'sinad','RBW',3); [r, totDistPow] = computeSINAD(plotType, Sxx/rbw, F, rbw); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [r, totDistPow] = psdSINAD(plotType, Pxx, F) if F(1)~=0 error(message('signal:sinad:MustBeOneSidedPxx')); end % ensure specified RBW is larger than a bin width df = mean(diff(F)); [r, totDistPow] = computeSINAD(plotType, Pxx, F, df); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [r, noisePow] = computeSINAD(plotType, Pxx, F, rbw) % save a copy of the original PSD estimates origPxx = Pxx; % bump DC component by 3dB and remove it. Pxx(1) = 2*Pxx(1); [~, ~, ~, iLeft, iRight] = signal.internal.getToneFromPSD(Pxx, F, rbw, 0); Pxx(iLeft:iRight) = 0; dcIdx = [iLeft; iRight]; % get an estimate of the actual frequency / amplitude, then remove it. [Pfund, ~, iFund, iLeft, iRight] = signal.internal.getToneFromPSD(Pxx, F, rbw); Pxx(iLeft:iRight) = 0; fundIdx = [iLeft; iRight]; % get an estimate of the noise floor by computing the median % noise power of the non-harmonic region estimatedNoiseDensity = median(Pxx(Pxx>0)); % extrapolate estimated noise density into dc/signal/harmonic regions Pxx(Pxx==0) = estimatedNoiseDensity; % prevent estimate from obscuring low peaks Pxx = min([Pxx origPxx],[],2); % compute the remaining harmonic and noise distortion. totalNoise = bandpower(Pxx, F, 'psd'); r = 10*log10(Pfund / totalNoise); noisePow = 10*log10(totalNoise); if ~strcmp(plotType,'none') plotSINAD(origPxx, F, rbw, plotType, iFund, dcIdx, fundIdx); title(getString(message('signal:sinad:SINADResult',sprintf('%6.2f',r)))); end function plotSINAD(Pxx, F, rbw, plotType, iFund, dcIdx, fundIdx) % scale Pxx by rbw Pxx = Pxx * rbw; % initialize distortion plot [hAxes, F, ~, colors] = initdistplot(plotType, F); % --- plot legend entries --- % plot fundamental line xData = F(fundIdx(1):fundIdx(2)); yData = 10*log10(Pxx(fundIdx(1):fundIdx(2))); line(xData, yData, 'Parent', hAxes, 'Color', colors(1,:)); % plot noise and distortion line xData = F; yData = 10*log10(Pxx); line(xData, yData, 'Parent', hAxes, 'Color', colors(2,:)); % plot DC xData = F(dcIdx(1):dcIdx(2)); yData = 10*log10(Pxx(dcIdx(1):dcIdx(2))); line(xData, yData, 'Parent', hAxes, 'Color', colors(3,:)); % --- use a solid grid slightly offset to accommodate text labels --- initdistgrid(hAxes); % --- replot on top of grid --- % plot noise and distortion line xData = F; yData = 10*log10(Pxx); line(xData, yData, 'Parent', hAxes, 'Color', colors(2,:)); % plot DC xData = F(dcIdx(1):dcIdx(2)); yData = 10*log10(Pxx(dcIdx(1):dcIdx(2))); line(xData, yData, 'Parent', hAxes, 'Color', colors(3,:)); % plot fundamental marker xData = F(iFund); yData = 10*log10(Pxx(iFund)); text(xData(1),yData(1),'F', ... 'VerticalAlignment','bottom', ... 'HorizontalAlignment','center', ... 'BackgroundColor','w', ... 'EdgeColor','k', ... 'Color', colors(1,:)); % plot fundamental line xData = F(fundIdx(1):fundIdx(2)); yData = 10*log10(Pxx(fundIdx(1):fundIdx(2))); line(xData, yData, 'Parent', hAxes, 'Color', colors(1,:)); legend(getString(message('signal:sinad:Fundamental')), ... getString(message('signal:sinad:NoiseAndDistortion')), ... getString(message('signal:sinad:DC')));