function varargout=instfreqplot(f,varargin)
%INSTFREQPLOT Plot of instantaneous frequency
% Usage: instfreqplot(f,op1,op2, ... );
% instfreqplot(f,fs,op1,op2, ... );
%
% INSTFREQPLOT(f) plots the instantaneous frequency of f using a DGT.
%
% INSTFREQPLOT(f,fs) does the same for a signal with sampling rate
% fs Hz.
%
% The instantaneous frequency contains extreme spikes in regions
% where the spectrogram is close to zero. These points are usually
% uninteresting and destroy the visibility of the plot. Use the 'thr'
% or 'clim' or 'climsym' options (see below) to remove these points.
%
% An example:
%
% instfreqplot(greasy,16000,'thr',.03,'climsym',100);
%
% will produce a nice instantaneous frequency plot of the GREASY
% signal.
%
% INSTFREQPLOT accepts the following optional arguments:
%
% 'tfr',v Set the ratio of frequency resolution to time resolution.
% A value v=1 is the default. Setting v>1 will give better
% frequency resolution at the expense of a worse time
% resolution. A value of 0<v<1 will do the opposite.
%
% 'wlen',s Window length. Specifies the length of the window
% measured in samples. See help of PGAUSS on the exact
% details of the window length.
%
% 'thr',r Keep the coefficients with a magnitude larger than r*
% times the largest magnitude. Set the instantaneous
% frequency of the rest of the coefficients to zero
%
% 'nf' Display negative frequencies, with the zero-frequency
% centered in the middle. For real signals, this will just
% mirror the upper half plane. This is standard for complex
% signals.
%
% 'tc' Time centering. Move the beginning of the signal to the
% middle of the plot. This is useful for visualizing the
% window functions of the toolbox.
%
% 'image' Use imagesc to display the spectrogram. This is the
% default.
%
% 'dgt' Use the 'dgt' method to compute the instantaneous
% frequency. This is the default.
%
% 'clim',clim Use a colormap ranging from clim(1) to clim(2). These
% values are passed to imagesc. See the help on imagesc.
%
% 'climsym',cval Use a colormap ranging from -cval to cval
%
% 'fmax',y Display y as the highest frequency.
%
% See also: sgram, resgram, phaseplot
%
% Url: http://ltfat.github.io/doc/gabor/instfreqplot.html
% Copyright (C) 2005-2023 Peter L. Soendergaard <peter@sonderport.dk> and others.
% This file is part of LTFAT version 2.6.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% AUTHOR: Peter L. Søndergaard
% REFERENCE: NA
% TESTING: NA
if nargin<1
error('Too few input arguments.');
end;
if sum(size(f)>1)>1
error('Input must be a vector.');
end;
% Define initial value for flags and key/value pairs.
definput.flags.wlen={'nowlen','wlen'};
definput.flags.tc={'notc','tc'};
definput.flags.method={'dgt','phase','abs'};
definput.flags.clim={'noclim','clim','climsym'};
if isreal(f)
definput.flags.posfreq={'posfreq','nf'};
else
definput.flags.posfreq={'nf','posfreq'};
end;
definput.keyvals.tfr=1;
definput.keyvals.wlen=0;
definput.keyvals.thr=[];
definput.keyvals.clim=[0,1];
definput.keyvals.climsym=1;
definput.keyvals.fmax=[];
definput.keyvals.xres=800;
definput.keyvals.yres=600;
definput.keyvals.fs=[];
[flags,kv,fs]=ltfatarghelper({'fs'},definput,varargin);
% Override the setting from tfplot, because INSTFREQPLOT does not support the
% dB-settings
flags.do_lin=1;
flags.do_db=0;
flags.do_dbsq=0;
dofs=~isempty(fs);
% Downsample
if ~isempty(kv.fmax)
if ~isempty(kv.fs)
resamp=kv.fmax*2/fs;
else
resamp=kv.fmax*2/length(f);
end;
f=fftresample(f,round(length(f)*resamp));
kv.fs=2*kv.fmax;
end;
Ls=length(f);
if flags.do_posfreq
kv.yres=2*kv.yres;
end;
try
[a,M,L,N,Ndisp]=gabimagepars(Ls,kv.xres,kv.yres);
catch
err=lasterror;
if strcmp(err.identifier,'LTFAT:noframe')
error(sprintf(['The signal is too long. INSTFREQPLOT cannot visualize all the details.\n' ...
'Try a shorter signal or increase the image resolution by calling:\n\n' ...
' sgram(...,''xres'',xres,''yres'',yres);\n\n' ...
'for larger values of xres and yres.\n'...
'The current values are:\n xres=%i\n yres=%i'],kv.xres,kv.yres));
else
rethrow(err);
end;
end;
b=L/N;
% Set an explicit window length, if this was specified.
if flags.do_wlen
kv.tfr=kv.wlen^2/L;
end;
g={'gauss',kv.tfr};
if flags.do_dgt
[coef,fgrad,dgtcoef]=gabphasegrad('dgt',f,g,a,M);
end;
if flags.do_phase
dgtcoef=dgt(f,g,a,M);
[coef,fgrad]=gabphasegrad('phase',angle(dgtcoef),a);
end;
if flags.do_abs
dgtcoef=dgt(f,g,a,M);
[coef,fgrad]=gabphasegrad('abs',abs(dgtcoef),g,a);
end;
if ~isempty(kv.thr)
% keep only the largest coefficients.
maxc=max(abs(dgtcoef(:)));
mask=abs(dgtcoef)<maxc*kv.thr;
coef(mask)=0;
end
% Cut away zero-extension.
coef=coef(:,1:Ndisp);
% Convert climsym into a normal clim
if flags.do_climsym
flags.do_clim=1;
kv.clim=[-kv.climsym,kv.climsym];
end;
if flags.do_posfreq
coef=coef(1:floor(M/2)+1,:);
plotdgtreal(coef,a,M,'argimport',flags,kv);
else
plotdgt(coef,a,'argimport',flags,kv);
end;
if nargout>0
varargout={coef};
end;