www.gusucode.com > signal 工具箱matlab源码程序 > signal/private/rpmmap.m
function [varargout] = rpmmap(x,fs,rpm,maptype,varargin) %RPMMAP Compute frequency or order maps vs. rpm values. % Copyright 2015 MathWorks, Inc. % parse input parameters [resolution,win,winparam,amplitude,scale,overlapPercent] = ... parseOptions(maptype,x,fs,rpm,varargin{:}); % cast to enforce precision rules (we already checked that the inputs are % numeric. rpm = double(rpm); fs = double(fs); resolution = double(resolution); overlapPercent = double(overlapPercent); % Define time vector time = (0:length(x)-1)/fs; % Convert x and rpm to column vectors x = x(:); rpm = rpm(:); time = time(:); if strcmp(maptype,'order') % Convert signal to rotational or order domain. % xp is the resampled signal (constant samples/cycle) and fsp is the % samples/cycle rate. phaseUp, rpmUp and timeUp vectors are at an % upsampled rate of 15*fs. [xp, fsp, phaseUp, rpmUp, timeUp] = toConstantSamplesPerCycle(x,fs,rpm,time); else xp = x; fsp = fs; end % Compute DFT window length % Define the minumum and maximum allowed window lengths minWindowLength = 4; maxWindowLength = length(xp); if isempty(resolution) % Use the default resolution value: (sampling frequency)/128 for % frequency maps and (sampling frequency)/256 for order maps. if strcmp(maptype,'order') resolution = fsp/256; else resolution = fsp/128; end % Make sure resolution is inside range of allowed values % If it is not, use the minimum or maximum allowed window length [resolution,winLengthBound] = ... validateResolutionValue(resolution,minWindowLength,maxWindowLength,fsp,win,winparam,false); else % Validate input resolution value [resolution,winLengthBound] = ... validateResolutionValue(resolution,minWindowLength,maxWindowLength,fsp,win,winparam,true); end if isempty(winLengthBound) % We have a valid resolution so compute window length [~,winLength] = signal.internal.getWinDurationForAGivenRBW(... resolution,win,winparam,fsp,true); else % We have an invalid resolution so use the upper or lower bound winLength = winLengthBound; end % Compute the number of overlapping samples nOverlap = min(ceil(overlapPercent/100*winLength),winLength-1); % Create DFT window and normalize it winfun = getWinFunction(win, winparam); win = winfun(winLength); win = win/sum(win); % Compute the resolution of the window and the FFT length winres = enbw(win)*fsp/winLength; nfft = max(256,length(win)); % Generate the frequency/order map if strcmp(maptype,'order') [map,ordSpec,phaseSpec] = spectrogram(xp,win,nOverlap,nfft,fsp,'onesided'); phaseSpec = phaseSpec(:); else [map,freqSpec,timeSpec] = spectrogram(xp,win,nOverlap,nfft,fsp,'onesided'); timeSpec = timeSpec(:); end % Apply amplitude and scale of the map map = mapAmplitudeScale(map,amplitude,scale,nfft); if strcmp(maptype,'order') % Limit to max order value as we oversampled in the order domain for % better results. Recall that Omax = fs/(2*max(rpm/60)) idx = find(ordSpec <= fs/(2*max(rpm/60))); if idx(end)+1 <= length(ordSpec) % Get one extra order point to ensure we have the max order value idx(end+1) = idx(end)+1; end ordSpec = ordSpec(idx); map = map(idx,:); end %--------------------------------------------------------------------- % Plot or assign outputs %--------------------------------------------------------------------- % Create plot of ordermap if no output arguments are specified % Otherwise, assign output variables to suppress command window output if nargout == 0 % Surf encodes data points at the edges and takes the color from the last % edge so we need to add an additional point so that surf does the right % thing. This is important especially when spectrogram has only one % estimate (e.g. window length = signal length). For the plot we set time % or phase values to be at: nwin/2-a/2, nwin/2+a/2, nwin/2+3a/2, % nwin/2+5a/2 ... where a is the number of new samples for each segment % (i.e. nwin-noverlap). For the case of zero overlap this corresponds to % 0, nwin, 2*nwin, ... a = winLength - nOverlap; if strcmpi(maptype,'order') phaseVectForPlot = [(winLength/2-a/2)/fsp; phaseSpec + ((a/2)/fsp)]; rpmVectForPlot = interp1(phaseUp, rpmUp, phaseVectForPlot,'linear','extrap'); timeVectForPlot = interp1(phaseUp, timeUp, phaseVectForPlot,'linear','extrap'); map = [map map(:,end)]; rpmmapplot(map, ordSpec, timeVectForPlot, rpmVectForPlot, ... winres,'order', scale, amplitude, rpm, time); else timeVectForPlot = [(winLength/2-a/2)/fsp; timeSpec + ((a/2)/fsp)]; rpmVectForPlot = interp1(time, rpm, timeVectForPlot,'linear','extrap'); map = [map map(:,end)]; rpmmapplot(map, freqSpec, timeVectForPlot, rpmVectForPlot, ... winres,'frequency', scale, amplitude, rpm, time); end else % Assign output variables varargout{1} = map; if strcmp(maptype,'order') % Order map outputs are [map, order, rpm, time, res] if nargout > 1 varargout{2} = ordSpec; if nargout > 2 % Interpolate RPM values based on upsampled phase and rpm values, % and on spectrogram phase output vector rpmSpecDerived = interp1(phaseUp,rpmUp,phaseSpec,'linear','extrap'); varargout{3} = rpmSpecDerived; if nargout > 3 % Interpolate time values based on upsampled phase and rpm values, % and on spectrogram phase output vector timeSpecDerived = interp1(phaseUp,timeUp,phaseSpec,'linear','extrap'); varargout{4} = timeSpecDerived; if nargout > 4 varargout{5} = winres; end end end end else if nargout > 1 varargout{2} = freqSpec; if nargout > 2 % Interpolate RPM values based on time and rpm values, and on % spectrogram time output vector rpmSpecDerived = interp1(time,rpm,timeSpec,'linear','extrap'); varargout{3} = rpmSpecDerived; if nargout > 3 varargout{4} = timeSpec; if nargout > 4 varargout{5} = winres; end end end end end end end %-------------------------------------------------------------------------- function [xp, fsp, phaseUp, rpmUp, timeUp] = toConstantSamplesPerCycle(x,fs,rpm,time) % Compute the maximum order that can be present in X with no aliasing. Max % frequency, fmax, is the signal's max frequency (max(rpm/60)) multiplied % by max order. For Nyquist to hold fs > 2*fmax = 2*max(rpm/60)*Omax. Omax = fs/(2*max(rpm/60)); % Define sampling rate (samples/cycle) in phase domain based on the maximum % order. Make sampling rate 4 times the Nyquist rate in the order domain. fsp = 4*(2*Omax); % Define upsample factor. Upsampling will improve accuracy when converting % from constant samples/second to constant samples/cycle. In the worst case % scenario, when time signal is critically sampled in time at Fmax/2, we % are increasing the Nyquist frequency by 15. upFactor = 15; % Upsample x and rpm if isa(x,'single') xUp = resample(double(x),upFactor,1); xUp = single(xUp); else xUp = resample(x,upFactor,1); end % Get upsampled time and rpm vectors timeUp = (0:length(xUp)-1).'/(upFactor*fs); rpmUp = interp1(time, rpm, timeUp, 'linear','extrap'); % Estimate the phase of each signal sample by integrating the instantaneous % signal frequency which is rpmUp/60. Divide by sampling rate which is % upFactor*fs; phaseUp = cumtrapz(rpmUp/(60*upFactor*fs)); % Interpolate signal x at constant phase increments (i.e. constant % samples/cycle). xp is uniformly sampled in the rotational domain --> same % samples per rotation for any rpm constPhase = (phaseUp(1):1/fsp:(phaseUp(end)))'; xp = interp1(phaseUp,xUp,constPhase,'linear','extrap'); end %-------------------------------------------------------------------------- function [resolution, win, winparam, amplitude, scale, overlapPercent, funName] = parseOptions(maptype,x,fs,rpm,varargin) % default values for n-v pairs if strcmp(maptype,'order') win = 'flattopwin'; funName = 'rpmordermap'; else win = 'hann'; funName = 'rpmfreqmap'; end winparam = []; amplitude = 'rms'; scale = 'linear'; resolution = []; overlapPercent = 50; % Check valid values for x,fs,rpm validateattributes(x,{'single','double'},... {'real','nonsparse','finite','nonnan','vector'},funName,'X'); validateattributes(fs,{'numeric'},... {'real','positive','nonsparse','finite','nonnan','scalar'},funName,'Fs'); validateattributes(rpm,{'numeric'},... {'real','positive', 'nonsparse','finite','nonnan','vector'},funName,'RPM'); % Check valid input dimensions if length(x) < 18 error(message('signal:rpmmap:MustBeMinLength')); end %check rpm is the same size and x if length(x)~=length(rpm) error(message('signal:rpmmap:MustBeSameLength')); end %check if a resolution parameter is specified and check if valid value if ~isempty(varargin) && isnumeric(varargin{1}) resolution = varargin{1}; validateattributes(resolution,{'numeric'},... {'real','positive','nonsparse','finite','nonnan','scalar'},funName,'resolution'); varargin(1) = []; end % If all is valid, then parse name value pairs if we have any if isempty(varargin) return; end %check that name-value inputs come in pairs and are all strings if rem(numel(varargin),2) error(message('signal:rpmmap:NVMustBeEven')); end idx =1; while idx < numel(varargin) name = validatestring(varargin{idx},{'Window','Amplitude','Scale','OverlapPercent'}); value = varargin{idx+1}; switch name case 'Window' winList = {'hann','hamming','flattopwin','kaiser','rectwin','chebwin'}; if iscell(value) win = validatestring(value{1}, winList); winparam = []; if numel(value) == 2 switch win case 'kaiser' winparam = value{2}; validateattributes(winparam,{'numeric'},... {'real', 'positive', 'finite','nonnan','scalar'},... funName,'Kaiser beta parameter'); case 'chebwin' winparam = value{2}; validateattributes(winparam,{'numeric'},... {'real', 'positive', 'finite','nonnan','scalar','>=',45},... funName,'chebwin sidelobe attenuation parameter'); otherwise error(message('signal:rpmmap:WindowCellArgumentInvalid',win)); end end else win = validatestring(value,winList); end case 'Amplitude' amplitude = validatestring(value,{'rms','peak','power'}); case 'Scale' scale = validatestring(value,{'linear','dB'}); case 'OverlapPercent' overlapPercent = value; validateattributes(overlapPercent,{'numeric'},... {'scalar','nonnegative','<=',100},funName,'OP'); end idx = idx + 2; end end %-------------------------------------------------------------------------- function [resolution,winLength] = validateResolutionValue(resolution,minLength,maxLength,fs,win,winParam,warn) %check that resolution falls between minimum and maximum allowable values %either return a valid resolution or a window length % Create window winfun = getWinFunction(win, winParam); % Compute resolution of the window minResValue = enbw(winfun(maxLength))*fs/maxLength; maxResValue = enbw(winfun(minLength))*fs/minLength; sminres = sprintf('%0.5g',minResValue); smaxres = sprintf('%0.5g',maxResValue); if resolution <= minResValue winLength = maxLength; resolution = []; if(warn) warning(message('signal:rpmmap:ResolutionMustBeInRangeUnder',... sminres,smaxres)); end elseif resolution >= maxResValue winLength = minLength; resolution = []; if(warn) warning(message('signal:rpmmap:ResolutionMustBeInRangeOver',... sminres,smaxres)); end else winLength = []; end end %-------------------------------------------------------------------------- function map = mapAmplitudeScale(map,amplitude,scale,nfft) % Map amplitude switch amplitude case 'peak' map = oneSidedSpectrum(abs(map),nfft); % Scale all components except dc map(2:end,:) = sqrt(2)*map(2:end,:); case 'rms' map = oneSidedSpectrum(abs(map),nfft); case 'power' map = oneSidedSpectrum(abs(map),nfft).^2; end % Map scale switch scale case 'linear' %no change is needed case 'dB' if isequal(amplitude,'power') map = 10*log10(map); else map = 20*log10(map); end end end %-------------------------------------------------------------------------- function mp = oneSidedSpectrum(map,nfft) % Get the correctly scaled one sided spectrum - map has magnitude values % not power values so scaling should be sqrt(2) instead of 2 (note that we % already have half the spectrum but we need to scale it correctly). if rem( nfft, 2 ) %fft is odd length % Don't double DC component mp = [ map( 1, : ); sqrt(2)*map( 2:end, : ) ]; else % Don't double DC component or unique Nyquist point mp = [ map( 1, : ); sqrt(2)*map( 2:end - 1, : ); map( end, : ) ]; end end %-------------------------------------------------------------------------- function winfun = getWinFunction(win,winParam) if strcmp(win,'kaiser') || strcmp(win,'chebwin') winname = str2func(win); winfun = @(x) winname(x,winParam); else winfun = str2func(win); end end