QAM_Matlab.rar_QAM璋冨埗_qam_qam调制 matlab_srrc.m

%============================================= % Modem Receiver Simulator/ SER Determination %--------------------------------------------- % John Walsh %--------------------------------------------- % This software and manual is based upon % work done by John MacLaren Walsh, Katie % Orlicki, Adam Pierce, Johnson Smith, Andy % Klein, and Sean Leventhal under the guidance % of Dr. C. Richard Johnson, Jr. %--------------------------------------------- % Version 1.0 % 1/20/2003 function [theoSER,uncSER,codBER,retMSE]=qpskRx(muStruct,recSigStruct) OverSamp=4; %number of samples in between each symbol ApproxPeriods=6; %1/2 number of periods to Approximate the SRRC with Tb=5*10^-3; %baud interval in seconds Ts=Tb/OverSamp; %sampling period fc=200; %carrier frequency (Hz) recSig=recSigStruct.recSig; training=recSigStruct.training; msg1=recSigStruct.msg1; msg2=recSigStruct.msg2; constl=recSigStruct.constl; msg1sc=recSigStruct.msg1sc; msg2sc=recSigStruct.msg2sc; Costas_Mu1=muStruct.Costas_Mu1; Costas_Mu2=muStruct.Costas_Mu2; Time_mu1=muStruct.Time_mu1; Time_mu2=muStruct.Time_mu2; EqMu=muStruct.EqMu; ddEqMu=muStruct.ddEqMu; debugFlag=muStruct.debugFlag; codingFlag=muStruct.codingFlag; scrambleFlag=muStruct.scrambleFlag; beta=muStruct.beta; timeRecFlag=1&debugFlag; %draw the timing parameter and constellation eqFlag=1&debugFlag; scramble_flag=1&debugFlag; mseFlag=1&debugFlag; costasFlag=1&debugFlag; if scrambleFlag==1 scrambler=recSigStruct.scrambler; scramLen=recSigStruct.scramLen; end ginv=[1 1;1 0;0 0;1 0;0 1]; syn=[0 0 0 0 0; 0 0 0 0 1; 0 0 0 1 0; 0 1 0 0 0; 0 0 1 0 0; 1 0 0 0 0; 1 1 0 0 0; 1 0 0 1 0]; G=[1 0 1 0 1; 0 1 0 1 1]; H=[1 0 1 0 0; 0 1 0 1 0; 1 1 0 0 1]; %begin reciever %first we attempt downconversion %need a lpf for the costas loop fl=200; fbe=2*[0 fc*Ts-1/(16*OverSamp) fc*Ts+1/(16*OverSamp) .5]; amps=[1 1 0 0]; b=firls(fl,fbe,amps); %cmu1=Costas_Mu1/((3/16)*mean(real(pshaped).^2)^2-(1/16)* ... % mean(real(pshaped).^4)); %cmu2=Costas_Mu2/((3/16)*mean(real(pshaped).^2)^2-(1/16)* ... % mean(real(pshaped).^4)); cmu1=Costas_Mu1/.0017; cmu2=Costas_Mu2/.0017; t=[1:length(recSig)]*Ts; theta=zeros(1,length(recSig)); theta2=zeros(1,length(recSig)); v1=zeros(1,length(recSig)); v2=zeros(1,length(recSig)); v3=zeros(1,length(recSig)); v4=zeros(1,length(recSig)); v1=filter(b,1,recSig.*sin(2*pi*fc*t)); v2=filter(b,1,recSig.*sin(2*pi*fc*t+pi/4)); v3=filter(b,1,recSig.*cos(2*pi*fc*t)); v4=filter(b,1,recSig.*cos(2*pi*fc*t+pi/4)); Bdelay=(length(b)-1)/2; z1=zeros(1); z2=zeros(1); z3=zeros(1); z4=zeros(1); for k=1:length(recSig)-1-length(b)/2 cosThetaK=cos(theta(k)); sinThetaK=sin(theta(k)); z1=v1(k+Bdelay)*cosThetaK+v3(k+Bdelay)*sinThetaK; z2=v2(k+Bdelay)*cosThetaK+v4(k+Bdelay)*sinThetaK; z3=v3(k+Bdelay)*cosThetaK-v1(k+Bdelay)*sinThetaK; z4=v4(k+Bdelay)*cosThetaK-v2(k+Bdelay)*sinThetaK; theta(k+1)=theta(k)+cmu1*(z2*z4*z3*z1); cosTheta12=cos(theta(k)+theta2(k)); sinTheta12=sin(theta(k)+theta2(k)); z1=v1(k+Bdelay)*cosTheta12+v3(k+Bdelay)*sinTheta12; z2=v2(k+Bdelay)*cosTheta12+v4(k+Bdelay)*sinTheta12; z3=v3(k+Bdelay)*cosTheta12-v1(k+Bdelay)*sinTheta12; z4=v4(k+Bdelay)*cosTheta12-v2(k+Bdelay)*sinTheta12; theta2(k+1)=theta2(k)+cmu2*(z4*z2*z3*z1); end downCSig=zeros(1,length(recSig)); downCSig=2*recSig.*exp(-j*(2*pi*fc*t+theta+theta2)); if costasFlag==1 figure; [Pxx,F]=psd(downCSig,512,1/Ts); plot(F,Pxx); title('Power Spectral Density of Signal at Output of Costas Loop'); xlabel('Frequency (Hx)'); ylabel('PSD'); figure; subplot(2,1,2); plot(theta2); title('Costas Loop \theta_2'); xlabel('Time'); subplot(2,1,1); plot(theta); title('Costas Loop \theta'); ylabel('Time'); drawnow; end %now we must low pass filter the down-converted signal downSig=filter(b,1,downCSig); %timing recovery stuff=(ApproxPeriods*OverSamp+1):OverSamp: ... (length(downSig)-ApproxPeriods*OverSamp); tau=zeros(1,length(stuff)); tau2=zeros(1,length(stuff)); timRecSig=zeros(1,length(stuff)); clear stuff; i=1; delta=.1; timRecPDelta=2*ApproxPeriods*OverSamp+1; timRecMDelta=2*ApproxPeriods*OverSamp+1; timRecPDelta2=2*ApproxPeriods*OverSamp+1; timRecMDelta2=2*ApproxPeriods*OverSamp+1; for k=(ApproxPeriods*OverSamp+1):OverSamp:(length(downSig)-ApproxPeriods*OverSamp) data=downSig((k-ApproxPeriods*OverSamp):(k+ApproxPeriods* ... OverSamp)); timRecSig(i)=srrc(ApproxPeriods,beta,OverSamp,tau(i)+tau2(i))*data.'; tempSig=srrc(ApproxPeriods,beta,OverSamp,tau(i))*data.'; timRecPDelta=srrc(ApproxPeriods,beta,OverSamp,tau(i)+delta)* ... data.'; timRecMDelta=srrc(ApproxPeriods,beta,OverSamp,tau(i)-delta)* ... data.'; timRecPDelta2=srrc(ApproxPeriods,beta,OverSamp,tau(i)+tau2(i)+delta)* ... data.'; timRecMDelta2=srrc(ApproxPeriods,beta,OverSamp,tau(i)+tau2(i)-delta)* ... data.'; dXdTau2=(abs(timRecPDelta2)-abs(timRecMDelta2))/delta; tau2(i+1)=tau2(i)-Time_mu2*abs(timRecSig(i))^3*dXdTau2; dXdTau=(abs(timRecPDelta)-abs(timRecMDelta))/delta; tau(i+1)=tau(i)-Time_mu1*abs(tempSig)^3*dXdTau; i=i+1; end if timeRecFlag==1 figure; subplot(2,1,1); plot(tau+tau2); title('Timing Offset'); subplot(2,1,2); plot(timRecSig((length(timRecSig)-200):end),'b.'); ... title('Convergent Eye Diagram'); axis square; drawnow; end %equalizer %first we have to find the training rrCorr=zeros(1,length(timRecSig)-length(training)+1); irCorr=zeros(1,length(timRecSig)-length(training)+1); for k=1:length(timRecSig)-length(training)+1 rrCorr(k)=sum(real(training).*real(timRecSig(k:(k+length(training)-1)))); irCorr(k)=sum(imag(training).*real(timRecSig(k:(k+length(training)-1)))); end eqLen=6; f=zeros(eqLen,1); f(eqLen/2)=1; delta=eqLen/2-1; %default delay if max(abs(irCorr))>max(abs(rrCorr)) [y,i]=max(abs(irCorr)); sigToEq=zeros(1,length(timRecSig(i:end))); if abs(min(irCorr))>abs(max(irCorr)) if length(timRecSig) > (i+length(training)+length(msg1sc)) sigToEq=-j*timRecSig(i:end); else sigToEq=-j*timRecSig((length(timRecSig)-length(training)- ... length(msg1sc)-delta-eqLen):end); end else if length(timRecSig) > (i+length(training)+length(msg1sc)) sigToEq=j*timRecSig(i:end); else sigToEq=j*timRecSig((length(timRecSig)-length(training)- ... length(msg1sc)-delta-eqLen):end); end end else [y,i]=max(abs(rrCorr)); sigToEq=zeros(1,length(timRecSig(i:end))); if abs(min(rrCorr))>abs(max(rrCorr)) if length(timRecSig) > (i+length(training)+length(msg1sc)) sigToEq=-1*timRecSig(i:end); else sigToEq=-1*timRecSig((length(timRecSig)-length(training)- ... length(msg1sc)-delta-eqLen):end); end else if length(timRecSig) > (i+length(training)+length(msg1sc)) sigToEq=timRecSig(i:end); else sigToEq=timRecSig((length(timRecSig)-length(training)- ... length(msg1sc)-delta-eqLen):end); end end end %now that we have found the training (and done some derotation) %let's equalize error=zeros(1,length(sigToEq)); eqOut=zeros(1,length(sigToEq)); if eqFlag==1 figure; subplot(2,1,1); h1=plot(f); subplot(2,1,2); h2=plot([1+j -1-j],'b.'); axis square; axis([-1.5 1.5 -1.5 1.5]); set(h1,'EraseMode','xor'); set(h2,'EraseMode','xor'); end for k=(eqLen+1):length(training) rr=sigToEq(k:-1:(k-eqLen+1)); error(k)=training(k-delta)-rr*f; f=f+EqMu*error(k)*rr'; eqOut(k)=rr*f; if (k>20) & eqFlag set(h2,'YData',imag(eqOut(k:-1:k-20)),'XData',real(eqOut(k:-1:k-20))); set(h1,'YData',f); end drawnow; end for k=length(training):length(sigToEq) rr=sigToEq(k:-1:(k-eqLen+1)); error(k)=(sqrt(2)/2*sign(real(rr*f))+j*sign(imag(rr*f))*sqrt(2)/2)-rr*f; f=f+ddEqMu*error(k)*rr'; eqOut(k)=rr*f; end if eqFlag==1 figure; subplot(2,1,1); plot(abs(error).^2); title('Erro