insttf_phase.m

Matlab时频分析工具箱,希望能对大家有所帮助啊

function [itime,ifreq]=insttf_phase(cphase,a)
%INSTFREQ_PHASE  Instantaneous frequency computed from phase
%           [itime, ifreq] = instfreq(cphase,a);
%
%   [itime,ifreq]=INSTTF_PHASE(cphase,a) computes the instantanous time and
%   frequency from the phase cphase of a DGT of the signal. The original
%   DGT from which the phase is obtained must have been computed using a
%   time-shift of a.
%
%   See the help text on INSTTF_DGT for a description of the output itime
%   and ifreq.
%
%   If the original signal is known, the instantaneous time and frequency
%   can be more accuratly computed by INSTTF_DGT
% 
%   SEE ALSO: INSTTF_DGT, INSTTF_ABS, REASSIGN
% 
%   REFERENCES:
%     J. Flanagan, D. Meinhart, R. Golden, and M. Sondhi. Phase Vocoder. The
%     Journal of the Acoustical Society of America, 38:939, 1965.
%     

% 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/>.
  
if ~isreal(cphase)
  error(['Input phase must be real valued. Use the "angle" function to ' ...
         'compute the argument of complex numbers.']);
end;
  
% --- linear method ---
[M,N,W]=size(cphase);
L=N*a;
b=L/M;

if 0

  % This is the classic phase vocoder algorithm by Flanagan.
  
  ifreq = cphase-circshift(cphase,[0,-1]);
  ifreq = ifreq- 2*pi*round(ifreq/(2*pi));
  ifreq = -ifreq/(2*pi)*L;
    
  % Phase-lock the angles.
  TimeInd = (0:(N-1))*a;
  FreqInd = (0:(M-1))/M;
  
  phaselock = FreqInd'*TimeInd;
  cphase = cphase+2*pi.*phaselock;
  
  itime = cphase-circshift(cphase,[1,0]);
  itime = itime- 2*pi*round(itime/(2*pi));
  itime = -itime/(2*pi)*L;
  
end;


if 1
  % This is the classic phase vocoder algorithm by Flanagan modified to
  % yield a second order centered difference approximation.
  
  % Forward approximation
  ifreq_1 = cphase-circshift(cphase,[0,-1]);
  ifreq_1 = ifreq_1 - 2*pi*round(ifreq_1/(2*pi));
  % Backward approximation
  ifreq_2 = circshift(cphase,[0,1])-cphase;
  ifreq_2 = ifreq_2 - 2*pi*round(ifreq_2/(2*pi));
  % Average
  ifreq = (ifreq_1+ifreq_2)/2;
  
  ifreq = -ifreq/(2*pi*a)*L;
  
  % Phase-lock the angles.
  TimeInd = (0:(N-1))*a;
  FreqInd = (0:(M-1))/M;
  
  phaselock = FreqInd'*TimeInd;
  cphase = cphase+2*pi.*phaselock;
  
  % Forward approximation
  itime_1 = cphase-circshift(cphase,[-1,0]);
  itime_1 = itime_1 - 2*pi*round(itime_1/(2*pi));
  % Backward approximation
  itime_2 = circshift(cphase,[1,0])-cphase;
  itime_2 = itime_2 - 2*pi*round(itime_2/(2*pi));
  % Average
  itime = (itime_1+itime_2)/2;
  
  itime = itime/(2*pi*b)*L;
  
end;