crtrl

 FUNCTION crtrl() implements the CRTRL algorithm

 Based on the paper of "Nonlinear adaptive prediction of complex-valued signals by complex-valued PRNN", 
 IEEE Transactions on Signal Processing, vol 53, no 5, 2005.


 INPUT:
 input: input signal which should be scaled according to the dynamic range of nonlinearity 
 target: desired signal 
 beta: value of slope of nonlinearity
 mode: choice of nonlinearity
 N: filter length
 p: length of input vector
 mode: choice of nonlinearity

 OUTPUT:
 y: output signal
 e: output error signal
 W: weight vector of adaptive filter 

 Complex Valued Nonlinear Adaptive Filtering toolbox for MATLAB
 Supplementary to the book:
 
 "Complex Valued Nonlinear Adaptive Filters: Noncircularity, Widely Linear and Neural Models"
 by Danilo P. Mandic and Vanessa Su Lee Goh
 
 (c) Copyright Danilo P. Mandic 2009
 http://www.commsp.ee.ic.ac.uk/~mandic
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    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 2 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 can obtain a copy of the GNU General Public License from
    http://www.gnu.org/copyleft/gpl.html or by writing to
    Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ...........................................

0001 % FUNCTION crtrl() implements the CRTRL algorithm
0002 %
0003 % Based on the paper of "Nonlinear adaptive prediction of complex-valued signals by complex-valued PRNN",
0004 % IEEE Transactions on Signal Processing, vol 53, no 5, 2005.
0005 %
0006 %
0007 % INPUT:
0008 % input: input signal which should be scaled according to the dynamic range of nonlinearity
0009 % target: desired signal
0010 % beta: value of slope of nonlinearity
0011 % mode: choice of nonlinearity
0012 % N: filter length
0013 % p: length of input vector
0014 % mode: choice of nonlinearity
0015 %
0016 % OUTPUT:
0017 % y: output signal
0018 % e: output error signal
0019 % W: weight vector of adaptive filter
0020 %
0021 % Complex Valued Nonlinear Adaptive Filtering toolbox for MATLAB
0022 % Supplementary to the book:
0023 %
0024 % "Complex Valued Nonlinear Adaptive Filters: Noncircularity, Widely Linear and Neural Models"
0025 % by Danilo P. Mandic and Vanessa Su Lee Goh
0026 %
0027 % (c) Copyright Danilo P. Mandic 2009
0028 % http://www.commsp.ee.ic.ac.uk/~mandic
0029 %
0030 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0031 %    This program is free software; you can redistribute it and/or modify
0032 %    it under the terms of the GNU General Public License as published by
0033 %    the Free Software Foundation; either version 2 of the License, or
0034 %    (at your option) any later version.
0035 %
0036 %    This program is distributed in the hope that it will be useful,
0037 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
0038 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0039 %    GNU General Public License for more details.
0040 %
0041 %    You can obtain a copy of the GNU General Public License from
0042 %    http://www.gnu.org/copyleft/gpl.html or by writing to
0043 %    Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0044 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0045 % ...........................................
0046 function [y, e, W, weightRecord] = crtrl(input, target, beta, mode, p, N, monte)
0047 
0048 pp = length(input);            % Number of input signal
0049 
0050 ita=0.1;
0051 
0052 s = zeros(p,1) + i*zeros(p,1);                         % The input vector
0053 DW = zeros(p+N+1, N) + i*zeros(p+N+1, N);                    % Change applied to the weight matrix for neurons 1&2
0054 
0055 
0056 
0057 % Weights vectors
0058 W = 0.1*randn(p+N+1,N) + i*0.1*randn(p+N+1,N);                     % The weight matrix
0059 PI_prev = zeros(N*N, p+N+1) + i*zeros(N*N, p+N+1);  PI_curr = zeros(N*N, p+N+1)+ i*zeros(N*N, p+N+1);  % Jacobian matrix
0060 
0061 
0062 weightRecord = zeros(p+N+1,pp) + i*zeros(p+N+1,pp);
0063 P_record = zeros(p+N+1,p+N+1,pp) + i*zeros(p+N+1,p+N+1,pp);
0064 
0065 E = zeros(pp,1) + i*zeros(pp,1);
0066 MSE = zeros(pp,1) + i*zeros(pp,1);
0067 % Output vectors from the neurons
0068 Y_curr = zeros(pp, N) + i*zeros(pp, N);
0069 
0070 U = zeros(p+N+1, 1) + i*zeros(p+N+1, 1);              % INput to the network(external,feedback and bias)
0071 V = zeros(1,N) + i*zeros(1,N);                        % weighted input vector to a neuron
0072 Errors = 0.001*(randn(pp, 1) + i*randn(pp, 1));
0073 Phi = zeros(1, N) + i*zeros(1,N);
0074 
0075 for Monte=1:monte
0076 for k=2:pp
0077     s = [input(k-1);s(1:p-1)];
0078     U(:) = [s;1;Y_curr(k-1,:).'];      % The input vector to the network
0079     V  = U.'*W;
0080         
0081     switch (mode)
0082         case 1
0083             %Y_curr(k,:)     = logsig(beta*V);
0084              Y_curr(k,:)     = 1./(1+exp(-beta*V));
0085         case 2
0086             %Y_curr(k,:)     = tansig( beta * V);
0087              Y_curr(k,:)     = tanh( beta * V);
0088     end
0089 
0090     Errors(k)     = target(k)-Y_curr(k, 1) ;   %taking the delay of observations to minus the y output
0091     
0092     switch(mode)
0093         case 1
0094         %Phi = dlogsig(beta*V,Y_curr(k,:));
0095         Phi = beta * (1./(1+exp(-beta*V))).*( 1 - 1./(1+exp(-beta*V)) );
0096         case 2
0097         %Phi = dtansig(beta * V,Y_curr(k,:));
0098         Phi = 1 - (tanh(beta * V)).^2;
0099     end
0100 
0101     for j=1:N      %Put N PI matrix into one Big PI, that's N x [(1+p+N) x N] size
0102         for n=1:N
0103             for l=1:p+1+N
0104                 temp = 0;
0105                 for m=1:N
0106                     temp = temp + conj(W(1+p+m, j))*PI_prev((m-1)*N+n, l);
0107                 end
0108                 wtemp(n,l) = temp;
0109                 if n==j
0110                     wtemp(n, l) = wtemp(n, l) + conj(U(l));
0111                 end
0112                 PI_curr((j-1)*N+n, l) = conj(Phi(j)) * wtemp(n, l);
0113             end
0114         end
0115     end
0116     PI_prev = PI_curr;
0117     DW = ita*Errors(k)*PI_curr(1:N,:).';
0118     W=W+DW;
0119     weightRecord(:, k) = W(:,1);
0120 end
0121 
0122 E = E +(1/2)*abs(Errors).^2;
0123 MSE = MSE+10*log10((1/2)*abs(Errors).^2);
0124 var_error(Monte) = var(Errors(1:end));
0125 var_signal(Monte) = var(input(1:end));
0126 Rp(Monte) = 10*log10(var_signal(Monte)/var_error(Monte));
0127 
0128 end
0129 %keyboard;
0130 y = Y_curr(:,1);
0131 e = Errors;
0132 
0133   error = norm(y(1:end).' - target(1:end)); % check on this
0134 
0135   fprintf('EKF error    = %6.2f',error)
0136   fprintf('\n')
0137   
0138   fprintf('Rp     = %6.2f',mean(Rp))
0139   fprintf('\n')
0140   
0141   % PLOT RESULTS:
0142   % ============
0143                                             
0144   figure(1)
0145 
0146   plot(1:length(y),abs(y),'r',1:length(input),abs(input),'b');
0147   ylabel('Prediction','fontsize',15)
0148   %axis([0 pp -1.5 1.5]);
0149   
0150   figure(3)
0151   plot(E)
0152   %save Emodel2_rtrl E
0153   
0154   figure(4)
0155   plot(MSE)
0156