roi.m

qam的matlab程序。m文件

%roi.m
%ring of initializations
%crj jr, revised 4/26/03

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);

n=1;     % largest delay
el=n+1;  % FIR equalizer dimension

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

%select training sequence length
p=4000;

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

schoice=0;
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=0.1;
om=1.4;
ng=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 equalizer

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')

%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;

algch=input('select algorithm: 0 - trained LMS, 1 - DD LMS, 2 - blind CMA =  ');
mu=input('enter stepsize (try 0.004) = ');

sdelt=delt; %training signal delay set to optimal
fdelt=foptmat(:,sdelt+1);
comdelt=conv(fdelt,tpars);
frcod=freqz(comdelt,1,ww);

for initind=1:8
thet(:,strt)=fdelt;
thet(1,strt)=thet(1,strt)*(1+0.3*cos((initind-1)*pi/4));
thet(2,strt)=thet(2,strt)+(0.3*sin((initind-1)*pi/4)*thet(1,strt));
yeq=zeros(size(st));
qyeq=zeros(size(yeq));

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
end

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

end
hold off

end