%Matlab program for use in conjunction with Homework 3
%on two-channel 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=5;
Md=N2-M2;
nc=-M1:M2;
beta=0.35;
noff=nc-0.5;
npi=noff*pi;
ndenom=1-(2*beta*noff).^2;
polemag=1;
ph1=pi-pi/16;
gc1=polemag*exp(j*ph1);
q1=gc1*(sin(npi)./npi).*(cos(beta*npi)./ndenom);
q1(1,M1+1)=q1(1,M1+1)+1;
ph2=pi-pi/8;
gc2=polemag*exp(-j*ph1);
q2=gc2*(sin(npi)./npi).*(cos(beta*npi)./ndenom);
q2(1,M1+1)=q2(1,M1+1)+1;
%
%create multipath distorted data, x[n]:
%
x1=zeros(1,Nb+N2+M1); x2=zeros(1,Nb+N2+M1);
for n=(M1+Md+1):(Nb+M1+Md-1);
x1(1,n-M1:n+M2)=x1(1,n-M1:n+M2)+a(1,n-M1-Md+1)*q1;
x2(1,n-M1:n+M2)=x2(1,n-M1:n+M2)+a(1,n-M1-Md+1)*q2;
end
x1=x1+0.1*(randn(size(x1))+j*randn(size(x1)));
x2=x2+0.1*(randn(size(x2))+j*randn(size(x2)));
%
%estimation of optimal filter coefficient vector
%based on block time-averaged estimates of Rxx and rdx
%
M=M1+1+N2;
Rxx=zeros(2*M,2*M);
rdx=zeros(2*M,1);
w=1; alpha=1/w;
for n=(M1+Md+2):(Nb+M1+Md);
x=[x1(1,n+M1:-1:n-N2).' ; x2(1,n+M1:-1:n-N2).'];
Rxx=Rxx+x*x';
rdx=rdx+conj(a(1,n-M1-Md))*x;
end
%
%initial filter coefficient vectors for LMS and RLS
%
hblock=inv(Rxx)*rdx;
hlms1=[zeros(M1,1) ; 1; zeros(N2,1)];
hlms2=[zeros(M1,1) ; 1; zeros(N2,1)];
hrls=[hlms1 ; hlms2];
hrls1=[zeros(M1,1) ; 1; zeros(N2,1)];
hrls2=[zeros(M1,1) ; 1; zeros(N2,1)];
rxx0=x*x'/(M1+M2+Nb);
Rxx=Rxx/Nb;
mu=2/(2*trace(Rxx));
Rinv=0.1*eye(2*M,2*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)=[x1(1,n+M1:-1:n-N2) x2(1,n+M1:-1:n-N2)]*conj(hblock);
%LMS adaptation
ylms(1,n-M1-Md)=x1(1,n+M1:-1:n-N2)*conj(hlms1)+ x2(1,n+M1:-1:n-N2)*conj(hlms2);
errlms=a(1,n-Md-M1)-ylms(1,n-M1-Md);
hlmsold1=hlms1;  hlmsold2=hlms2;
hlms1=hlmsold1+mu*conj(errlms)*x1(1,n+M1:-1:n-N2).';
hlms2=hlmsold2+mu*conj(errlms)*x2(1,n+M1:-1:n-N2).';
%RLS adaptation
hrlsold1=hrls1;  hrlsold2=hrls2;
errrls=a(1,n-Md-M1)-x1(1,n+M1:-1:n-N2)*conj(hrlsold1)-x2(1,n+M1:-1:n-N2)*conj(hrlsold2);
Rinvold=Rinv;
f=Rinvold*[x1(1,n+M1:-1:n-N2).' ; x2(1,n+M1:-1:n-N2).'];
murls=[x1(1,n+M1:-1:n-N2) x2(1,n+M1:-1:n-N2)]*conj(f);
Kvec=f/(w+murls);
Rinv=alpha*(Rinvold-Kvec*f');
hrls=[hrlsold1 ; hrlsold2]+conj(errrls)*Kvec;
hrls1=hrls(1:M);  hrls2=hrls(M+1:2*M);
yrls(1,n-M1-Md)=x1(1,n+M1:-1:n-N2)*conj(hrls1)+ x2(1,n+M1:-1:n-N2)*conj(hrls2);
end
%
%Plot output signal constellations
%
Nl=32;
subplot(1,2,1)
plot(real(x1(1,Md+M1:Nb+Md+M1)),imag(x1(1,Md+M1:Nb+Md+M1)),'r+','MarkerSize',18,'Linewidth',3);
axis([-2 2 -2 2]);
title('Constellation 1');
subplot(1,2,2)
plot(real(x2(1,Md+M1:Nb+Md+M1)),imag(x2(1,Md+M1:Nb+Md+M1)),'r+','MarkerSize',18,'Linewidth',3);
axis([-2 2 -2 2]);
title('Constellation 2');
pause
clf
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);
Hmag1=20*log10(abs(fftshift(fft(q1,256))));
Hmag2=20*log10(abs(fftshift(fft(q2,256))));
hblock1=hblock(1:M);  hblock2=hblock(M+1:2*M);
Gmag1=20*log10(abs(fftshift(fft(hblock1,256))));
Gmag2=20*log10(abs(fftshift(fft(hblock2,256))));
Pmax=max([max(Hmag1) max(Hmag2) max(Gmag1) max(Gmag2)]);
plot(freqomega,Hmag1,'k','LineWidth',4);
axis([-pi pi -30 Pmax]);
hold on
plot(freqomega,Hmag2,'r','LineWidth',4);
%plot(freqomega,Gmag1,'m','LineWidth',4);
%plot(freqomega,Gmag2,'g','LineWidth',4);
g=conv(q1,conj(hblock1))+conv(q2,conj(hblock2));
Pmag=20*log10(abs(fftshift(fft(g,256))));
plot(freqomega,Pmag,'b','LineWidth',4);
%legend('Channel 1','Channel 2','G1(omega)', 'G2(omega)','Overall System')
legend('Channel 1','Channel 2','Overall System')
title('Frequency Responses'),...
xlabel('Omega'), ylabel('Magnitude (dB)'), grid
hold off