%Matlab program for use in conjunction with Homework no. 6
%on adaptive equalization for digital communications.
clear all
clf
set(0,'defaultaxesfontsize',22);
Nb=128;
%
%generate random QPSK information symbols
%
br=ones(1,Nb);
temp=rand(1,Nb);
br(find(temp<.5))=-1;
bi=ones(1,Nb);
temp=rand(1,Nb);
bi(find(temp<.5))=-1;
a=br+j*bi;
%
%create multipath distorted waveform -- channel impulse response q[n]:
%
M1=5; M2=5;
N2=12;
Md=N2-M2;
ph=pi-pi/8;
gc=.707*exp(j*ph);
nc=-M1:M2;
beta=0.35;
noff=nc-0.5;
npi=noff*pi;
ndenom=1-(2*beta*noff).^2;
q=gc*(sin(npi)./npi).*(cos(beta*npi)./ndenom);
q(1,M1+1)=q(1,M1+1)+1;
%
%create multipath distorted data, x[n]:
%
x=zeros(1,Nb+N2+M1);
for n=(M1+Md):(Nb+M1+Md-1);
x(1,n-M1:n+M2)=x(1,n-M1:n+M2)+a(1,n-M1-Md+1)*q;
end
x=x+0.1*(randn(size(x))+j*randn(size(x)));
%
%estimation of optimal filter coefficient vector
%based on block time-averaged estimates of Rxx and rdx
%
M=M1+1+N2;
Rxx=zeros(M,M);
rdx=zeros(M,1);
w=1; alpha=1/w;
for n=(M1+Md+2):(Nb+M1+Md);
Rxx=Rxx+x(1,n+M1:-1:n-N2).'*conj(x(1,n+M1:-1:n-N2));
rdx=rdx+conj(a(1,n-M1-Md))*x(1,n+M1:-1:n-N2).';
end
%
%initial filter coefficient vectors for LMS and RLS
%
hblock=inv(Rxx)*rdx;
hlms=[zeros(M1,1) ; 1; zeros(N2,1)];
hrls=hlms;
rxx0=x*x'/(M1+M2+Nb);
Rxx=Rxx/Nb;
mu=2/(2*trace(Rxx));
Rinv=0.1*eye(M,M);
%
%real-time simulation of LMS and RLS adaptations
%iterating over successive data blocks of length M
%
for n=(M1+Md+1):(Nb+M1+Md);
y(1,n-M1-Md)=x(1,n+M1:-1:n-N2)*conj(hblock);
%LMS adaptation
ylms(1,n-M1-Md)=x(1,n+M1:-1:n-N2)*conj(hlms);
errlms=a(1,n-Md-M1)-ylms(1,n-M1-Md);
hlmsold=hlms;
hlms=hlmsold+mu*conj(errlms)*x(1,n+M1:-1:n-N2).';
%RLS adaptation
hrlsold=hrls;
errrls=a(1,n-Md-M1)-x(1,n+M1:-1:n-N2)*conj(hrlsold);
Rinvold=Rinv;
f=Rinvold*x(1,n+M1:-1:n-N2).';
murls=x(1,n+M1:-1:n-N2)*conj(f);
Kvec=f/(w+murls);
Rinv=alpha*(Rinvold-Kvec*f');
hrls=hrlsold+conj(errrls)*Kvec;
yrls(1,n-M1-Md)=x(1,n+M1:-1:n-N2)*conj(hrls);
end
%
%Plot output signal constellations
%
Nl=32;
plot(real(x(1,Md+M1:Nb+Md+M1)),imag(x(1,Md+M1:Nb+Md+M1)),'r+','MarkerSize',18,'Linewidth',3);
axis([-2 2 -2 2]);
legend('Distored Bit Values');
title('Received Signal Constellation');
pause
plot(real(y),imag(y),'kx','MarkerSize',18,'Linewidth',3);
axis([-2 2 -2 2]);
legend('Block Estimate');
title('Block Estimated Signal Constellation');
pause
plot(real(yrls(1,Nl:Nb)),imag(yrls(1,Nl:Nb)),'mo','MarkerSize',18,'Linewidth',3);
axis([-2 2 -2 2]);
legend('RLS Estimate');
title('RLS Estimated Signal Constellation');
pause
plot(real(ylms(1,Nl:Nb)),imag(ylms(1,Nl:Nb)),'r+','MarkerSize',18,'Linewidth',3);
axis([-2 2 -2 2]);
legend('LMS Estimate');
title('LMS Estimated Signal Constellation');
pause
%
%Plot frequency response of channel, inverse filter,
%and overall frequency response
%
freqomega=linspace(-pi,pi,256);
Hmag=20*log10(abs(fftshift(fft(q,256))));
%H0=max(hmag);
%Hmag=Hmag-H0;
Imag=20*log10(abs(fftshift(fft(conj(hblock),256))));
%Imag=Imag-H0;
Pmax=max(Imag);
plot(freqomega,Hmag,'k','LineWidth',4);
axis([-pi pi -30 Pmax]);
hold
plot(freqomega,Imag,'r','LineWidth',4);
g=conv(q,conj(hblock));
Pmag=20*log10(abs(fftshift(fft(g,256))));
plot(freqomega,Pmag,'b','LineWidth',4);
legend('Channel','Inverse Channel','Convolution')
title('DTFT of Channel and Inverse'),...
xlabel('Omega'), ylabel('Magnitude (dB)'), grid