【Hammerstein模型的级联】快速估计构成一连串哈默斯坦模型的结构元素研究（Matlab代码实现）

function hhat = Hammerstein_ID(input_sig,output,duration,f1,f2,fs,N,opt_meth,opt_filt)
%---------------------------------------------------------
%
% hhat = Hammerstein_ID(input_sig,output,f1,f2,fs,N,opt_meth,opt_filt)
%
% Estimates the Kernels "h" of the cascade of Hammerstein model of order N fed with
% the input signal "input" and where the corresponding output signal "output"
% has been recorded. "input" has to be an exponential sine sweep going from
% f1 to f2.
%
% Input parameters:
%   input_sig  : input exponential sine sweep
%   output : output of the system fed with the input signal
%   f1 : starting frequency of the sweep
%   f2 : end frequency of the sweep
%   fs : sampling frequency
%   N  : Order of the model to be identified
%   opt_meth : Method to use for the estimation (string expected)
%       - 'Reb': Method proposed by R閎illat et al. in [1]
%       - 'Nov': Method proposed by Novak et al in [2]
%   opt_filt : Specifies the method to use to compute the inverse filter 
%       (string expected). By default 'TFB_linear' is chosen.
%       - 'TFB_square': FTT based filter with a square window and
%         regularization (see [1])
%       - 'TFB_linear': FTT based filter with a square window with continuous
%         linear borders and regularization (see [1])
%       - 'TFB_gevrey': FTT based filter with a square window with infinitely
%         continuous gevrey borders and regularization (see [1])
%       - 'Nov' : Filter based on the analytical formulation using aymptotic 
%         signals (see [2]).
%
% Output:
%   h : 2D matrix containing the pseudo impulse responses (temporal domain) 
%   of the estimated kernels.
display('--> Hammerstein Identification in progress ...')
% Check arguments
if nargin<6
    display('    => ERROR : Incorrect number of arguments')
    return
elseif nargin<7
    display('    => No method option and filtering option specified. ')
    display('    => Method option = ''Reb'' chosen by default.')
    display('    => Filtering option = ''TFB_linear'' chosen by default.')
    opt_meth = 'Reb' ;
    opt_filt = 'TFB_linear';
elseif nargin<8
    display(['    => Method ' opt_meth ' chosen'])
    display('    => No filtering option specified. ')
    if strcmp(opt_meth,'Reb')
        opt_filt = 'TFB_linear';
        display('    => Filtering option = ''TFB_linear'' chosen by default.')
    elseif strcmp(opt_meth,'Nov')
        opt_filt = 'Nov';
        display('    => Filtering option = ''Nov'' chosen by default.')
    else 
        display('    => ERROR : Unknown method option')
        display('    => Select ''Reb'' or ''Nov''')
        return
    end
else
    if ( strcmp(opt_meth,'Nov') || strcmp(opt_meth,'Reb'))
        display(['    => Method ' opt_meth ' chosen'])
    else
        display('    => ERROR : Unknown method option')
        display('    => Select ''Reb'' or ''Nov''')
        return
    end
    if ( strcmp(opt_filt,'TFB_square') || strcmp(opt_filt,'TFB_linear') || strcmp(opt_filt,'TFB_gevrey') || strcmp(opt_filt,'Nov'))
        display(['    => Filtering ' opt_filt ' chosen'])
    else
        display('    => ERROR : Unknown filtering option')
        display('    => Select ''TFB_square'', ''TFB_linear'', ''TFB_gevrey'' or ''Nov''')
        return
    end
end
% Equivalent pulsations
w1 = f1/fs*2*pi;
w2 = f2/fs*2*pi;
% Convolution of the response with the inverse of the sweep
if strcmp(opt_meth,'Reb')
    inverse_input_sig = compute_inverse_filter(input_sig,f1,f2,fs,opt_filt) ;
    gToCut = convq(output,inverse_input_sig);
elseif strcmp(opt_meth,'Nov')
    % Nonlinear convolution in the spectral domain
    gToCut = nonlinear_convolution(output,duration,f1,f2,fs);        
    gToCut = [gToCut; gToCut];
end
% Computation of the delay of the pseudo RI
if strcmp(opt_meth,'Reb')
    T = length(input_sig); % Actual length of the sweep (in samples)
    deltaT = T*log(1:N)/log(w2/w1);
elseif strcmp(opt_meth,'Nov')
    T = length(output); % Actual length of the output (in samples)
    L = 1/f1*round( (duration*f1)/(log(f2/f1)) ); 
    deltaT = L*log(1:N)*fs;
end