In this project, we create a link-level Monte-Carlo simulator. Then, we analyze the performance of convolutional codes over different wireless channels: additive white Gaussian noise(AWGN) and Rayleigh fading.
We use following blocks to implement our simulator: source generator, encoder/decoder, inter-leaver/deinterleaver, modulator/demodulator, and wireless channel.
Development tool: Matlab 2012a with Communication toolbox
function ObjectActive(Btn_handles, Line_handles, status)
if strcmp(status , 'on')
set(Btn_handles, 'Enable','on');
set(Line_handles, 'ForegroundColor','r');
elseif strcmp(status , 'off')
set(Btn_handles, 'Enable','off');
set(Line_handles, 'ForegroundColor','k');
end
%--------------------------------------------------------------------------
% Set the delay time of the objects in the Model Panel
function ActionDelay(status, sec)
if status == 1
pause(sec);
elseif status == 0
pause(0);
end
%--------------------------------------------------------------------------
%Start Simulation
function Btn_Simulation_Callback(hObject, eventdata, handles)
% -----------------Start Initialization-----------------
ObjectActive(handles.Btn_Model_Bernoulli, handles.Line1_1, 'off');
ObjectActive(handles.Btn_Model_Encoder, handles.Line1, 'off');
ObjectActive(handles.Btn_Model_Interleaver, handles.Line2, 'off');
ObjectActive(handles.Btn_Model_Modulator, handles.Line3, 'off');
ObjectActive(handles.Btn_Model_AWGN, handles.Line4, 'off');
ObjectActive(handles.Btn_Model_Rayleigh, handles.Line4, 'off');
ObjectActive(handles.Btn_Model_Demodulator, handles.Line5, 'off')
ObjectActive(handles.Btn_Model_Deinterleaver, handles.Line6, 'off');
ObjectActive(handles.Btn_Model_Decoder, handles.Line7, 'off');
ObjectActive(handles.Btn_Model_ErrorRateCalc, handles.Line8, 'off');
dDelayTime = 0;
DelayChecked = get(handles.CheckBox_Delay, 'value');
ActionDelay(DelayChecked, dDelayTime);
global gFrameSize_F;
FrameNum_K = str2double(get(handles.Edit_FrameNum,'string')); % Number of Frames (K)
gFrameSize_F = str2double(get(handles.Edit_FrameSize,'string')); % Frame Size (F): bits
ShiftBit_k = str2double(get(handles.Edit_NumBitShit,'string')); % Number of Bits Shifted (k): default k=1
EncoderOutputBit_N = str2double(get(handles.Edit_NumEncoderOutputBit,'string')); % Number of Encoder Output Bits (N)
CodeRate_Rc = ShiftBit_k / EncoderOutputBit_N; % Code Rate (Rc=k/N)
global gTimesOfOutputBits;
gTimesOfOutputBits = EncoderOutputBit_N / ShiftBit_k;
set(handles.Edit_CodeRate,'string', CodeRate_Rc);
GeneratorPoly = str2num(get(handles.Edit_GeneratorPoly,'string'));
MemoryNum_M = str2double(get(handles.Edit_NumMemory,'string')); % Number of Memory Elements (M)
MatrixSize = size(GeneratorPoly);
% Verify Trellis Structure
if MatrixSize(1,2) ~= gTimesOfOutputBits
set(handles.Edit_NumEncoderOutputBit, 'ForegroundColor','r');
set(handles.Edit_GeneratorPoly, 'ForegroundColor','r');
error_msg = 'Check the values of "Number of Encoder Output Bits" and "Generator Polynomials"';
msgbox(error_msg,'Error','error');
return;
else
set(handles.Edit_NumEncoderOutputBit, 'ForegroundColor','k');
set(handles.Edit_GeneratorPoly, 'ForegroundColor','k');
end
SNR_dB = str2double(get(handles.Edit_SNR_dB,'string')); % AWGN Parameter: Signal to Noise Ratio (SNR)
MobileSpeed_Low = str2double(get(handles.Edit_LowMobility,'string')); % Rayleigh Parameter: Mobile Speeds - Low (v)
MobileSpeed_High = str2double(get(handles.Edit_HighMobility,'string')); % Rayleigh Parameter: Mobile Speeds - High (v)
CarrierFrequency = str2double(get(handles.Edit_CarrierFrequency,'string')); % Rayleigh Parameter: Carrier Frequency (fc)
BER_Target = str2double(get(handles.Edit_BER_Target,'string')); % Target BER
SelectedVitebiDecisionType = (get(handles.ListBox_VitebiDecisionType,'Value')); % Viterbi Convolutional Encoder Parameter: Decision Type
if SelectedVitebiDecisionType == 1
VitebiDecisionType = 'hard';
fprintf('Vitebi Decision Type: Hard Decision Decoding (HDD)\n');
elseif SelectedVitebiDecisionType == 2
VitebiDecisionType = 'soft';
fprintf('Vitebi Decision Type: Soft Decision Decoding (SDD)\n');
end
SelectedWirelessChannelType = (get(handles.WirelessChannel_AWGN,'Value')); % Selecting Wireless Channel Type
if SelectedWirelessChannelType == 1
WirelessChannelType = 1;
fprintf('Wireless Channel Type: AWGN\n');
elseif SelectedWirelessChannelType == 0
WirelessChannelType = 2;
fprintf('Wireless Channel Type: Rayleigh Fading\n');
SelectedRayleighSpectrumType = (get(handles.ListBox_SpectrumType,'Value')); % Rayleigh Parameter: Doppler Spectrum Type
if SelectedRayleighSpectrumType == 1
RayleighSpectrumType = 'flat';
fprintf('Rayleigh Spectrum Type: Flat\n');
elseif SelectedRayleighSpectrumType == 2
RayleighSpectrumType = 'jakes';
fprintf('Rayleigh Spectrum Type: Jakes\n');
end
end
% -----------------End Initialization-----------------
% Generate data and apply fading channel.
M = 2; % DBPSK modulation order
hModulator = comm.BPSKModulator; % Create a DPSK modulator
hDemod = comm.BPSKDemodulator; % Create a DPSK demodulator
Bernoulli = randi([0 1],gFrameSize_F,1);
dpskSig = step(hModulator, Bernoulli); % DPSK modulate the signal
chan = rayleighchan(1/10000,100);
fadedSig = filter(chan,dpskSig); % Apply the channel effects
% Compute error rate for different values of SNR.
SNR = 0:2:20; % Range of SNR values, in dB.
numSNR = length(SNR);
berVec = zeros(3, numSNR);
% Create an AWGNChannel and ErrorRate calculator System object
hChan = comm.AWGNChannel('NoiseMethod', 'Signal to noise ratio (SNR)');
hErrorCalc = comm.ErrorRate;
for n = 1:numSNR
hChan.SNR = SNR(n);
rxSig = step(hChan,fadedSig); % Add Gaussian noise
DeModData = step(hDemod,rxSig);
reset(hErrorCalc)
% Compute error rate.
berVec(:,n) = step(hErrorCalc,Bernoulli,DeModData);
end
BER = berVec(1,:);
% Compute theoretical performance results, for comparison.
BERtheory = berfading(SNR,'dpsk',2,1);
figure
% Plot BER results.
semilogy(SNR,BERtheory,'b-',SNR,BER,'r*');
legend('Theoretical BER','Empirical BER');
xlabel('SNR (dB)'); ylabel('BER');
title('Binary DPSK over Rayleigh Fading Channel');
% Step 1. Generating Bernoulli Source
ObjectActive(handles.Btn_Model_Bernoulli, handles.Line1_1, 'on');
global gBernoulli
gBernoulli = Bernoulli_Binary_Generator(gFrameSize_F);
ActionDelay(DelayChecked, dDelayTime);
% Step 2. Convolutional Encoder
ObjectActive(handles.Btn_Model_Encoder, handles.Line1, 'on');
trel = poly2trellis(MemoryNum_M , GeneratorPoly);
PolySize = size(GeneratorPoly);
global gConvolutionalCode;
gConvolutionalCode = ConvolutionalEncoder(gFrameSize_F,gBernoulli,EncoderOutputBit_N,ShiftBit_k,MemoryNum_M,trel);
ActionDelay(DelayChecked, dDelayTime);
% 1 -> -1 and 0 -> 1
% Step 3. Block Interleaver
ObjectActive(handles.Btn_Model_Interleaver, handles.Line2, 'on');
permVec = randperm(gFrameSize_F*gTimesOfOutputBits)';
global gInverleaverData;
gInverleaverData = Interleaver('BlockInterleaver', gConvolutionalCode', permVec);
ActionDelay(DelayChecked, dDelayTime);
% Step 4. BPSK Modulator
ObjectActive(handles.Btn_Model_Modulator, handles.Line3, 'on');
global gModData;
gModData = BPSK_Modulator(gInverleaverData);
ActionDelay(DelayChecked, dDelayTime);
% Step 5. Wireless Channel: AWGN and Rayleigh
global gNoiseCode;
gNoiseCode = WirelessChannel_Rayleigh(RayleighSpectrumType,CarrierFrequency,MobileSpeed_Low,gModData); % not finished yet
% Step 6. BPSK Demodulator
ObjectActive(handles.Btn_Model_Demodulator, handles.Line5, 'on');
global DeModData;
DeModData = BPSK_Demodulator(gNoiseCode);
ActionDelay(DelayChecked, dDelayTime);
% Step 7. Block Deinterleaver
ObjectActive(handles.Btn_Model_Deinterleaver, handles.Line6, 'on');
global gDegInverleaverData;
gDegInverleaverData = Deinterleaver('BlockInterleaver',DeModData,permVec);
ActionDelay(DelayChecked, dDelayTime);
% Step 8. Viterbi Decoder
ObjectActive(handles.Btn_Model_Decoder, handles.Line7, 'on');
global gDecoded;
gDecoded = ViterbiDecoder(gDegInverleaverData, VitebiDecisionType, 'cont', 2,trel); % Issue: soft doesn't work
ActionDelay(DelayChecked, dDelayTime);
% Step 9. Error Rate Calc
ObjectActive(handles.Btn_Model_ErrorRateCalc, handles.Line8, 'on');
rate = ErrorRateCalc(gBernoulli,gDecoded)
ActionDelay(DelayChecked, dDelayTime);
global gresult_msg;
fprintf('Error Rate: %d\n', rate);
%--------------------------------------------------------------------------
% Step 1. Generating Bernoulli Source
% Sequence: 0 and 1
function Bernoulli = Bernoulli_Binary_Generator(frame_size)
%global Bernoulli;
Bernoulli = randi([0 1],1,frame_size);
%--------------------------------------------------------------------------
% Step 2. Convolutional Encoder
function ConvolutionalCode = ConvolutionalEncoder(leng, source, Output_N, shift_k, Memory_Element_M, Trellis)
hConvEnc = comm.ConvolutionalEncoder(Trellis);
ConvolutionalCode = step(hConvEnc,source')'; % Encode the message.
%--------------------------------------------------------------------------
% tblen is a positive integer scalar that specifies the traceback depth.
% If the code rate is 1/2, a typical value for tblen is about 5 times the constraint length of the code.
% http://www.mathworks.com/help/comm/ref/vitdec.html
function Decoded = ViterbiDecoder(NoiseSource, DecissionType, Opmode, TracebackDepth, Trellis)
if strcmp(DecissionType , 'hard')
hVitDec = comm.ViterbiDecoder(Trellis, 'InputFormat', 'hard','TracebackDepth', TracebackDepth, 'TerminationMethod', 'Truncated');
Decoded = step(hVitDec, NoiseSource); % Decode.
elseif strcmp(DecissionType , 'soft')
% To prepare for soft-decision decoding, map to decision values.
[x,qcode] = quantiz(1-2*NoiseSource,[-.75 -.5 -.25 0 .25 .5 .75],[7 6 5 4 3 2 1 0]);
Decoded = vitdec(qcode,Trellis,TracebackDepth,Opmode,DecissionType, 3); % Decode.
end
%--------------------------------------------------------------------------
% Step 5. Wireless Channel: AWGN
function NoiseCode = WirelessChannel_AWGN(leng,source, SNR)
hAWGN = comm.AWGNChannel('NoiseMethod', 'Signal to noise ratio (SNR)','SNR',SNR);
%hAWGN = comm.AWGNChannel('NoiseMethod', 'Signal to noise ratio (SNR)','SNR',SNR);
NoiseCode = step(hAWGN, source);
%--------------------------------------------------------------------------
% Step 5. Wireless Channel: Rayleigh
function RayleighData = WirelessChannel_Rayleigh(DopplerType, Frequency, Speed, NoiseSignal)
% ts is the sample time of the input signal, in seconds. fd is the maximum Doppler shift, in hertz.
fm = DopplerShift((Speed*1000)/3600,Frequency*10^6);
if strcmp(DopplerType,'flat') == 1
channel = rayleighchan(1/Frequency,fm); % Create object.
channel.DopplerSpectrum = doppler.flat;
elseif strcmp(DopplerType,'jakes') == 1
channel = rayleighchan(1/Frequency,fm); % Create object.
end
data=filter(channel,NoiseSignal);
% http://www.mathworks.com/matlabcentral/answers/13373-ber-calculation-bpsk-simulation-with-rayleigh-channel
hAWGN = comm.AWGNChannel('NoiseMethod', 'Signal to noise ratio (SNR)','SNR',10);
%hAWGN = comm.AWGNChannel('NoiseMethod','Signal to noise ratio (SNR)','SNR',SNR);
% NoiseCode = step(hAWGN, source);
RayleighData = step(hAWGN,data); % Add Gaussian noise
%RayleighData = step(hDemod, rxSig); % Demodulate
%--------------------------------------------------------------------------
function intData = Interleaver(InterleaverType, SourceData, PermutationVector)
if strcmp(InterleaverType,'BlockInterleaver')
hInt = comm.BlockInterleaver(PermutationVector);
intData = step(hInt, SourceData);
end
%--------------------------------------------------------------------------
function deIntData = Deinterleaver(InterleaverType, SourceData, PermutationVector)
if strcmp(InterleaverType,'BlockInterleaver')
hDeInt = comm.BlockDeinterleaver(PermutationVector);
deIntData = step(hDeInt, SourceData);
end
%--------------------------------------------------------------------------
function gModData = BPSK_Modulator(Data)
% Create a BPSK modulator System object
hModulator = comm.BPSKModulator;
%constellation(hModulator);
% Modulate the data
gModData = step(hModulator, Data);
%--------------------------------------------------------------------------
function DeModData = BPSK_Demodulator(noisySignal)
hDemod = comm.BPSKDemodulator;
DeModData = step(hDemod,noisySignal);
%--------------------------------------------------------------------------
function fm = DopplerShift(v,fc)
fm = (v/(3*10^8))*fc; %Hz
%--------------------------------------------------------------------------
function rate = ErrorRateCalc(SourceCode, DecodedCode)
hErrorCalc = comm.ErrorRate;
reset(hErrorCalc)
[number,rate] = biterr(DecodedCode',SourceCode);
%--------------------------------------------------------------------------
function DrawSquareWave(title_name,x_axis_label, y_axis_label, bit_array, min_x, max_x, min_y, max_y)
figure
t= 0:max_x;
bit_array(max_x+1) = bit_array(max_x); % draw the last bit
stairs(t, bit_array,'LineWidth',2);
axis([min_x max_x min_y max_y])
%grid on
xlabel(x_axis_label);
ylabel(y_axis_label);
title(title_name);
%--------------------------------------------------------------------------