reassign.m

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

function sr=reassign(s,itime,ifreq,a,p5);
%REASSIGN  Reassign time-frequency distribution.
%   Usage:  sr = reassign(s,itime,ifreq,a);
%           sr = reassign(s,itime,ifreq,a,'aa');
%
%   REASSIGN(s,itime,ifreq,a) will reassign the values of the positive
%   time-frequency distribution s using the instantaneous time and frequency
%   itime and ifreq. The lattice is determined by the time shift a and
%   the number of channels deduced from the size of s.
%
%   itime and ifreq can be obtained by the routines INSTTF_DGT,
%   INSTTF_ABS and INSTTF_PHASE.
%
%   REASSIGN(s,itime,ifreq,'aa') will do the same using an anti-aliasing
%   algorithm to do the reassignment.
%
%   The standard way of calling this routine to generate a reassigned
%   spectrogram from a signal f is:
% 
%      [itime, ifreq, c] = insttf_dgt(f,'gauss',a,M);
%      sr = reassign(abs(c).^2,itime,ifreq,a);
%   
%   SEE ALSO: RESGRAM, INSTFTF_DGT, INSTTF_ABS, INSTTF_PHASE
% 
%   REFERENCES:
%     F. Auger and P. Flandrin. Improving the readability of time-frequency
%     and time-scale representations by the reassignment method.
%     IEEETrans.SignalProcess., 43(5):1068-1089, 1995.

% 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. Soendergaard, 2008.
  
error(nargchk(4,5,nargin));

[M,N,W]=size(s);
L=N*a;
b=L/M;

freqpos=fftindex(M);  
for w=1:W
  ifreq(:,:,w)=ifreq(:,:,w)/b+repmat(freqpos,1,N);
end;

timepos=fftindex(N);
for w=1:W
  itime(:,:,w)=itime(:,:,w)/a+repmat(timepos.',M,1);
end;

ifreq=round(ifreq);
itime=round(itime);

ifreq=mod(ifreq,M);
itime=mod(itime,N);  
  
%t = repmat((0:L-1),M,1);
%O = repmat((0:L-1).',1,N);

sr=zeros(M,N,W);

if nargin==5
  switch(lower(p5))
   case {'aa'}
    
    % --- use antialiasing ----
    
    fk=mod(floor(ifreq),M)+1;
    ck=mod(ceil(ifreq),M)+1;
    fn=mod(floor(itime),N)+1;
    cn=mod(ceil(itime),N)+1;
    
    alpha = itime-floor(itime);
    beta  = ifreq-floor(ifreq);
    m1 =(1-alpha).*(1-beta).*s;
    m2 =(1-alpha).*beta.*s;
    m3 =alpha.*(1-beta).*s;
    m4 =alpha.*beta.*s;
    for ii=1:M
      for jj=1:N
        sr(fk(ii,jj),fn(ii,jj))=sr(fk(ii,jj),fn(ii,jj))+m1(ii,jj);
        sr(ck(ii,jj),fn(ii,jj))=sr(ck(ii,jj),fn(ii,jj))+m2(ii,jj);
        sr(fk(ii,jj),cn(ii,jj))=sr(fk(ii,jj),cn(ii,jj))+m3(ii,jj);
        sr(ck(ii,jj),cn(ii,jj))=sr(ck(ii,jj),cn(ii,jj))+m4(ii,jj);
        
      end;
    end;
  end;
else
  
  % --- no antialiasing, just round ---
  itime=itime+1;
  ifreq=ifreq+1;
  
  for ii=1:M
    for jj=1:N      
      sr(ifreq(ii,jj),itime(ii,jj)) = sr(ifreq(ii,jj),itime(ii,jj))+s(ii,jj);
    end;
  end;  
  
end;