clear all; close all; clc;

% Коэффициент ослабления лучей
r = [0.057662, 0.176809, 0.407163, 0.303585, 0.258782, ...
0.061831, 0.150340, 0.051534, 0.185074, 0.400967, ...
0.295723, 0.350825, 0.262909, 0.225894, 0.170996, ...
0.149723, 0.240140, 0.116587, 0.221155, 0.259730;];

% Задержка лучей
tau = [1.003019, 5.422091, 0.518650, 2.751772, 0.602895, ...
1.016585, 0.143556, 0.153832, 3.324866, 1.935570, ...
0.429948, 3.228872, 0.848831, 0.073883, 0.203952, ...
0.194207, 0.924450, 1.381320, 0.640512, 1.368671;] * 1e-6;


SNR = 20; % Отношение сигнал шум в дБ
NumBlock = 100; % Колличество символов прошедших через канал
SamplesRate = (2.5e6/3)*8/5; % Частота передискретизированного сигнала
NumCarriers = 533; % Колличество несущих
FFTLength = 3000; % Колличество отсчетов полезной части сигнала
GuardLength = 375; % Колличество отсчетов защитной части сигнала


% Расчет Доплеровского смещения
FreqCarrier = 108e6; % Несущая сигнала
Speed = 200e3/3600; % Скорость движения приемника относительно передатчика
c = 300e6; % Скорость света
DopplerShift = FreqCarrier*Speed/c; % Максимальный доплеровский сдвиг

SampleTime = 1/SamplesRate; % Период дискретизации сигнала
DelayVector = tau;
GainVector = 10*log10®; % Коэффициент ослабления лучей в дБ
scale = 1./sqrt(sum(r.^2));

ChanObj = rayleighchan(SampleTime,DopplerShift,DelayVector,GainVector);
ChanObj.ResetBeforeFiltering = 1;
ChanObj.NormalizePathGains = 1;
ChanObj.StorePathGains = 0;
ChanObj.DopplerSpectrum = doppler.gaussian(DopplerShift/2);
ChanObj.StoreHistory = 1;

M=4;
k=log2(M);
H_QAM4=modem.qammod('M',M,'PhaseOffset', 0, 'SymbolOrder',...
'binary', 'InputType', 'bit');
hDeMod = modem.qamdemod(H_QAM4);

data = [];
for n = 1:NumBlock
Input_Data = randint(k*NumCarriers,1,2);
QAM_tx = modulate(H_QAM4,Input_Data);

Symbol = ifft([QAM_tx(ceil(NumCarriers/2)+1:NumCarriers);complex(zeros(FFTLength-NumCarriers,1),zeros(FFTLength-NumCarriers,1));...
QAM_tx(1:ceil(NumCarriers/2))]);

SymbolCP = [Symbol(FFTLength-GuardLength+1:FFTLength);Symbol];

OutputBlock = filter(ChanObj, SymbolCP);
OutputBlock = awgn(OutputBlock, SNR+10*log10(k)-10*log10(1),'measured');

Symbol = fft(OutputBlock(GuardLength+1:FFTLength+GuardLength));
QAM_rx(1,1:ceil(NumCarriers/2)) = Symbol(FFTLength - ceil(NumCarriers/2)+1:FFTLength);
QAM_rx(1,ceil(NumCarriers/2):NumCarriers) = Symbol(1:ceil(NumCarriers/2));

end

plot(linspace(0,SamplesRate,3375), 20*log10(abs(fft(SymbolCP)))); grid on;
%plot(ChanObj);