www.gusucode.com > signal 工具箱matlab源码程序 > signal/toi.m
function [oip3, fundPow, fundFreq, imodPow, imodFreq] = toi(varargin) %TOI Third Order Intercept point % OIP3 = TOI(X) computes the output third order intercept point, in dB, % of the real sinusoidal two-tone input signal X. The computation is % performed over a periodogram of the same length as the input using a % Kaiser window. % % OIP3 = TOI(X, Fs) specifies the sampling rate, Fs. The default value % of Fs is 1. % % OIP3 = TOI(Pxx, F, 'psd') specifies the input as a one-sided PSD, Pxx, % of a real signal. F is a vector of frequencies that corresponds to the % vector of Pxx estimates. % % OIP3 = TOI(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. % % [OIP3, FUNDPOW, FUNDFREQ, IMODPOW, IMODFREQ] = TOI(...) also returns the % power, FUNDPOW, and frequencies, FUNDFREQ of the two fundamental % sinusoids, and the power IMODPOW, and frequencies, IMODFREQ, of the % lower and upper intermodulation products. FUNDPOW, FUNDFREQ, IMODPOW, % and IMODFREQ are a two-element row vector containing the lower and upper % fundamentals or intermodulation products, respectively. % % TOI(...) with no output arguments plots the spectrum of the signal and % annotates the lower and upper fundamentals (F1, F2) and intermodulation % products (2 F1 - F2, 2 F2 - F1). % % % Example 1: % % create a 5 kHz and 6kHz sinusoid sampled at 48 kHz % Fi1 = 5e3; Fi2 = 6e3; Fs = 48e3; N = 1000; % x = sin(2*pi*Fi1/Fs*(1:N)) + sin(2*pi*Fi2/Fs*(1:N)); % % % model nonlinearity via polynomial and some noise % y = polyval([.5e-3 1e-7 .1 3e-3], x) + 0.00001*randn(1,N); % % % plot and annotate the spectrum % toi(y, Fs); % % % Example 2: % % create a 5 kHz and 6kHz sinusoid sampled at 48 kHz % Fi1 = 5e3; Fi2 = 6e3; Fs = 48e3; N = 1000; % x = sin(2*pi*Fi1/Fs*(1:N)) + sin(2*pi*Fi2/Fs*(1:N)); % % % model nonlinearity via polynomial and some noise % y = polyval([.5e-3 1e-7 .1 3e-3], x) + 0.00001*randn(1,N); % % % use a Kaiser window % w = kaiser(numel(y),38); % % % compute via a power spectrum % [Sxx, F] = periodogram(y,w,N,Fs,'power'); % [myTOI, Pfund, Ffund, Pim3, Fim3] = toi(Sxx, F, enbw(w,Fs), 'power'); % % % plot and annotate the spectrum % toi(Sxx, F, enbw(w,Fs), 'power'); % % See also THD SFDR SINAD SNR. % 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' [oip3, fundPow, fundFreq, imodPow, imodFreq] = psdTOI(plotType, varargin{:}); case 'power' [oip3, fundPow, fundFreq, imodPow, imodFreq] = powerTOI(plotType, varargin{:}); case 'time' [oip3, fundPow, fundFreq, imodPow, imodFreq] = timeTOI(plotType, varargin{:}); end %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [oip3, fundPow, fundFreq, imodPow, imodFreq] = timeTOI(plotType, x, fs) % force column vector before checking attributes if max(size(x)) == numel(x) x = x(:); end validateattributes(x,{'numeric'},{'real','finite','vector'}, ... 'toi','x',1); if nargin > 2 validateattributes(fs, {'numeric'},{'real','finite','scalar','positive'}, ... 'toi','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); [Pxx, F] = periodogram(x,w,n,fs,'psd'); [oip3, fundPow, fundFreq, imodPow, imodFreq] = computeTOI(plotType, Pxx, F, enbw(w,fs)); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [oip3, fundPow, fundFreq, imodPow, imodFreq] = powerTOI(plotType, Sxx, F, rbw) if F(1)~=0 error(message('signal:toi:MustBeOneSidedSxx')); end % ensure specified RBW is larger than a bin width df = mean(diff(F)); validateattributes(rbw,{'numeric'},{'real','finite','positive','scalar','>=',df}, ... 'toi','RBW',3); [oip3, fundPow, fundFreq, imodPow, imodFreq] = computeTOI(plotType, Sxx/rbw, F, rbw); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [oip3, fundPow, fundFreq, imodPow, imodFreq] = psdTOI(plotType, Pxx, F) if F(1)~=0 error(message('signal:toi:MustBeOneSidedPxx')); end % ensure specified RBW is larger than a bin width df = mean(diff(F)); [oip3, fundPow, fundFreq, imodPow, imodFreq] = computeTOI(plotType, Pxx, F, df); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function [oip3, fundPow, fundFreq, imodPow, imodFreq] = computeTOI(plotType, Pxx, F, rbw) % plot one-sided power spectrum if ~strcmp(plotType,'none') [~, fscale] = psdplot(Pxx, F, rbw, plotType); end % bump DC component by 3dB and remove it. Pxx(1) = 2*Pxx(1); [~, ~, ~, iDCLeft, iDCRight] = signal.internal.getToneFromPSD(Pxx, F, rbw, 0); Pxx(iDCLeft:iDCRight) = 0; % get an estimate of the dominant frequency / amplitude, then temporarily remove it. [fundPow(1), fundFreq(1), idxTone, iLeft, iRight] = signal.internal.getToneFromPSD(Pxx, F, rbw); [psdFundPow(1), psdFundFreq(1)] = idx2psddb(Pxx, F, idxTone); tmpPxx = Pxx(iLeft:iRight); Pxx(iLeft:iRight) = 0; % get an estimate of the next dominant frequency / amplitude. [fundPow(2), fundFreq(2), idxTone] = signal.internal.getToneFromPSD(Pxx, F, rbw); [psdFundPow(2), psdFundFreq(2)] = idx2psddb(Pxx, F, idxTone); % restore signal (in case it interferes with the IMOD product) Pxx(iLeft:iRight) = tmpPxx; % ensure f2 > f1 if fundFreq(2) < fundFreq(1) fundFreq([1 2]) = fundFreq([2 1]); fundPow([1 2]) = fundPow([2 1]); psdFundPow([1 2]) = psdFundPow([2 1]); psdFundFreq([1 2]) = psdFundFreq([2 1]); end % get an estimate of the lower third imod product [imodPow(1), imodFreq(1), idxTone] = signal.internal.getToneFromPSD(Pxx, F, rbw, 2*fundFreq(1)-fundFreq(2)); [psdImodPow(1), psdImodFreq(1)] = idx2psddb(Pxx, F, idxTone); % get an estimate of the upper third imod product [imodPow(2), imodFreq(2), idxTone] = signal.internal.getToneFromPSD(Pxx, F, rbw, 2*fundFreq(2)-fundFreq(1)); [psdImodPow(2), psdImodFreq(2)] = idx2psddb(Pxx, F, idxTone); % If the second primary tone is the 2nd harmonic of the first primary tone, % then the lower intermodulation product is at DC. We return NaN under this % condition if isnan(imodFreq(1)) || imodFreq(1) <= F(iDCRight) imodPow(1) = NaN; imodFreq(1) = NaN; end % convert to dB fundPow = 10*log10(fundPow); imodPow = 10*log10(imodPow); % report back the output ip3 carrierPow = mean(fundPow); suppressPow = mean(fundPow) - mean(imodPow); oip3 = carrierPow + suppressPow/2; if ~strcmp(plotType,'none') if strcmp(plotType,'power') psdplottoipeaks(fscale, fundPow, fundFreq, imodPow, imodFreq); else % 'psd' psdplottoipeaks(fscale, psdFundPow, psdFundFreq, psdImodPow, psdImodFreq); end title(getString(message('signal:toi:TOIResult',sprintf('%6.2f',oip3)))); end