%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