genphaseerrorberdata.m

这是一个关于完成QAM调制的Matlab示例程序

function [] = GenPhaseErrorBerData(gammaRangeAwgn,gammaRangeRayleigh,phaseRange,saveFile)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This function plots the BER for QAM modulation, with phase bias in 
% the receiver.  
%
% gammaRangeAwgn = the Es/No range for AWGN plotting, in dB.  gammaRangeAwgn(1) = min value,
%                  gammaRangeAwgn(2) = max value, gammaRangeAwgn(3) = granularity
% gammaRangeRayleigh = the average Es/No range for Rayleigh plotting, in dB.  Analogous to 
%                      gammaRangeAwgn
% phaseRange = the phase range for both AWGN and Rayleigh plots.  phaseRange(1) = min value,
%              phaseRange(2) = max value.  phaseRange(3) = granularity.  Values in radians.
% saveFile = filename where the computed values shall be stored.  Allows for later
%            analysis without re-running simulations.  If saveFile = '', then 
%            data is not saved.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NUM_SYMB_LIMIT = 100000;
BER_LIMIT = 0.00001

gammaSamplesAwgn = floor((gammaRangeAwgn(2)-gammaRangeAwgn(1))/gammaRangeAwgn(3)) + 1;
gammaSamplesRayleigh = floor((gammaRangeRayleigh(2)-gammaRangeRayleigh(1))/gammaRangeRayleigh(3))+1;
phaseSamples = floor((phaseRange(2)-phaseRange(1))/phaseRange(3))+1;
% in the following matrix, the third index indicates if it is a theoretical or simulated value,
% the fourth index indicates if it is AWGN or Rayleigh on AWGN.  Index = 1 means the former,
% index = 2 means the latter.
ber = ones(max([gammaSamplesAwgn,gammaSamplesRayleigh]),phaseSamples,2,2);

% Compute AWGN BER values from theoretical equation
phaseIndex = 1;
for phase = phaseRange(1):phaseRange(3):phaseRange(2)
    gammaIndex = 1;
    for gamma = gammaRangeAwgn(1):gammaRangeAwgn(3):gammaRangeAwgn(2)
        ber(gammaIndex,phaseIndex,1,1) = 0.25*(erfc(sqrt(10^(gamma/10))*sin((pi/4)-phase))+ ...
            erfc(sqrt(10^(gamma/10))*sin((pi/4)+phase)));
        gammaIndex = gammaIndex + 1;
    end
    phaseIndex = phaseIndex + 1;
end    

% Compute Rayleigh BER values from theoretical equation
phaseIndex = 1;
for phase = phaseRange(1):phaseRange(3):phaseRange(2)
    gammaIndex = 1;
    for gamma = gammaRangeRayleigh(1):gammaRangeRayleigh(3):gammaRangeRayleigh(2)
        ber(gammaIndex,phaseIndex,1,2) = BerRayleighPhaseError(((10^(gamma/10))/2),phase);
        gammaIndex = gammaIndex + 1;
    end
    phaseIndex = phaseIndex + 1;
end    

% Compute AWGN BER values from simulations
phaseIndex = 1;
for phase = phaseRange(1):phaseRange(3):phaseRange(2)
    gammaIndex = 1;
    for gamma = gammaRangeAwgn(1):gammaRangeAwgn(3):gammaRangeAwgn(2)
        if (ber(gammaIndex,phaseIndex,1,1) > BER_LIMIT)
           % Use the theoretical BER for a guess of the appropriate number of symbols
           numSymb = round(100/ber(gammaIndex,phaseIndex,1,1))
            if (numSymb > NUM_SYMB_LIMIT)
                numIters = round(numSymb/NUM_SYMB_LIMIT);
                numSymb = NUM_SYMB_LIMIT;
                %break;
            else
                numIters = 1;
            end
            ber(gammaIndex,phaseIndex,2,1) = 0;
            for jj = 1:numIters
               [tempBer,ser] = basicModDemod(gamma,0,numSymb,phase);
               ber(gammaIndex,phaseIndex,2,1) = ber(gammaIndex,phaseIndex,2,1) + tempBer;
            end
            ber(gammaIndex,phaseIndex,2,1) = ber(gammaIndex,phaseIndex,2,1)/numIters;           
            %[ber(gammaIndex,phaseIndex,2,1),ser] = basicModDemod(gamma,0,numSymb,phase);
        end
        % Save data
		if (~isempty(saveFile))
            save(saveFile);
        end
        gammaIndex = gammaIndex + 1;
    end
    phaseIndex = phaseIndex + 1;
end    


% Compute Rayleigh BER values from simulations
phaseIndex = 1;
for phase = phaseRange(1):phaseRange(3):phaseRange(2)
    gammaIndex = 1;
    for gamma = gammaRangeRayleigh(1):gammaRangeRayleigh(3):gammaRangeRayleigh(2)
        if (ber(gammaIndex,phaseIndex,1,2) > BER_LIMIT)
            % Use the theoretical BER for a guess of the appropriate number of symbols
            numSymb = round(100/ber(gammaIndex,phaseIndex,1,2))
            if (numSymb > NUM_SYMB_LIMIT)
                numIters = round(numSymb/NUM_SYMB_LIMIT);
                numSymb = NUM_SYMB_LIMIT;
                %break;
            else
                numIters = 1;
            end
            ber(gammaIndex,phaseIndex,2,2) = 0;
            for jj = 1:numIters
               [tempBer,ser] = basicModDemod(gamma,1,numSymb,phase);
               ber(gammaIndex,phaseIndex,2,2) = ber(gammaIndex,phaseIndex,2,2) + tempBer;
            end
            ber(gammaIndex,phaseIndex,2,2) = ber(gammaIndex,phaseIndex,2,2)/numIters;   
        end
        % Save data
		if (~isempty(saveFile))
            save(saveFile);
        end
        gammaIndex = gammaIndex + 1;
    end
    if 0
        for gamma = gammaRangeRayleigh(1):gammaRangeRayleigh(3):gammaRangeRayleigh(2)
            numSymb = round(100/ber(gammaIndex,phaseIndex,1,2))
            if (numSymb > NUM_SYMB_LIMIT)
                break;
            end
            [ber(gammaIndex,phaseIndex,2,2),ser] = basicModDemod(gamma,1,numSymb,phase);
            gammaIndex = gammaIndex + 1;
        end
    end %0
    phaseIndex = phaseIndex + 1;
end    

% Save data
if (~isempty(saveFile))
    save(saveFile);
end

PhaseErrorBerPlotterOnly;

if 0
% Perform plotting
figure;
% Plot AWGN values
%subplot(2,1,1);
%hold on
phaseIndex = 1;
%axis([gammaRangeAwgn(1) gammaRangeAwgn(3) ...
%        max(1e-8,min(min(min(ber(1:gammaSamplesAwgn,:,1,1))),min(min(ber(1:gammaSamplesAwgn,:,2,1))))) ...
%        max(max(max(ber(1:gammaSamplesAwgn,:,1,1))),max(max(ber(1:gammaSamplesAwgn,:,2,1))))]);
for phase = phaseRange(1):phaseRange(3):phaseRange(2)
    semilogy([gammaRangeAwgn(1):gammaRangeAwgn(3):gammaRangeAwgn(2)],ber(1:gammaSamplesAwgn,phaseIndex,1,1),'b-',...
        [gammaRangeAwgn(1):gammaRangeAwgn(3):gammaRangeAwgn(2)],ber(1:gammaSamplesAwgn,phaseIndex,1,1),'bo',...
        [gammaRangeAwgn(1):gammaRangeAwgn(3):gammaRangeAwgn(2)],ber(1:gammaSamplesAwgn,phaseIndex,2,1),'c-.',...
        [gammaRangeAwgn(1):gammaRangeAwgn(3):gammaRangeAwgn(2)],ber(1:gammaSamplesAwgn,phaseIndex,2,1),'c+');
    hold on;
    phaseIndex = phaseIndex + 1;
end
%hold off
% Plot Rayleigh values
%subplot(2,1,2);
%hold on
figure
phaseIndex = 1;
for phase = phaseRange(1):phaseRange(3):phaseRange(2)
    semilogy([gammaRangeRayleigh(1):gammaRangeRayleigh(3):gammaRangeRayleigh(2)],ber(1:gammaSamplesRayleigh,phaseIndex,1,2),'b-',...
        [gammaRangeRayleigh(1):gammaRangeRayleigh(3):gammaRangeRayleigh(2)],ber(1:gammaSamplesRayleigh,phaseIndex,1,2),'bo',...
        [gammaRangeRayleigh(1):gammaRangeRayleigh(3):gammaRangeRayleigh(2)],ber(1:gammaSamplesRayleigh,phaseIndex,2,2),'c-.',...
        [gammaRangeRayleigh(1):gammaRangeRayleigh(3):gammaRangeRayleigh(2)],ber(1:gammaSamplesRayleigh,phaseIndex,2,2),'c+');
    hold on;
    phaseIndex = phaseIndex + 1;
end
%hold off
end % 0