function[Avg_BER]=c10_MCBPSKRayFlatrun(N,M,z)
BlockSize=M;                             %block size
RealizationNumber=N;                    %Number of realizations
NoiseSigma=sqrt(1/z);                   %scale noise level
DetectedSymbols=zeros(1,BlockSize);         %initialize vector

Errors=zeros(1,RealizationNumber);                                   %initialize error counter
RxInput=zeros(RealizationNumber,BlockSize);
IntegratorOutput=zeros(RealizationNumber,BlockSize);

%Generate fading channel gain factor for each sample
   for j=1 : RealizationNumber
        gain(j)=abs(randn(1)+i*randn(1));
   end;
    
    %Generate transmitted symbols
    %
    [TxOutput,TxSymbols]=random_binary(BlockSize,1);

    
    %Generate thermal noise for each realization
    NoiseSamples=NoiseSigma*randn(N,BlockSize);  %generate ch. noise
    
    for j=1:RealizationNumber
        RxInput(j,:)=(gain(j)*TxOutput(1,:))+NoiseSamples(j,:); %Add signal and noise
    end;
    
    %Pass Received signal through matched filter
    %
    for j=1:RealizationNumber
        IntegratorOutput(j,:)=filter(gain(j),1,RxInput(j,:)); 
        for k=1:BlockSize
            DetectedSymbols(j,k)=(1-sign(IntegratorOutput(j,k)))/2;
            if(DetectedSymbols(j,k)~=TxSymbols(k))
               Errors(j)=Errors(j)+1;
            end;
        end;
    end;
    %Sample matched filter output every SamplePerSymbol sa,ples, compare
    %to transmitted bit, and count errors.
    %

BER=Errors/BlockSize;              %calculate BER
Avg_BER=sum(BER(:))/RealizationNumber;

%End