www.gusucode.com > signal 工具箱matlab源码程序 > signal/fwht.m

    function y = fwht(x,N,ordering)
%FWHT Fast Discrete Walsh-Hadamard Transform
%   Y = FWHT(X) is the discrete Walsh-Hadamard transform of vector X. The 
%   transform coefficients are stored in Y. If X is a matrix, the function 
%   operates on each column.  
%
%   Y = FWHT(X,N) is the N-point discrete Walsh-Hadamard transform of the 
%   vector X where N must be a power of two. X is padded with zeros if it 
%   has less than N points and truncated if it has more. The default value
%   of N is equal to the length of the vector X if it is a power of two or
%   the next power of two greater than the signal vector length. The
%   function errors out if N is not equal to a power of two.
%
%   Y = FWHT(X,N,ORDERING) or FWT(X,[],ORDERING) specifies the order of the
%   Walsh-Hadamard transform coefficients. ORDERING can be 'sequency', 
%   'hadamard' or 'dyadic'. Default ORDERING type is 'sequency'.
%
%   
%   EXAMPLES:
%
%   % Example 1: Walsh-Hadamard transform of a signal made up of Walsh 
%                % functions
%                w1 = [1 1 1 1 1 1 1 1];
%                w2 = [1 1 1 1 -1 -1 -1 -1];
%                w3 = [1 1 -1 -1 -1 -1 1 1];
%                w4 = [1 1 -1 -1 1 1 -1 -1];
%                x = w1 + w2 + w3 + w4; % signal formed by adding Walsh 
%                                       % functions  
%                y = fwht(x); % first four values of y should be non-zero 
%                             % equal to one
%
%   % Example 2: Walsh-Hadamard transform - 'hadamard' function and ordering
%                w = hadamard(8); % Walsh functions in Hadamard order
%                x = w(:,1) + w(:,2) + w(:,3) + w(:,4);
%                y = fwht(x,[],'hadamard'); % first four values equal to one
%
%   For more information see the <a href="matlab:web(fullfile(docroot,'signal/examples/discrete-walsh-hadamard-transform.html'))">Discrete Walsh-Hadamard Transform Example</a> 
%   or enter "doc fwht" at the MATLAB command line.
%   See also IFWHT, FFT, IFFT, DCT, IDCT, DWT, IDWT.

%   Copyright 2008-2013 The MathWorks, Inc.

% error out if number of input arguments is not between 1 and 3
narginchk(1,3)
if isempty(x)
    y = [];
    return
end
% check optional inputs' specifications and/or make default assignments
if nargin < 3
    ordering = 'sequency';  % default ordering is sequency
    if nargin < 2
        N = defaulttransformlength(size(x));
    end
end

% Cast to enforce precision rules. 
N = signal.internal.sigcasttofloat(N,'double','fwht','N','allownumeric');
signal.internal.sigcheckfloattype(x,'','fwht','X');

if isempty(N)
  N = defaulttransformlength(size(x));
end
validatetransformlength(N);
validateordering(ordering);

% do pre-processing on input signal if necessary
[x,tFlag] = preprocessing(x,N,ordering);

% calculate first stage coefficients and store in x
for i = 1:2:N-1
    x(i,:) = x(i,:) + x(i+1,:);
    x(i+1,:) = x(i,:) - 2 * x(i+1,:);
end

L = 1;
% same data type as x to enforce precision rules
y = zeros(size(x),class(x)); %#ok<ZEROLIKE> 
for nStage = 2:log2(N) % log2(N) = number of stages in the flow diagram
    % calculate coefficients for the ith stage specified by nStage
    M = 2^L;
    J = 0; K = 1;
    while (K < N)
        for j = J+1:2:J+M-1
            if strcmpi(ordering,'sequency') == 1
                y(K,:) = x(j,:) + x(j+M,:);
                y(K+1,:) = x(j,:) - x(j+M,:);
                y(K+2,:) = x(j+1,:) - x(j+1+M,:);
                y(K+3,:) = x(j+1,:) + x(j+1+M,:);
            else
                y(K,:) = x(j,:) + x(j+M,:);
                y(K+1,:) = x(j,:) - x(j+M,:);
                y(K+2,:) = x(j+1,:) + x(j+1+M,:);
                y(K+3,:) = x(j+1,:) -  x(j+1+M,:);
            end
            K = K + 4;
        end
        J = J + 2*M;
    end
    % store coefficients in x at the end of each stage
    x = y;
    L = L + 1;
end
% perform scaling of coefficients
y = x ./ N;
if tFlag
    y = transpose(y);
end

%--------------------------------------------------------------------------

function validatetransformlength(N)
% check if transform length is specified correctly - positive scalar and
% power of two
% sigDim = [rows cols] of input signal

if ~isempty(find(char(num2str(log2(N)))=='.', 1)) || (N < 1)
    error(message('signal:fwht:InvalidTransformLength'));
end
 
%--------------------------------------------------------------------------

function N = defaulttransformlength(sigDim) 
% this function assigns default transform length as signal length if it is
% a power of two or the next power of two greater than signal length
% sigDim = [rows cols] of input signal

if min(sigDim) == 1
    N = max(sigDim);
else
    N = sigDim(1);
end
if ~isempty(find(char(num2str(log2(N)))=='.', 1))
    p = ceil(log2(N));
    N = 2^p;
end

%--------------------------------------------------------------------------

function validateordering(ordering) 
% this functions checks if the specified ordering is a string input and one
% of sequency, hadamard or dyadic'
if ~ischar(ordering) || ~any(strcmpi(ordering,{'sequency','hadamard','dyadic'}))
    error(message('signal:fwht:InvalidOrderingType', 'ORDERING', 'sequency', 'hadamard', 'dyadic'));
end
    
%--------------------------------------------------------------------------

function [x,tFlag] = preprocessing(x,N,ordering)
% this function performs zero-padding, truncation or input bit-reversal if
% necessary. NROWS amd MCOLS specify the output orientation which is kept
% same as that of input.

tFlag = false;
if size(x,1) == 1
    x = x(:); % column vectorizing input sequence
    tFlag = true;
end
n = size(x,1);    
if n < N
    % zero-pad
    x(n+1:N,:) = 0;
elseif n > N
    % truncate
    x = x(1:N,:);
end
% Re-arrange input in bit-reversed order if ordering is hadamard
if strcmpi(ordering,'hadamard') == 1
    x = bitrevorder(x);
end

%--------------------------------------------------------------------------