receiver_mainst.asv

该系统主要使用matlab实现了stbc2*1系统的发送和接收

%Spectrum of QPSK
%
M1=64;        %M1 = number of points per lobe.
fstep=fb/M1; %fstep = freq. step.  
Nf=floor((fend-fstart)/fstep); %Nf = total number of freq. points computed.
for i = 1:Nf
f(i)=(i-1)*fstep+fstart;
        if f(i) == fc S(i)=Ac^2*tb;
        else
         S(i)=(Ac^2*tb)*(sin(pi*(f(i)-fc)*2*tb)/(pi*(f(i)-fc)*2*tb))^2;
        end
end
for i = 1:Nf
dBS(i) = 10*log(S(i));
if dBS(i) < -200 dBS(i)=-200; end
end
figure
subplot(2,1,1)
plot(f,S)
grid
%xlabel('Frequency')
ylabel('Magnitude')
title('QPSK Spectrum (Linear)')
subplot(2,1,2)
plot(f,dBS)
grid
xlabel('Frequency')
ylabel('Magnitude')
title('QPSK Spectrum (dB)')
end

%BEGIN ROUTINE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; tic;
%Loop over the SNR range
for snr=Lsnr:Usnr,
    %Write to file to show routines' progress
    fid = fopen([fname '.txt'],'a');
    fprintf(fid,'Lower SNR: %d, Current SNR: %d, Upper SNR: %d, Elapsed Time: %6.0f \n',Lsnr,snr,Usnr,toc);
    fclose(fid);
    
    %LOCAL INITIALIZATIONS (current SNR loop)
    errShat = 0; 
    ndatabits = 0;
    t = 0;
    
    %Initialize particle filter
    xP = zeros(nRx,nTx,Npart);  %Current estimate of H matrix  
    Xhat = zeros(size(HH));
    for r = 1:ratio(2), %Over number of blocks
        H = HH(:,:,t+delay);
        xP = repmat(H,[1 1 Npart]);  %Initialize particles to truth
        xP = xP + sqrt(10^(-snr/10)*nTx/2)*(randn(size(xP)) + sqrt(-1)*randn(size(xP))); %Add noise proportional to SNR, (from equivalent PAT estimate)
        Xhat(:,:,t+1:t+delay) = repmat(H,[1 1 delay]); %Initialize state estiamte to truth
        t = t + delay;
    end
    wP = 1/Npart * ones(Npart,1);  %Initialize weights
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %START OF RECEIVER
    for f = 2:nframe,   %loop over number of frames
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %TRAINING BLOCK%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        MLest = zeros(size(H));
        for r = 1:ratio(1), %How many training blocks do I send
            %Next delay channel realizations
            Hd = repmat(HH(:,:,t),[1 1 delay]);
            Hml = zeros(nRx,nTx,Npart); %Maximum liklihood H
            
            %Simulate transmission across channel using Space-time block
            %code
            [Sv,rx,G,msg]=stbcG2_tx(snr,Hd,delay);
              
            switch TrainKind
                case 'perfect' %Assume perfect training
                    xhat = Hd(:,:,end); %Take most recent channel realizations
                    xP = repmat(xhat,[1 1 Npart]);
                    wP = 1/Npart*ones(Npart,1); %Force equal weights
                    
                case 'ML' %Maximum Liklihood Estimate
                     MLest = MLest*(r-1)/r + 1/r*rx*G'*inv(G*G');
                     xhat = MLest;
                     xP = repmat(MLest,[1 1 Npart]);
                     wP = ones(Npart,1)/Npart;
                    
                case 'PF',
                    %ML Channel estimate
                     MLest = MLest*(r-1)/r + 1/r*rx*G'*inv(G*G');
                     %Weight the ML estimate more heavily at low SNR,
                     %Otherwise PF may diverge
                     delta = snr/30; if delta>1, delta=1; end;
                     %Perform particle correction.
                     xP = delta*xP + (1-delta)*repmat(MLest,[1 1 Npart]);
                     wP = wP/sum(wP); %Force equal weights
                    
                    %PREDICT - Current state estimate
                     [xP] = model_AR(xP,DS,delay);
                    
                    %Gradient Correction
                    if grad, [xP] = predict_grad(xP,G,rx); 
                    end
                
                    %FILTER - Current `better' state estiamte 
                    [wP,xP] = Pfilter(wP,xP,snr,rx,G);
                    [xP,wP,nip] = resampleA(xP,wP,MLest);
                    [xhat]=particle_choose(xP,wP);
                                 
                    %Update AR coefficients (if necessary)
                    switch model
                        case 'DAR'
                            xold = Xhat(:,:,t);
                            [kX,kP] = model_DAR_filter(kX,kR,kRR,kQQ,xhat,xold);
                    end
            end%switch TrainKind
            Xhat(:,:,t+1:t+delay) = repmat(xhat,[1 1 delay]); %Upate Xhat
            t = t + delay;  %Increment time
        end%ratio(1)  
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        %FILTERING BLOCK%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for r = 1:diff(ratio),
            %Next delay channel realizations
            Hd = repmat(HH(:,:,t),[1 1 delay]);
            %Hd = repmat((randn(nRx,nTx)+i*randn(nRx,nTx))/sqrt(2),[1 1 delay]); %For testing only.
            %Simulate transmission
            [Sv,rx,G,msg]=stbcG2_tx(snr,Hd,delay);
            
            %Determine kind of filtering mode state to use
            switch FilterKind
                case 'perfect'
                    xhat = Hd(:,:,delay);
                    [Svest,Gest,shat]=stbcG2_rx(rx,xhat);
                    
                case 'static',
                    xhat = MLest;
                    [Svest,Gest,shat]=stbcG2_rx(rx,xhat);
                    
                case 'PF', %i.e., perform particle filtering
                    
                   %PREDICT - Current state estimate
                    [xP] = model_AR(xP,DS,delay);
                                      
                    %Estimate Symbol
                    [xhat]=particle_choose(xP,wP);
                    [Svest,Gest,shat]=stbcG2_rx(rx,xhat);
                              
                    %FILTER - Current `better' state estiamte 
                    if grad, [xP] = predict_grad(xP,Gest,rx); 
                    end
                    [wP,xP] = Pfilter(wP,xP,snr,rx,Gest);
                    [xP,wP,nip] = resampleA(xP,wP,MLest);
                    [xhat]=particle_choose(xP,wP);
                            
                    %Estimate Symbol
                    [Svest,Gest,shat]=stbcG2_rx(rx,xhat);
                                                
            end%case kind  
            %Form Decoded Symbol Vector and Estimated Channel
            errShat = errShat + sum(shat~=msg); %or Sv and Svest
            Xhat(:,:,t+1:t+delay) = repmat(xhat,[1 1 delay]); %Save state estimate
            ndatabits = ndatabits + symbperblock;
            t = t + delay;  %Increment time
        end %forward time loop   End Filtering Block
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    end%nframe
    eval( ['Hhat.dB' num2str(snr) ' = Xhat'] );
    
    %Get statistics of results of filtering
    fprintf('CPU time (seconds), %6.1f\n',toc)
    %Get SER
    SER(snr-Lsnr+1)=errShat/ndatabits; 
    %Get MSE
    mse = zeros(nRx,nTx);
    for n = 1:nRx,
        for m = 1:nTx,
            path1 = squeeze(HH(n,m,ratio(2):end));
            path2 = squeeze(Xhat(n,m,ratio(2):end));
            mse(n,m) = mean( real(path1-path2).^2 + imag(path1-path2).^2 );
        end
    end
    MSE(snr-Lsnr+1) = mean(mean(mse));
end%snr

fid = fopen([fname '.txt'],'a');
fprintf(fid,'FINISHED!!!   Total Elapsed Time: %6.0f\n',toc)    
fclose(fid);

%Display Results
[MSE SER]
fprintf('Total time, %6.0f\n',toc)
figure(1);
subplot(1,2,1)
semilogy(Lsnr:Usnr,MSE,plotopt);grid on; hold on;
title('Average MSE / Channel')
xlabel('SNR')
ylabel('MSE')
subplot(1,2,2)
semilogy(Lsnr:Usnr,SER,plotopt);grid on; hold on;
title('Symbol Error Rate')
xlabel('SNR')
ylabel('SER')
figure(2); clf;
subplot(2,1,1);
plot(squeeze(real(Xhat(1,1,:))),'r'); hold on
plot(squeeze(real(HH(1,1,:))),'k'); hold on
xlabel('Real Component')
ylabel('Voltage')
subplot(2,1,2);
plot(squeeze(imag(Xhat(1,1,:))),'r'); hold on
plot(squeeze(imag(HH(1,1,:))),'k'); hold on
xlabel('Imaginary Component')
ylabel('Voltage')
legend('PF Estimate','Channel')