FUNCTION split_complex_rtrl() implements the Split-Complex RTRL algorithm
Based on the paper of SuLee, "Nonlinear adaptive prediction of complex-valued signals by complex-valued PRNN", IEEE Transactions on Signal Processing, vol 53, no 5, 2005.
INPUT:
input: Signal which should be scaled according to the dynamic range of nonlinearity
OUTPUT:
Elog: Error in dB
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 split_complex_rtrl() implements the Split-Complex RTRL algorithm 0002 % 0003 % Based on the paper of SuLee, "Nonlinear adaptive prediction of complex-valued signals by complex-valued PRNN", IEEE Transactions on Signal Processing, vol 53, no 5, 2005. 0004 % INPUT: 0005 % input: Signal which should be scaled according to the dynamic range of nonlinearity 0006 % 0007 % OUTPUT: 0008 % Elog: Error in dB 0009 % 0010 % 0011 % Complex Valued Nonlinear Adaptive Filtering toolbox for MATLAB 0012 % Supplementary to the book: 0013 % 0014 % "Complex Valued Nonlinear Adaptive Filters: Noncircularity, Widely Linear and Neural Models" 0015 % by Danilo P. Mandic and Vanessa Su Lee Goh 0016 % 0017 % (c) Copyright Danilo P. Mandic 2009 0018 % http://www.commsp.ee.ic.ac.uk/~mandic 0019 % 0020 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0021 % This program is free software; you can redistribute it and/or modify 0022 % it under the terms of the GNU General Public License as published by 0023 % the Free Software Foundation; either version 2 of the License, or 0024 % (at your option) any later version. 0025 % 0026 % This program is distributed in the hope that it will be useful, 0027 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0028 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0029 % GNU General Public License for more details. 0030 % 0031 % You can obtain a copy of the GNU General Public License from 0032 % http://www.gnu.org/copyleft/gpl.html or by writing to 0033 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0034 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0035 % ........................................... 0036 function Elog =split_complex_rtrl(input) 0037 0038 %Defining parameters Danilo's book 0039 p=4; %input taps 0040 N=2;%number of neurons 0041 alpha=0.1;%learning rate 0042 len=5000; %length of input sample 0043 0044 %Initialization of weight parameter 0045 w1=zeros(p+N+1,1)+i*zeros(p+N+1,1); % weights of neuron1 +1 for bias yaa... 0046 w2=zeros(p+N+1,1)+i*zeros(p+N+1,1); % weights of neuron2 0047 W=zeros(p+N+1,N)+i*zeros(p+N+1,N); %total weights matrix 0048 0049 %Initialization of weight change 0050 dW1=zeros(1,p+N+1)+i*zeros(1,p+N+1); %weight change +1 for bias also 0051 dW2=zeros(1,p+N+1)+i*zeros(1,p+N+1); %weight change 0052 DW=zeros(N,p+N+1)+i*zeros(N,p+N+1); 0053 0054 %PI for split case 0055 PII_old=zeros(p+N+1,N,N);%previous sample imaginary 0056 PII_new=zeros(p+N+1,N,N);%current sample imaginary 0057 PRR_old=zeros(p+N+1,N,N);%previous sample real 0058 PRR_new=zeros(p+N+1,N,N);%current sample real 0059 PRI_old=zeros(p+N+1,N,N);%previous sample imaginary 0060 PRI_new=zeros(p+N+1,N,N);%current sample imaginary 0061 PIR_old=zeros(p+N+1,N,N);%previous sample real 0062 PIR_new=zeros(p+N+1,N,N);%current sample real 0063 0064 %Output Matrix 0065 Y_old1=zeros(len+1,1)+i*zeros(len+1,1); % Previous output matrix 1 0066 Y_out1=zeros(len,1)+i*zeros(len,1); % output matrix 1 0067 Y_old2=zeros(len+1,1)+i*zeros(len+1,1); % Previous output matrix 2 0068 Y_out2=zeros(len,1)+i*zeros(len,1); % output matrix 2 0069 0070 %Preparing the basic format fo the input signal 0071 x=zeros(p,1)+i*zeros(p,1);%input signal to the tap 0072 E_dB=zeros(1,len); 0073 Elog=zeros(1,len); 0074 0075 for monte=1:100 0076 monte; 0077 0078 %weight initialization 0079 W_init1real=0.01*rand(p+N+1,1); %initialise the weight1 randomly 0080 W_init1imag=0.01*rand(p+N+1,1); 0081 w1=W_init1real+i*W_init1imag; 0082 W_init2real=0.01*rand(p+N+1,1); %initialise the weight2 randomly 0083 W_init2imag=0.01*rand(p+N+1,1); 0084 w2=W_init2real+i*W_init2imag; 0085 0086 %Load Data of Complex Colored Input 0087 d=input(1:len); 0088 xin(1)=0; 0089 xin(2:len)=d(1:len-1); 0090 0091 %Activity of the Neurons 0092 for k=1:len 0093 x=[xin(k);x(1:p-1)]; %input of the system due to input 0094 Uin=[Y_old1(k);Y_old2(k);1;x]; %the main input to the system, 1 represents the bias 0095 0096 %First Neuron Activity 0097 Vout1=w1.'*Uin; 0098 Sr1=real(w1).'*real(Uin)-imag(w1).'*imag(Uin); 0099 Si1=imag(w1).'*real(Uin)+real(w1).'*imag(Uin); 0100 sig_function1R = 1 ./ ( 1 + exp((-(Sr1))));%Output of the neuron 1 real 0101 sig_function1I= 1 ./ ( 1 + exp( (-(Si1))));%Output of neuron 1 imag 0102 sig_function1= sig_function1R + i*sig_function1I; 0103 sig_function_der1R = sig_function1R .* ( 1 - sig_function1R );%the derivative,f' 0104 sig_function_der1I = sig_function1I .* ( 1 - sig_function1I ); 0105 Y_out1(k)=sig_function1;%store in the output matrix of 1st neuron 0106 Y_old1(k+1)=Y_out1(k); 0107 0108 %Output of 1st Neuron 0109 u1(k)=sig_function1R; 0110 v1(k)=sig_function1I; 0111 0112 %Second Neuron Activity 0113 Vout2=w2.'*Uin; 0114 Sr2=real(w2).'*real(Uin)-imag(w2).'*imag(Uin); 0115 Si2=imag(w2).'*real(Uin)+real(w2).'*imag(Uin); 0116 sig_function2R = 1 ./ ( 1 + exp( -(Sr2)));%Output of the neuron 1 real 0117 sig_function2I= 1 ./ ( 1 + exp( -(Si2)));%Output of neuron 1 imag 0118 sig_function2= sig_function2R + i*sig_function2I; 0119 sig_function_der2R = sig_function2R .* ( 1 - sig_function2R );%the derivative,f' 0120 sig_function_der2I = sig_function2I .* ( 1 - sig_function2I ); 0121 Y_out2(k)=sig_function2;%store in the output matrix of 1st neuron 0122 Y_old2(k+1)=Y_out2(k); 0123 0124 %Error Calculation 0125 e(:,k) = d(:,k) - Y_out1(k); 0126 e_real(k)=real(d(k))-u1(k); 0127 e_imag(k)=imag(d(k))-v1(k); 0128 0129 %2nd output error calculation 0130 0131 %MSE Error Calculation 0132 E(k)=(1/2)*((e_real(k)).^2+(e_imag(k)).^2); 0133 E_dB(k)=10*log10(E(k));%error value at k step 0134 0135 %Matrix ready before calculating Pij 0136 sig_function_derR=[sig_function_der1R,sig_function_der2R]; 0137 sig_function_derI=[sig_function_der1I,sig_function_der2I]; 0138 %misinterpret formula for splitcomplex; will be back again! 0139 %Calculating Pij 0140 for l=1:N%k 0141 for t=1:N%destination 0142 for j=1:p+N+1%source 0143 tempRR=0; 0144 tempII=0; 0145 tempRI=0; 0146 tempIR=0; 0147 for r=1:N%it is typical to sum both Pij value 0148 tempRR=tempRR+((real(W(r,l))).*(PRR_old(j,t,r))-imag(W(r,l)).*(PIR_old(j,t,r))); 0149 tempII=tempII+((real(W(r,l))).*(PII_old(j,t,r))+imag(W(r,l)).*(PRI_old(j,t,r))); 0150 tempRI=tempRI+((real(W(r,l))).*(PRI_old(j,t,r))-imag(W(r,l)).*(PII_old(j,t,r))); 0151 tempIR=tempIR+((real(W(r,l))).*(PIR_old(j,t,r))+imag(W(r,l)).*(PRR_old(j,t,r))); 0152 end 0153 if l==t 0154 tempRR=tempRR+real(Uin(j)); 0155 tempII=tempII+real(Uin(j)); 0156 tempRI=tempRI-imag(Uin(j)); 0157 tempIR=tempIR+imag(Uin(j)); 0158 else 0159 tempRR=tempRR; 0160 tempII=tempII; 0161 tempRI=tempRI; 0162 tempIR=tempIR; 0163 0164 end 0165 PRR_new(j,t,l)=sig_function_derR(l).*tempRR; 0166 PII_new(j,t,l)=sig_function_derI(l).*tempII; 0167 PRI_new(j,t,l)=sig_function_derR(l).*tempRI; 0168 PIR_new(j,t,l)=sig_function_derI(l).*tempIR; 0169 0170 end 0171 end 0172 end 0173 PRR_old=PRR_new; 0174 PII_old=PII_new; 0175 PRI_old=PRI_new; 0176 PIR_old=PIR_new; 0177 0178 %weight change 0179 DW=alpha.*(e_real(k).*PRR_new(:,:,1)+e_imag(k).*PIR_new(:,:,1))+i*(e_real(k).*PRI_new(:,:,1)+e_imag(k).*PII_new(:,:,1)); 0180 w1=w1+DW(:,1); 0181 w2=w2+DW(:,2); 0182 W=[w1,w2]; 0183 0184 end%k 0185 0186 Elog=Elog+E_dB; 0187 0188 end%monte 0189 Elog=Elog/monte; 0190 0191 0192 figure(1) 0193 plot(1:len,Elog,'b-') 0194 %title('Split CRTRL Algoritm') 0195 %ylabel('Error(dB)') 0196 %xlabel('Number of iterations') 0197 %grid; 0198 %legend('CRTRL') 0199 0200 %figure(2) 0201 %plot(E) 0202 %title('Split CRTRL Algorithm') 0203 %ylabel('Error') 0204 %xlabel('Number of iterations') 0205 %grid; 0206 %legend('CRTRL') 0207 0208 0209 0210