How can I simulate OFDM with Rayleigh fading by using Communications System Toolbox?

I would like to adapt the example code in the following link such that the AWGN channel is preceded by a Rayleigh fading channel.

https://www.mathworks.com/help/comm/gs/qpsk-and-ofdm-with-matlab-system-objects-1.html

My code gives a bit-error rate of approximately 0.5, regardless of the signal-to-noise ratio. Could you please tell me what I am doing wrong?

M = 4;                 % Modulation alphabet
k = log2(M);           % Bits/symbol
numSC = 128;           % Number of OFDM subcarriers
cpLen = 32;            % OFDM cyclic prefix length
maxBitErrors = 1e2;    % Maximum number of bit errors
maxNumBits = 1e7;      % Maximum number of bits transmitted

qpskMod = comm.QPSKModulator('BitInput',true);
qpskDemod = comm.QPSKDemodulator('BitOutput',true);

ofdmMod = comm.OFDMModulator('FFTLength',numSC,'CyclicPrefixLength',cpLen);
ofdmDemod = comm.OFDMDemodulator('FFTLength',numSC,'CyclicPrefixLength',cpLen);

rayleighChannel = comm.RayleighChannel;

awgnChannel = comm.AWGNChannel('NoiseMethod','Variance', 'VarianceSource','Input port');

errorRate = comm.ErrorRate('ResetInputPort',true);

ofdmDims = info(ofdmMod);

numDC = ofdmDims.DataInputSize(1);

frameSize = [k*numDC 1];

EbNoVec = (0:20)';
snrVec = EbNoVec + 10*log10(k) + 10*log10(numDC/numSC);

berVec = zeros(length(EbNoVec),3);
errorStats = zeros(1,3);

for m = 1:length(EbNoVec)
    snr = snrVec(m);

      while errorStats(2) <= maxBitErrors && errorStats(3) <= maxNumBits
          dataIn = randi([0,1],frameSize);              % Generate binary data
          qpskTx = qpskMod(dataIn);                     % Apply QPSK modulation
          txSig = ofdmMod(qpskTx);                      % Apply OFDM modulation
          fadedSig = rayleighChannel(txSig);
          powerDB = 10*log10(var(fadedSig));               % Calculate Tx signal power
          noiseVar = 10.^(0.1*(powerDB-snr));           % Calculate the noise variance
          rxSig = awgnChannel(fadedSig,noiseVar);          % Pass the signal through a noisy channel
          qpskRx = ofdmDemod(rxSig);                    % Apply OFDM demodulation
          dataOut = qpskDemod(qpskRx);                  % Apply QPSK demodulation
          errorStats = errorRate(dataIn,dataOut,0);     % Collect error statistics
      end

      berVec(m,:) = errorStats;                         % Save BER data
      errorStats = errorRate(dataIn,dataOut,1);         % Reset the error rate calculator
  end

berTheory = berawgn(EbNoVec,'psk',M,'nondiff');
berTheory_2 = berfading(EbNoVec,'psk',M,1);

figure
semilogy(EbNoVec,berVec(:,1),'*')
hold on
semilogy(EbNoVec,berTheory_2)
legend('Simulation','Theory','Location','Best')
xlabel('Eb/No (dB)')
ylabel('Bit Error Rate')
grid on
hold off

Answer 1

[Ahmed Aladi]

clear all;
M = 4;                 % Modulation alphabet
k = log2(M);           % Bits/symbol
numSC = 128;           % Number of OFDM subcarriers
cpLen = 32;            % OFDM cyclic prefix length
maxBitErrors = 1e3;    % Maximum number of bit errors
maxNumBits = 1e7;      % Maximum number of bits transmitted
qpskMod = comm.QPSKModulator('BitInput',true);
qpskDemod = comm.QPSKDemodulator('BitOutput',true);
ofdmMod = comm.OFDMModulator('FFTLength',numSC,'CyclicPrefixLength',cpLen);
ofdmDemod = comm.OFDMDemodulator('FFTLength',numSC,'CyclicPrefixLength',cpLen);
awgnChannel = comm.AWGNChannel('NoiseMethod','Variance', 'VarianceSource','Input port');
errorRate = comm.ErrorRate('ResetInputPort',true);
ofdmDims = info(ofdmMod);
numDC = ofdmDims.DataInputSize(1);
frameSize = [k*numDC 1];
EbNoVec = (0:20)';
snrVec = EbNoVec + 10*log10(k) + 10*log10(numDC/numSC);
berVec = zeros(length(EbNoVec),3);
errorStats = zeros(1,3);
for m = 1:length(EbNoVec)
    snr = snrVec(m);
      while errorStats(2) <= maxBitErrors && errorStats(3) <= maxNumBits
          dataIn = randi([0,1],frameSize);              % Generate binary data
          qpskTx = qpskMod(dataIn);                     % Apply QPSK modulation
          txSig = ofdmMod(qpskTx);                      % Apply OFDM modulation
          %fadedSig = rayleighChannel(txSig);
          %% Rayleigh Fading Channel
          taps=3;  %number of channel taps
          h = [randn(1,taps)+j*randn(1,taps) zeros(1,size(txSig,1)-taps)].' ; % generate channel taps and pad with zeros
          h_t = toeplitz(h); %  channel matrix  
          fadedSig=h_t*txSig; % Apply Rayleigh fading
          powerDB = 10*log10(var(fadedSig));               % Calculate Tx signal power
          noiseVar = 10.^(0.1*(powerDB-snr));           % Calculate the noise variance
          rxSig = awgnChannel(fadedSig,noiseVar);          % Pass the signal through a noisy channel
          rxSig_eq = inv(h_t'*h_t)*(h_t')*rxSig; % zero forcing equalization
         % Apply OFDM demodulation
          qpskRx = ofdmDemod(rxSig_eq);  
          dataOut = qpskDemod(qpskRx);                  % Apply QPSK demodulation
          errorStats = errorRate(dataIn,dataOut,0);     % Collect error statistics
      end
      berVec(m,:) = errorStats;                         % Save BER data
      errorStats = errorRate(dataIn,dataOut,1);         % Reset the error rate calculator
  end
berTheory = berawgn(EbNoVec,'psk',M,'nondiff');
berTheory_2 = berfading(EbNoVec,'psk',M,1);
figure
semilogy(EbNoVec,berVec(:,1),'*')
hold on
semilogy(EbNoVec,berTheory_2)
legend('Simulation','Theory','Location','Best')
xlabel('Eb/No (dB)')
ylabel('Bit Error Rate')
grid on
hold off

[Yasantha Samarawickrama]

Hi Ahmed! Many thanks for the solution. It works fine with the tap based Rayleigh channel. However, when i use comm.RayleighChannel block it gives i Bit error rate around 0.5 like random decoding. can't we use comm.RayleighChannel block in MATLAB so then we can customize channel properties better? Thanks in advance!

[Ahmed Aladi]

Enable PathGainsOutputPort then get the path gains and use them for the equalization.

clear all;
M = 4;                 % Modulation alphabet
k = log2(M);           % Bits/symbol
numSC = 128;           % Number of OFDM subcarriers
cpLen = 32;            % OFDM cyclic prefix length
maxBitErrors = 5e3;    % Maximum number of bit errors
maxNumBits = 1e7;      % Maximum number of bits transmitted

qpskMod = comm.QPSKModulator('BitInput',true);
qpskDemod = comm.QPSKDemodulator('BitOutput',true);

ofdmMod = comm.OFDMModulator('FFTLength',numSC,'CyclicPrefixLength',cpLen);
ofdmDemod = comm.OFDMDemodulator('FFTLength',numSC,'CyclicPrefixLength',cpLen);

rayleighChannel = comm.RayleighChannel;
rayleighChannel.PathGainsOutputPort =true;

%rayleighChannel.PathDelays= [0 1e-3 3e-3 ];
%rayleighChannel.AveragePathGains=[0 -1 -2 ];

awgnChannel = comm.AWGNChannel('NoiseMethod','Variance', 'VarianceSource','Input port');

errorRate = comm.ErrorRate('ResetInputPort',true);

ofdmDims = info(ofdmMod);

numDC = ofdmDims.DataInputSize(1);

frameSize = [k*numDC 1];

EbNoVec = (0:20)';
snrVec = EbNoVec + 10*log10(k) + 10*log10(numDC/numSC);

berVec = zeros(length(EbNoVec),3);
errorStats = zeros(1,3);

for m = 1:length(EbNoVec)
    snr = snrVec(m);

      while errorStats(2) <= maxBitErrors && errorStats(3) <= maxNumBits
          dataIn = randi([0,1],frameSize);              % Generate binary data
          qpskTx = qpskMod(dataIn);                     % Apply QPSK modulation
          txSig = ofdmMod(qpskTx);                      % Apply OFDM modulation
          [fadedSig pathGains] = rayleighChannel(txSig); % Apply fading and get the path gains
          powerDB = 10*log10(var(txSig));               % Calculate Tx signal power
          noiseVar = 10.^(0.1*(powerDB-snr));           % Calculate the noise variance
          rxSig = awgnChannel(fadedSig,noiseVar);          % Pass the signal through a noisy channel
          rxSig_eq=rxSig./sum(pathGains, 2);             % Equalize
          qpskRx = ofdmDemod(rxSig_eq);                    % Apply OFDM demodulation
          dataOut = qpskDemod(qpskRx);                  % Apply QPSK demodulation
          errorStats = errorRate(dataIn,dataOut,0);     % Collect error statistics
      end

      berVec(m,:) = errorStats;                         % Save BER data
      errorStats = errorRate(dataIn,dataOut,1);         % Reset the error rate calculator
  end

berTheory = berawgn(EbNoVec,'psk',M,'nondiff');
berTheory_2 = berfading(EbNoVec,'psk',M,1);

figure
semilogy(EbNoVec,berVec(:,1),'*')
hold on
semilogy(EbNoVec,berTheory_2)
legend('Simulation','Theory','Location','Best')
xlabel('Eb/No (dB)')
ylabel('Bit Error Rate')
grid on
hold off

[Yasantha Samarawickrama]

Hi Ahmed! Many thanks for your reply. It works fine for Frequency flat fading channels using comm.Rayleighchannel block. However, when the channel becomes Frequency selective, again the BER is high like random decoding. So do you use a different type of Equalization when channel is Frequency selective? Please include a code if you have got. Thanks in advance!

[Maggie Shammaa]

Have you found a solution to this problem please??

[alza jabar]

hi, please the place of mmse and zf equalizers before or after the ofdm demodulation?

[Pablo Sanz]

Hi, any solution to the frequency selective problem?

Answer 2

[N/A]

I think the problem is when you generate the rayleigh channel, try to add some inputs such as the sampling time/frequency, the maximum doppler shift and you can add also the delays and there pdp,

Here is an exampe of how I simulate an EVA channel :

DS = 1e-9;

Tau = [0 30 150 310 370 710 1090 1730 2510]*DS;

Pdp = [0 -1.5 -1.4 -3.6 -0.6 -9.1 -7 -12 -16.9];

fd = 20/3.6*fc/c; % Maximum Doppler Shift

chan = rayleighchan(Ts, fd, Tau, Pdp);

with Ts being the sampling time,

Answer 3

[N/A]

Did you find the problem with your program ?

Answer 4

[mustafa sami]

Did you find the problem with your program ?

Answer 5

[Yasantha Samarawickrama]

I am encountering the same problem? Can you please share if you found the solution for this issue?

Answer 6

[Yue Shang]

Some thoughts on the solution using comm.RayleighChannel

1. The frequency-flat fading channel just applies a gain to each sample to the channel. The gain is the second output of the comm.RayleighChannel object's step method, i.e., pathGains. So the time-domain equalization can be done as rxSig_eq=rxSig./pathGains for ZF or rxSig_eq = rxSig.*conj(pathGains)./(abs(pathGains).^2+noiseVar) for MMSE.
2. Set the MaximumDopplerShift property of comm.RayleighChannel to a larger number such as 0.1 for the channel to vary more quickly so you have the channel's "average" behavior in a shorter period of time. As a result, the maxBitErrors can be reduced for a shorter simulation.
3. For this frequency-flat fading channel, the time-domain equation is simpler. For a frequency-selective channel with multipath (when PathDelays and AveragePathGains properties of comm.RayleighChannel are vectors), using a frequency-domain equalization after OFDM demodulation is usually simpler. You can refer to the Indoor MIMO-OFDM Communication Link using Ray Tracing example and in particular the helperIndoorRayTracingRxProcessing function in it for how this can be done.