www.gusucode.com > signal 工具箱matlab源码程序 > signal/snr.m
function [r, noisePow] = snr(varargin) %SNR Signal to Noise Ratio % R = SNR(X, Y) computes the signal to noise ratio (SNR) in dB, by % computing the ratio of the summed squared magnitude of the signal, X, % to the summed squared magnitude of the noise, Y, where Y has the same % dimensions as X. Use this form of SNR when your input signal is not % sinusoidal and you have an estimate of the noise. % % R = SNR(X) computes the signal to noise ratio (SNR), 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 % and excludes the power of first six harmonics (including the % fundamental). % % R = SNR(X, Fs, N) computes the signal to noise ratio (SNR) in dBc, of % the real sinusoidal input signal, X, with sampling rate, Fs, and number % of harmonics, N, to exclude from computation when computing SNR. The % default value of Fs is 1. The default value of N is 6 and includes the % fundamental frequency. % % R = SNR(Pxx, F, 'psd') specifies the input as a one-sided PSD estimate, % Pxx, of a real signal. F is a vector of frequencies that corresponds % to the vector of Pxx estimates. The computation of noise excludes the % first six harmonics (including the fundamental). % % R = SNR(Pxx, F, N, 'psd') specifies the number of harmonics, N, to % exclude when computing SNR. The default value of N is 6 and includes % the fundamental frequency. % % R = SNR(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 = SNR(Sxx, F, RBW, N, 'power') specifies the number of harmonics, N, % to exclude when computing SNR. The default value of N is 6 and % includes the fundamental frequency. % % [R, NOISEPOW] = SNR(...) also returns the total noise power of the non- % harmonic components of the signal. % % [...] = SNR(...,'aliased') also excludes harmonics of the fundamental % aliased into the Nyquist range. Use this option when the input signal % is undersampled. If unspecified or 'omitaliases' is used instead, % harmonics of the fundamental frequency beyond Nyquist are treated as % noise. This works for all signatures listed above except SNR(X, Y). % % SNR(...) with no output arguments plots the spectrum of the signal and % annotates the fundamental, DC component, harmonics and noise. The DC % component is removed before computing SNR. This works for all % signatures listed above except SNR(X, Y). % % % Example 1: % % Compute the SNR of a 2 second 20ms rectangular pulse sampled at % % 10 kHz in the presence of gaussian noise % % Tpulse = 20e-3; Fs = 10e3; % x = rectpuls((-1:1/Fs:1),Tpulse); % y = 0.00001*randn(size(x)); % s = x + y; % pulseSNR = snr(x,s-x) % % % Example 2: % % Plot the SNR of a 2.5 kHz distorted sinusoid sampled at 48 kHz % load('sineex.mat','x','Fs'); % snr(x,Fs) % % % Example 3: % % Generate the periodogram of a 2.5 kHz distorted sinusoid sampled % % at 48 kHz and measure the SNR (in dB) % load('sineex.mat','x','Fs'); % w = kaiser(numel(x),38); % [Sxx, F] = periodogram(x,w,numel(x),Fs,'power'); % % % Measure SNR on the power spectrum % rbw = enbw(w,Fs); % sineSNR = snr(Sxx,F,rbw,'power') % % % annotate the spectrum % snr(Sxx,F,rbw,'power') % % See also SINAD THD SFDR TOI. % Copyright 2013 The MathWorks, Inc. narginchk(1,5); % handle canonical definition of SNR if nargin == 2 && isequal(size(varargin{1}), size(varargin{2})) [r, noisePow] = sampleSNR(varargin{1}, varargin{2}); return end % look for psd or power window compensation flags [estType, varargin] = psdesttype({'psd','power','time'},'time',varargin); % allow undersampled harmonics when 'aliased' is specified [harmType, varargin] = getmutexclopt({'aliased','omitaliases'},'omitaliases',varargin); % check for unrecognized strings chknostropts(varargin{:}); % plot if no arguments are specified plotType = distplottype(nargout, estType); switch estType case 'psd' [r, noisePow] = psdSNR(plotType, harmType, varargin{:}); case 'power' [r, noisePow] = powerSNR(plotType, harmType, varargin{:}); case 'time' [r, noisePow] = timeSNR(plotType, harmType, varargin{:}); end %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [r, noisePow] = sampleSNR(x, y) signalPow = rssq(x(:)).^2; noisePow = rssq(y(:)).^2; r = 10 * log10(signalPow / noisePow); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [r, noisePow] = timeSNR(plotType, harmType, x, fs, nHarm) % force column vector before checking attributes if max(size(x)) == numel(x) x = x(:); end validateattributes(x,{'numeric'},{'real','finite','vector'}, ... 'snr','x',1); if nargin > 3 validateattributes(fs, {'numeric'},{'real','finite','scalar','positive'}, ... 'snr','Fs',2); else fs = 1; end if nargin > 4 validateattributes(nHarm,{'numeric'},{'integer','finite','positive','scalar','>',1}, ... 'snr','N',3); else nHarm = 6; 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'); % specify sample rate for aliased harmonics (assume same as Fs) if strcmp(harmType, 'aliased') aliasFs = fs; else aliasFs = NaN; end [r, noisePow] = computeSNR(plotType, Pxx, F, rbw, nHarm, aliasFs); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [r, noisePow] = powerSNR(plotType, harmType, Sxx, F, rbw, nHarm) if F(1)~=0 error(message('signal:snr:MustBeOneSidedSxx')); end % ensure specified RBW is larger than a bin width df = mean(diff(F)); validateattributes(rbw,{'double'},{'real','finite','positive','scalar','>=',df}, ... 'snr','RBW',3); if nargin > 5 validateattributes(nHarm,{'double'},{'integer','finite','positive','scalar','>',1}, ... 'snr','N',4); else nHarm = 6; end % assume frequency transform is from an even-length input when % undersampling harmonics (i.e. F(end) = Fs/2) if strcmp(harmType, 'aliased') aliasFs = 2*F(end); else aliasFs = NaN; end [r, noisePow] = computeSNR(plotType, Sxx/rbw, F, rbw, nHarm, aliasFs); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [r, noisePow] = psdSNR(plotType, harmType, Pxx, F, nHarm) if F(1)~=0 error(message('signal:snr:MustBeOneSidedPxx')); end % use the average bin width df = mean(diff(F)); if nargin > 4 validateattributes(nHarm,{'double'},{'integer','finite','positive','scalar','>',1}, ... 'snr','N',4); else nHarm = 6; end % assume frequency transform is from an even-length input when % undersampling harmonics (i.e. F(end) = Fs/2) if strcmp(harmType, 'aliased') aliasFs = 2*F(end); else aliasFs = NaN; end [r, noisePow] = computeSNR(plotType, Pxx, F, df, nHarm, aliasFs); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [r, noisePow] = computeSNR(plotType, Pxx, F, rbw, nHarm, aliasFs) % save a copy of the original PSD estimates origPxx = Pxx; % pre-allocate harmonic table psdHarmPow = NaN(nHarm,1); psdHarmFreq = NaN(nHarm,1); harmIdx = NaN(nHarm, 2); % 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, Ffund, iFund, iLeft, iRight] = signal.internal.getToneFromPSD(Pxx, F, rbw); [psdHarmPow(1), psdHarmFreq(1)] = idx2psddb(Pxx, F, iFund); Pxx(iLeft:iRight) = 0; harmIdx(1, :) = [iLeft; iRight]; % remove harmonic content for i=2:nHarm toneFreq = i*Ffund; if isfinite(aliasFs) toneFreq = aliasToNyquist(i*Ffund, aliasFs); end [harmPow, ~, iHarm, iLeft, iRight] = signal.internal.getToneFromPSD(Pxx, F, rbw, toneFreq); [psdHarmPow(i), psdHarmFreq(i)] = idx2psddb(Pxx, F, iHarm); % obtain local maximum value in neighborhood of bin if ~isnan(harmPow) % remove the power of this tone Pxx(iLeft:iRight) = 0; harmIdx(i, :) = [iLeft; iRight]; end end % 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 noise distortion. totalNoise = bandpower(Pxx, F, 'psd'); r = 10*log10(Pfund / totalNoise); noisePow = 10*log10(totalNoise); if ~strcmp(plotType,'none') plotSNR(origPxx, F, rbw, plotType, psdHarmFreq, psdHarmPow, dcIdx, harmIdx); title(getString(message('signal:snr:SNRResult',sprintf('%6.2f',r)))); end function f = aliasToNyquist(f, fs) f = mod(f,fs); if f > fs/2 f = fs-f; end function plotSNR(Pxx, F, rbw, plotType, psdHarmFreq, psdHarmPow, dcIdx, harmIdx) % scale Pxx by rbw Pxx = Pxx * rbw; psdHarmPow = psdHarmPow + 10*log10(rbw); % initialize distortion plot [hAxes, F, fscale, colors] = initdistplot(plotType, F); % --- plot legend entry items --- % plot fundamental xData = F(harmIdx(1,1):harmIdx(1,2)); yData = 10*log10(Pxx(harmIdx(1,1):harmIdx(1,2))); line(xData, yData, 'Parent', hAxes, 'Color', colors(1,:)); % plot noise 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 over the grid --- % plot fundamental marker xData = psdHarmFreq(1)*fscale; yData = psdHarmPow(1); text(xData(1),yData(1),'F', ... 'VerticalAlignment','bottom', ... 'HorizontalAlignment','center', ... 'BackgroundColor', 'w', ... 'EdgeColor', 'k', ... 'Color', colors(1,:)); % plot harmonic markers xData = psdHarmFreq(2:end)*fscale; yData = psdHarmPow(2:end); for i=1:numel(xData) text(xData(i),yData(i),num2str(i+1), ... 'VerticalAlignment','bottom', ... 'HorizontalAlignment','center', ... 'BackgroundColor', 'w', ... 'EdgeColor', 'k', ... 'Color', colors(3,:)); end % plot noise line xData = F; yData = 10*log10(Pxx); line(xData, yData, 'Parent', hAxes, 'Color', colors(2,:)); % plot fundamental xData = F(harmIdx(1,1):harmIdx(1,2)); yData = 10*log10(Pxx(harmIdx(1,1):harmIdx(1,2))); line(xData, yData, 'Parent', hAxes, 'Color', colors(1,:)); % plot dc on top xData = F(dcIdx(1):dcIdx(2)); yData = 10*log10(Pxx(dcIdx(1):dcIdx(2))); line(xData, yData, 'Parent', hAxes,'Color', colors(3,:)); % plot harmonics for i=2:size(harmIdx,1) if ~any(isnan(harmIdx(i,:))) xData = F(harmIdx(i,1):harmIdx(i,2)); yData = 10*log10(Pxx(harmIdx(i,1):harmIdx(i,2))); line(xData, yData, 'Parent', hAxes,'Color', colors(3,:)); end end legend(getString(message('signal:snr:Fundamental')), ... getString(message('signal:snr:Noise')), ... getString(message('signal:snr:DCHarmonics')), ... 'Location','best');