nspplot.m

一种新的时频分析方法的matlab源程序。

function [spec,tscale,fscale] = nspplot(freq,amp,twin,fwin);
% COMPONENTS/NSPPLOT generates and optionally plots a nonlinear spectrum
% using precalculated frequency and amplitude
%
% Usage:
%  [spec,tscale,fscale]=nsp(freq,amp[,twin][,fwin])
%
% Input:
%  freq       - the Components data structure containing the instantaneous
%               frequency of each component
%  amp        - the Components data structure containing the instantaneous
%               amplitude of each component
%  twin       - Time window vector: [start, end, res]
%               If unspecified or [], NSPPLOT defaults to spanning data with 400 "buckets"
%               If scalar, twin is taken to be 'res', the number of
%               "buckets" in time, and 'start' and 'end' span the signal.
%               If a two-element vector, 'res' defaults to 400.
%  fwin       - Frequency window vector [lower, upper, res]
%               If unspecified or [], NSPPLOT defaults to showing all frequencies
%               present, sampled into 400 "buckets".
%               If scalar, fwin is taken to be 'res', the number of
%               frequency "buckets," and all frequencies are shown.
%               If a two-element vector, 'res' defaults to 400.
% Output: (all outputs optional)
%  spec       - A 2D matrix containing the spectrum. Time increases with
%               increasing column number and frequency increases with
%               increasing row number
%  tscale     - A row vector giving the time values corresponding to the
%               x-axis indices 
%  fscale     - A column vector giving the frequency values corresponding
%               to the y-axis indices
%
% Request no outputs and the spectrum will be plotted in the current
% figure.
%
% Examples:
%
% * Get a quick look at the spectrum (normalized Hilbert):
%  >> nspplot(freq,amp);
%
% * Zoom in to look at time 0-1 s and frequency 10-20 Hz:
%  >> nspplot(freq,amp,[0,1],[10,20]);
%
% * Get the spectrum data to plot later:
%  >> [spec,tscale,fscale] = nspplot(freq,amp);
%  >> img(tscale,fscale,spec);

% Kenneth C. Arnold (for NASA GSFC), 2004-08-06

if nargin < 2; error('MATLAB:minrhs', 'Not enough arguments to NSPPLOT'); end
if nargin < 3; twin=[]; end
if nargin < 4; fwin=[]; end

% Default values (FIXME: make these exactly equal default MATLAB figure
% sizes.)
DFLT_TRES=400;
DFLT_FRES=400;

%----- Get frequency bounds
nfwin = max(size(fwin));
if nfwin <= 1
    minf = min(freq);
    if minf < 0
        warning('hht:negfreq', 'nsp: negative frequency encountered, setting minimum frequency to 0');
        minf = 0;
    end
    if nfwin == 1
        fres = fwin;
    else
        fres = DFLT_FRES; % Default frequency resolution
    end
    fwin = [minf max(freq) fres];
elseif nfwin == 2
    fwin = [fwin DFLT_FRES];
end

%----- Get time bounds
ntwin = max(size(twin));
if ntwin <= 1
    if ntwin == 1
        tres = twin;
    else
        tres = DFLT_TRES; % Default time resolution
    end
    twin = [freq.stime freq.etime tres];
elseif ntwin == 2
    twin = [twin DLFT_TRES];
end

fw=fwin(2)-fwin(1);
tw=twin(2)-twin(1);

%----- Construct the ploting matrix
spec=zeros(fwin(3),twin(3));
nc=get(freq,'nc');
for i=1:nc
    % Skip this component if it will not show in the window
    % (messy because old MATLABs lacked short-circuit evaluation.)
    if freq.d(i).stime > twin(2), continue, end % starts after window
    if freq.d(i).etime < twin(1), continue, end % ends before window
    if min(freq.d(i).c) > fwin(2), continue, end % min frequency above window
    if max(freq.d(i).c) < fwin(1), continue, end % max frequency below window
    
    npts = max(size(freq.d(i).c));
    treal = linspace(freq.d(i).stime, freq.d(i).etime, npts);
    dt = get(freq,'dt',i);
    sidx=floor((twin(1)-freq.d(i).stime)/dt)+1;
    eidx=ceil((twin(2)-freq.d(i).stime)/dt)+1;
    if sidx < 1; sidx=1; end
    if eidx > npts, eidx=npts; end
    fc = freq.d(i).c;
    ac = amp.d(i).c;
    for x=sidx:eidx
        f = fc(x);
        freqidx = round((fwin(3)-1)*(f-fwin(1))/fw+1);
        if freqidx < 1, continue, end
        if freqidx <= fwin(3)
            % deweight amplitude by number of points per wave to get a
            % spectral *density*. FIXME: This is missing a constant factor.
            a = ac(x) * f;
            tx = round((twin(3)-1)*(treal(x)-twin(1))/tw+1);
            spec(freqidx,tx)=spec(freqidx,tx)+a;
        end
    end
end

%----- Create the output scales
tscale=linspace(twin(1),twin(2),twin(3));
fscale=linspace(fwin(1),fwin(2),fwin(3))';

%----- Plot if no output arguments are requested
if nargout == 0
    img(tscale,fscale,spec);
    clear spec;
end