spreadop.m

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

function h=spreadop(f,coef)
%SPREADOP  Spreading operator
%   Usage: h=spreadop(f,c);
%
%   SPREADOP(f,c) applies the spreading operator with symbol c to the
%   input f. c must be square.
%
%   SPREAOP(f,c) computes the following for c of size LxL:
%
%             L-1 L-1 
%    h(l+1) = sum sum c(m+1,n+1)*exp(-2*pi*i*l*m/L)*f(l-n+1)
%             n=0 m=0
%
%   where l=0,...,L-1 and l-n is computed modulo L.
%
%   The combined symbol of two spreading operators can be found by
%   using TCONV. Consider two symbols c1 and c2 and define f1 and f2 by:
% 
%     h  = tconv(c1,c2)
%     f1 = spreadop(spreadop(f,c2),c1);
%     f2 = spreadop(f,h);
% 
%   then f1 and f2 are equal.
%
%   SEE ALSO:  TCONV, SPREADFUN, SPREADINV, SPREADADJ

% 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/>.

error(nargchk(2,2,nargin));

if ndims(coef)>2 || size(coef,1)~=size(coef,2)
    error('Input symbol coef must be a square matrix.');
end;

L=size(coef,1);

% Change f to correct shape.
[f,Ls,W,wasrow,remembershape]=comp_sigreshape_pre(f,'DGT',0);

f=postpad(f,L);

h=zeros(L,W);

if issparse(coef) && nnz(coef)<L

    % The symbol is so sparse that the straighforward definition is
    % the fastest way to apply it.

    [mr,nr,cv]=find(coef);
    h=zeros(L,W);

    % We need mr and nr to be zero-indexed
    mr=mr-1;
    nr=nr-1;

    % This is the basic idea of the routine below
    %for ii=1:length(mr)
    %  for l=0:L-1
    %	h(l+1,:)=h(l+1,:)+cv(ii)*exp(-2*pi*i*l*mr(ii)/L)*f(mod(l-nr(ii),L)+1,:);
    %  end;
    %end;

    l=(-2*pi*i*(0:L-1)/L).';
    for ii=1:length(mr)
        bigmod=repmat(exp(l*mr(ii)),1,W);
        h=h+cv(ii)*(bigmod.*circshift(f,nr(ii)));
    end;

else

    if 0

        % This version explicitly constructs the matrix representation T
        % and then applies this matrix as the final step.
        coef=fft(coef);

        T=zeros(L);

        for ii=0:L-1
            for jj=0:L-1
                T(ii+1,jj+1)=coef(ii+1,mod(ii-jj,L)+1);
            end;
        end;

        h=T*f;

    end;

    if 1

        % This version only touches coef one column at a time, and it suited
        % if coef is sparse.

        for n=0:L-1
            cf=fft(coef(:,n+1));
            for l=0:L-1
                for w=1:W
                    h(l+1,w)=h(l+1,w)+cf(l+1)*f(mod(l-n,L)+1,w);
                end;
            end;
        end;


    end;

end;