%FDBDeconv  Frequency-domain blind deconvolution.
%
%   Marcel Joho, 24.11.2000 (joho@ieee.org)
%

clear
%%% parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C=512; Nw=50; Tx=250; iter=500; plot_isi=20;
Tu=Tx-Nw; Ty=Tx-3*Nw; L=2*Tu+1; L2=2*Ty+1;

%%% initial conditions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
len_s=iter*L+L;
x=zeros(C,1); u=x; y=u; u_out=zeros(L*iter,1);

Unit=fft([1;zeros(C-1,1)],C); W=Unit;

%%% a(z) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=[ -1-4i 1-5i -11-2i -17-11i -1+20i]/10;	% a(z)

Na=length(a); a=a(:)/sqrt(sum(a.*conj(a)));
a=[a;zeros(C-Na,1)]; A=fft(a,C);

%%% projection matrices %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ng=Na+Nw-1;
gg=[(C-Ng+1:C) (1:Ng+1)]';
ww=[(C-Nw+1:C) (1:Nw+1)]'; Pw=zeros(C,1); Pw(ww)=1;
xx=[(C-Tx+1:C) (1:Tx+1)]';
uu=[(C-Tu+1:C) (1:Tu+1)]'; Pu=zeros(C,1); Pu(uu)=1;
yy=[(C-Ty+1:C) (1:Ty+1)]'; Py=zeros(C,1); Py(yy)=1;

%%% optimal deconvolution filter w(z) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W_opt=1./A; w_opt=ifft(W_opt);
figure(1); plot((-Nw:Nw),[real(w_opt(ww)) imag(w_opt(ww))]);
ylabel('w_{opt}(z)'); xlabel('tap n'); drawnow ; pause(1);

%%% source signal %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%s=randn(len_s,2);s=sign(s).*s.^2;s_seq=s*[1;i]/sqrt(6);  % super-G
s_seq=(sign(randn(len_s,1))+i*sign(randn(len_s,1)))/sqrt(2);% sub-G

n_seq=0.032*(randn(len_s,1) + i*randn(len_s,1)) / sqrt(2);
x_seq=filter(a(1:Na),[1],s_seq) + n_seq;

%%% algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:iter
  t=k*L; k

  x(xx)=x_seq(t-Tx:t+Tx); X=fft(x);
  U=W.*X; u=Pu.*ifft(U); u_out((k-1)*L+1:k*L) = u(uu);

  %%% nonlinearity y=g(u) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  ur=real(u); ui=imag(u);
% y=sign(ur) + i*sign(ui);		% super-G (re/im)
% y=sign(u);				% super-G (abs)
% y=ur.*abs(ur).^2 + i*ui.*abs(ui).^2;	%   sub-G (re/im)
  y=u.*abs(u).^2;			%   sub-G (abs)

  eb=Pu.*(u-y); Eb=fft(eb);		% blind error
  y=Py.*y; Y=fft(y);

  %%% update equations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
							% EASI
% W  = W + 0.05/L2 * (L2*Unit-U.*conj(U)+U.*conj(Y)-Y.*conj(U)).*W;
  W  = W + 0.05 * (Unit - Y.*conj(U)/L2).*W;		% nat.grad.I
% W  = W + 0.05 * (W - Y.*conj(U).*W/L2);		% nat.grad.II
% W  = W + 0.05/L2 * (L2*conj(W).^(-1) - Y.*conj(X));	% Infomax

% mu = 0.3/L; mu_2 = mu./(Unit-mu*conj(U).*Eb); 	% BLMS-1Wx
% W  = W + mu_2.*Eb.*conj(U).*W; 			% BLMS-1Wx
% W  = W + 0.3/L*Eb.*conj(U).*W; 	 		% BLMS-1Wx-c
% mu = 0.05; mu_2 = mu./((1-mu)*Unit+mu/L*conj(U).*U);	% BLMS-2Wx-a
% W  = 1/(1-mu) * (Unit - mu_2/L2.*Y.*conj(U)).*W; 	% BLMS-2Wx-a
% mu = 0.05; W = 1./((1-mu)*Unit + mu/L*Y.*conj(U)).*W;	% BLMS-2Wx-b
% W  = W + 0.05/L * Eb.*conj(X);  			% BLMS-3Ws
% W  = W + 0.05/L2 * (L2*conj(W).^(-1) - Y.*conj(X));	% BLMS-4Ws

  %%% filter projection operation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  w=Pw.*ifft(W); W=fft(w);

  %%% performance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  G=W.*A; g=ifft(G);
  Jisi(k)=sum(g.*conj(g))/max(g.*conj(g))-1;

  %%% plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if rem(k,plot_isi)==0,
    figure(1); plot(10*log10(Jisi));
               ylabel('J_{ISI}  [dB]');xlabel('k')
    figure(2); plot((-Nw:Nw),[real(w(ww)) imag(w(ww))]);
               ylabel('w(z)');xlabel('tap n')
    figure(3); plot((-Ng:Ng),[real(g(gg)) imag(g(gg))]);
               ylabel('g(z)');xlabel('tap n')
    figure(4); plot(real(u(uu)),imag(u(uu)),'.');
    drawnow
  end
end
