oae.m

qam的matlab程序。m文件

%oae.m
%optimal and adaptive equalization
%crj jr, revised 6/29/01

clear all
close all


%FIR channel impulse response vector
rawtp=[1 0.5 0.3].';   
%rawtp=[1 1 1 0.2 -0.4 2 -1].';
%rawtp=[-0.2 .1 .3 1 1.2 .4 -.3 -.2 .3 .1 -.1].';   

chdim=length(rawtp);

%select FIR equalizer dimension
n=input('enter largest equalizer delay (try 3)= ');
el=n+1;

%maximum delay candidate (> ch dim + eq max delay);
alph=chdim+n+1;

%select training sequence length
p=input('enter training/source sequence length (try 4000) = ');

%set up time vector
time=[1:p];

schoice=input('enter 0 for binary or 1 for 4-PAM = ');
if schoice <0.5
%set up white, zero-mean, BPSK source
s=sign(rand(size(time))-.5);
dc=1;   %dispersion constant
else
%set up white, zero-mean, PAM source
yo=(rand(size(time))-.5);
dc=41/25;   %dispersion constant
for indx=1:p
if abs(yo(indx)) > 0.25
s(indx)=(3/sqrt(5))*sign(yo(indx));
else
s(indx)=(1/sqrt(5))*sign(yo(indx));
end
end
end
st=s';


%make Toeplitz input matrix for channel
toprowc=fliplr(s(1:chdim));
lcolc=s(chdim:p).';
Ach=toeplitz(lcolc,toprowc);

%create received signal r(chdim) to r(p)
%NOTE: rsig (and r and d) index j to time index i via j=i-chdim+1 !!! 
%NOTE: To use this r in forming Rbar alph-n needs to be > chdim !!!
sg=input('enter sinusoidal interferer gain (try 0.1) = ');
om=input('enter sinusoidal interferer frequency (try 1.4) = ');
ng=input('enter uniform channel noise -/+ max magnitude (try 0.04) ='); 
tpars=sqrt((1-ng*(ng/3)-0.5*sg*sg)/(rawtp'*rawtp))*rawtp;
rsig=Ach*tpars;
d=zeros(size(rsig));
for i=1:length(rsig),
d(i,1)=sg*sin(om*(i-1))+ng*2*(rand-0.5);
end
r=(rsig+d)';

%calculate MMSE equqlizer

cosvec=cos(om*[0:(el-1)])'; 
V=toeplitz(cosvec);   %interference matrix

isnr=(ng^2)/3;  %inverse of source to noise variance ratio 

C=convmtx(tpars,el);  %channel convolution matrix

CC=C'*C+isnr*eye(el)+0.5*(sg^2)*V;  
     %convlution matrix transposed times itself + isnr times I

foptmat=pinv(CC)*C';   %pseudoinverse of (C'C+ isnr*I)

[mmse,optind]=min(diag(eye(size(C*foptmat))-C*foptmat)); 
     %value and index of min entry of I - C*foptmat
mmse_w_opt_delay=mmse
delt=optind-1;
optimum_delay=delt


fopt=foptmat(:,optind);   
     %optimum equalizer over equalizer parameters 
     %and possible delays

combo=conv(fopt,tpars); %channel and MMSE equalizer combination

%create MMSE equalizer output
%compose Rbar
toprowrbar=fliplr(r(alph-chdim-n+2:alph-chdim+2));
lcolrbar=r(alph-chdim+2:p-chdim+1);
Rbar=toeplitz(lcolrbar,toprowrbar);
y=Rbar*fopt;

figure (1)
subplot(2,2,1)
plot([1:length(r)],r,'.')
title('received signal')
subplot(2,2,2)
plot([1:length(y)],y,'.')
title('optimal equalizer output')
subplot(2,2,3)
plot([0:length(combo)-1],combo,'o')
title('combined chan and opt eq imp resp')
subplot(2,2,4)
qy=zeros(size(y));
if schoice <0.5
qy=sign(y);
else
for indd=1:length(y)
if abs(y(indd))>2/sqrt(5)
qy(indd)=sign(y(indd))*3/sqrt(5);
else
qy(indd)=sign(y(indd))/sqrt(5);
end
end
end
plot([1:length(y)],qy.'-s(alph-delt+1:p-delt),'.')
title('decision device recovery error')
percent_dec_errs=100*sum(abs(sign(qy.'-s(alph-delt+1:p-delt))))/length(y)
figure(2)
subplot(2,2,1)
axis equal
axis square
zplane(roots(tpars),[])
title('fir channel zeros')
subplot(2,2,3)
axis equal
axis square
zplane(roots(combo),[])
title('chan-optimum eq combo zeros')
ww=[0:pi/100:pi];
[frtp]=freqz(tpars,1,ww);
[frfo]=freqz(fopt,1,ww);
[frco]=freqz(combo,1,ww);
subplot(2,2,2)
axis normal 
plot(ww,20*log10(abs(frtp)),'.',ww,20*log10(abs(frfo)),'--',...
                                ww,20*log10(abs(frco)),'-')
title('Freq Resp Magnitude') %chan: dot, eq: dash, combo:solid
ylabel('db')
xlabel('radians')
subplot(2,2,4)
plot(ww,unwrap(angle(frtp)),'.',ww,unwrap(angle(frfo)),'--',...
                                ww,unwrap(angle(frco)),'-')
title('Freq Resp Phase') %chan: dot, eq: dash, combo:solid
ylabel('radians')
xlabel('radians')

adapt=input('enter 0 to stop or 1 to adaptively equalize = ');

%adaptive equalizer simulation loop

itno=p-length(tpars)-length(fopt);
thet=zeros(n+1,itno);
sspe=zeros(1,itno);
e=zeros(size(st));
strt=1+max(delt,length(tpars))+n;
if adapt >0.5
sdelt=input('enter training signal delay (try optimal) = ');
df=input('enter initial eq disp factor about mmse soln (try 0.5) = ');
fdelt=foptmat(:,sdelt+1);
comdelt=conv(fdelt,tpars);
frcod=freqz(comdelt,1,ww);
thet(:,strt)=fdelt-df*2*(rand(size(thet(:,strt)))-0.5*ones(size(thet(:,strt))));
sspe(strt)=(fdelt-thet(:,strt))'*(fdelt-thet(:,strt));
else
end
yeq=zeros(size(st));
qyeq=zeros(size(yeq));
cnt=0;
while adapt >0.5
cnt=cnt+1;
algch=input('select algorithm: 0 - trained LMS, 1 - DD LMS, 2 - blind CMA =  ');

mu=input('enter stepsize (try 0.004) = ');

for i=strt:itno-1,
yeq(i)=thet(:,i)'*flipud(r(i-n-chdim+1:i-chdim+1)');
if schoice<0.5
qyeq(i)=sign(yeq(i));
else
if abs(yeq(i))>2/sqrt(5)
qyeq(i)=sign(yeq(i))*3/sqrt(5);
else
qyeq(i)=sign(yeq(i))/sqrt(5);
end
end

e(i)=st(i-sdelt,1)-yeq(i);  
if algch < 0.5 
thet(:,i+1)=thet(:,i)+mu*e(i)*flipud(r(i-n-chdim+1:i-chdim+1)'); %LMS
elseif algch < 1.5
thet(:,i+1)=thet(:,i)+mu*(qyeq(i)-yeq(i))...
                     *flipud(r(i-n-chdim+1:i-chdim+1)'); %DD LMS
else
thet(:,i+1)=thet(:,i)+mu*yeq(i)*(dc-yeq(i)*yeq(i))...
                               *flipud(r(i-n-chdim+1:i-chdim+1)'); %CMA
end
sspe(i+1)=(fdelt-thet(:,i+1))'*(fdelt-thet(:,i+1));
end

%plot summed squared parameter errors
figure(3+cnt)
subplot(2,2,1)
axis('normal')
semilogy([1:itno-strt+1],sspe(strt:itno),'.')
xlabel('iterations')
ylabel('summed squared parameter error')

%plot squared prediction errors
subplot(2,2,2)
axis('normal')
plot([1:itno-strt],e(strt:itno-1).*e(strt:itno-1),'.')
xlabel('iterations')
ylabel('squared prediction error')

%plot equalizer output
subplot(2,2,3)
axis('normal')
plot([1:itno-strt],yeq(strt:itno-1),'.')
xlabel('iterations')
ylabel('adaptive equalizer output')

%plot combined channel and first and last adapted equalizer frequency response
subplot(2,2,4)
axis('normal')
ww=[0:pi/100:pi];
[fradin]=freqz(conv(tpars,thet(:,strt)),1,ww);
[fras]=freqz(conv(tpars,thet(:,itno)),1,ww);
plot(ww,20*log10(abs(fras)),'.', ww,20*log10(abs(frcod)),'-',...
                                 ww,20*log10(abs(fradin)),'--')
title('Combo Freq Resp Magnitude') %adap fin: dot, adap init: dash, opt:solid
ylabel('db')
xlabel('radians')

figure(3)
hold on
axis equal
plot(thet(1,strt:itno-1),thet(2,strt:itno-1))
plot(thet(1,strt),thet(2,strt),'o')
plot(thet(1,itno),thet(2,itno),'s')
xlabel('f0')
ylabel('f1')
hold off

adapt=input('enter 0 to stop or 1 to adaptively equalize = ');
end