%%
% This script mimic a real time simulation and describes a procedure to 
% simulate the whdlOFDMReceiver model with a custom transmitter waveform 
% and adding only the fading channel impairment and verifies the model 
% simulation channel estimates with the frequency response derived from 
% chanTaps in whdlexamples.OFDMReceiverInit.m script.
%
% You can choose a custom frame configuration in this script and a custom 
% channel type and run the script to simulate the model, initializing it 
% with custom frame configuration.
%
% The chanTaps are derived from the comm.RayleighChannel.

%% Do not modify this section

% Create a function handle for trBlksize for dataBits size per frame.
% For more information on how to calculate trBlksize see the
% whdlexamples.OFDMTx script.

trBlkSize = @(modOrder,codeRateIdx)...
    ((72-12) * 32 * log2(modOrder) *...
    ((codeRateIdx ~= 3) .* ((codeRateIdx+1)./(codeRateIdx+2)) + ...
    (codeRateIdx == 3) .* ((codeRateIdx+2)./(codeRateIdx+3)))) - ...
    (7-1)-32;

%% Custom Frame Configuration 

%====================================
% Transmitter parameter configuration
%====================================

% You can choose modOrder, codeRateIndex and numFrames for each
% structure configuration. You can add more structures for more configurations.
% For more information on creating txParam, see whdlexamples.OFDMTx.m

txParam    = struct('modOrder',4,'codeRateIndex',3,'numFrames',5);

% Generate data bits
% Data bits are generated at random. You may choose custom data bits of
% your choice

for n = 1 : length(txParam)
    txParam(n).txDataBits = randsrc(...
        trBlkSize(txParam(n).modOrder,txParam(n).codeRateIndex)*txParam(n).numFrames,1,[0 1]);
end
clear n;

%====================================
% Applying fading channel impairment
%====================================

% You can choose the type of fading channel to be applied. Other
% impairments will interfere with the fading channel response.

channelMode       = 'EVA';           % {'EPA','EVA'}

%% Do Not Modify Below Contents

% Disable all the other impairments to avoid their interference with fading
% channel estimates
EsN0dB            = inf;          % in dB
CFO               = 0;            % in Hz
CPO               = 0;            % in radians

modelname = 'whdlOFDMReceiver';
load_system(modelname);

%=======================================
% Derive fading channel impulse response
%=======================================
Fs = 1.92e6;
[pathDelays,averagePathGains] = OFDMRxFadingChannelConditions(channelMode);
rng(0);
ob = comm.RayleighChannel('SampleRate',Fs,...
    'PathDelays',pathDelays,...
    'AveragePathGains',averagePathGains,...
    'MaximumDopplerShift',0);
%         'Visualization','Impulse response');
chanTaps = ob([ones(1,1);zeros(size(info(ob).ChannelFilterCoefficients,2)-1,1)]);
chanTaps = chanTaps/sqrt(sum(abs(chanTaps).^2));
clear Fs;
        
%==============================================
% Impairment Control configuration
%==============================================
% This configuration will avoid interference due to residual frequency
% estimate on fading channel estimates
set_param([modelname '/Frequency Correction Type'],...
    'CorrType','Use externally provided Frequency Offset');


set_param([modelname '/Channel Equalization'],...
    'mode','Active');
set_param([modelname '/CPE Correction'],...
    'mode','Active');

%===========================
% simulate whdlOFDMReceiver
%===========================
sim(modelname);
simFrequencyEstimate = simFrequencyEstimate.Data(:);
simFramesSynchronized = simFramesSynchronized.Data(:);
simTimingOffset = simTimingOffset.Data(:);
simSyncSignalCorrData = simSyncSignalCorrData.Data(:);
simSyncSignalCorrThresold = simSyncSignalCorrThresold.Data(:);
simChannelEstimates = simChannelEstimates.Data(:);

%=====================================
% Display number of header CRC failed
%=====================================
FramesDecoded = simFramesSynchronized(simDetectedFramesIndices).';
fprintf('\n Number of header CRC failed = %d per %d\n',...
    uint32(numel(simFramesSynchronized)-numel(simDetectedFramesIndices)),...
    uint32(numel(simFramesSynchronized)));

%=====================================
% Display number of bit errors
%=====================================
% Bit errors are calculated only for frames whose CRC passes
matUserData = [];
for frameNumber = FramesDecoded
    matUserData = [matUserData;UserData{frameNumber}]; %#ok<AGROW>
end
nbiterrors = nnz(matUserData - simUserData);
fprintf('\n Number of bit errors = %d per %d\n',uint32(nbiterrors),uint32(numel(simUserData)));

%===========================================
% Display Fading Channel Frequency Response
%===========================================
% The transmitter waveform is passed through Tx filter, fading channel
% filter, and Rx filter. The estimates from simulink is the product of
% frequency responses of all the 3 filters. The fading channel frequency 
% response is derived by dividing simulink estimates with the frequency
% responses of Tx and Rx filters.

% Deriving the channel frequency response computed in 
% whdlexamples.OFDMReceiverInit from chanTaps 
FiltersGroupDelay = simTimingOffset(1) - 159; % timing shift
circshiftFactor = fftshift(exp(1i*2*pi*(0:127).'*(FiltersGroupDelay/128)));
ActiveSubcarrierIndices = [28+1:28+36 28+36+1+1:28+36+1+36];
matChannelFrequencyResponse = channel.FadingChannelFrequencyResponse;
% Apply the effect of timing shift and save it in a variable
FR1 = matChannelFrequencyResponse(ActiveSubcarrierIndices).*circshiftFactor(ActiveSubcarrierIndices);

% Deriving Tx and Rx filter frequency responses
Param = whdlexamples.OFDMRxParameters;
RxFilterCoeff = Param.rxFIRCoeff(:);clear Param;
RxFilterFrequencyResponse = fftshift(fft([RxFilterCoeff;zeros(128-length(RxFilterCoeff),1)]));
Param = whdlexamples.OFDMTxParameters;
TxFilterCoeff = Param.firCoeff(:);clear Param;
TxFilterFrequencyResponse = fftshift(fft([TxFilterCoeff;zeros(128-length(TxFilterCoeff),1)]));

% Deriving the channel frequency response from simulink estimates and save 
% it in a variable
FR2 = reshape(simChannelEstimates,72,[])./(TxFilterFrequencyResponse(ActiveSubcarrierIndices).*RxFilterFrequencyResponse(ActiveSubcarrierIndices));

% % Compute error in the simulation channel estimates and estimates 
% % computed from chanTaps 
% for i = 1:size(FR2,2)
%     err = FR1 - FR2(:,i);
%     errmax = [max(abs(real(err))) max(abs(imag(err)))] %#ok<NASGU>
% end
% clear i err errmax

%==========================
% Display Impulse Response
%==========================
figure;
subplot(211);
stem(abs(chanTaps),'o-');
xlabel('Channel Coefficients');
ylabel('Magnitude');
subplot(212);
stem(angle(chanTaps)/pi,'o-');
xlabel('Channel Coefficients');
ylim([-1 1]);ylabel('Normalized Phase (x \pi)');
sgtitle('Impulse Response');

%==========================
% Display Frequency Response
%==========================
figure;
subplot(211);
plot(abs(FR1),'ob-');
hold on;
plot(abs(FR2(:,1)),'.r-');
plot(abs(FR1) - abs(FR2(:,1)),'*k-');
legend('Derived from Impulse Response','Derived from Simulink Model','Error','Location','best');
xlabel('Sub carriers');
ylabel('Magnitude');
subplot(212);
plot(angle(FR1)/pi,'ob-');
hold on;
plot(angle(FR2(:,1))/pi,'.r-');
plot(angle(FR1)/pi - angle(FR2(:,1))/pi,'*k-');
legend('Derived from Impulse Response','Derived from Simulink Model','Error','Location','best');
xlabel('Sub carriers');
ylim([-1 1]);ylabel('Normalized Phase (x \pi)');
sgtitle('Frequency Response');

close_system(modelname,0);

clear ActiveSubcarrierIndices averagePathGains channelMode circshiftFactor...
    derivedChannelFrequencyResponse FiltersGroupDelay FR1 FR2...
    FramesDecoded matUserData nbiterrors ob pathDelays RxFilterCoeff...
    RxFilterFrequencyResponse TxFilterCoeff TxFilterFrequencyResponse...
    matChannelFrequencyResponse modelname
