function A=hmapdn(f)
%
% Copyright(c) D.E.Newland 1992
% All right reserved
%
N=length(f);
n=round(log(N)/log(2));
a=hwtdn(f);
%
for j=1:2^(n-2)%為什么進行諧波變換后數據量減少為1/4
    A(1,j)=(abs(a(1))^2);%log
    A(n+1,j)=(abs(a(2^(n-1)+1))^2);%log
end
for j=2:n
    for k=1:2^(j-2)
        for m=1:2^(n-j)
            A(j,(k-1)*2^(n-j)+m)=(abs(a(2^(j-2)+k))^2+abs(a(N+2-2^(j-2)-k))^2);%log
        end
    end
end