clear;

clc;

s = rng(211);            % 设置RNG状态的可重复性

numFFT = 1024;           % FFT点的数量

numRBs = 50;             % 资源块的数量

rbSize = 12;             % 每个资源块的副载波数量

cpLen = 72;              % 样本中的循环前缀长度

bitsPerSubCarrier = 6;   % 2: QPSK, 4: 16QAM, 6: 64QAM, 8: 256QAM

snrdB =30;              % 以dB为单位的SNR

toneOffset = 2.5;        % 音调偏移或超额带宽(在副载波中)

L = 513;                 % 滤波器长度(=滤波器次序+ 1)，奇数

numDataCarriers = numRBs*rbSize;    % 子带中数据子载波的数量

halfFilt = floor(L/2);

n = -halfFilt:halfFilt;

% Sinc函数原型滤波器

pb = sinc((numDataCarriers+2*toneOffset).*n./numFFT);

% Sinc截断窗口

w = (0.5*(1+cos(2*pi.*n/(L-1)))).^0.6;

% 归一化的低通滤波器系数

fnum = (pb.*w)/sum(pb.*w);

% 过滤脉冲响应

h = fvtool(fnum, 'Analysis', 'impulse', ...

'NormalizedFrequency', 'off', 'Fs', 15.36e6);

h.CurrentAxes.XLabel.String = 'Time (\mus)';

h.FigureToolbar = 'off';

% 使用dsp过滤器对象进行过滤

filtTx = dsp.FIRFilter('Structure', 'Direct form symmetric', ...

'Numerator', fnum);

filtRx = clone(filtTx); % Rx的匹配过滤器

% QAM符号映射器

qamMapper = comm.RectangularQAMModulator( ...

'ModulationOrder', 2^bitsPerSubCarrier, 'BitInput', true, ...

'NormalizationMethod', 'Average power');

% 为谱图设置一个数字

hFig = figure('Position', figposition([46 50 30 30]), 'MenuBar', 'none');

axis([-0.5 0.5 -200 -20]);

hold on;

grid on

xlabel('Normalized frequency');

ylabel('PSD (dBW/Hz)')

title(['F-OFDM, ' num2str(numRBs) ' Resource blocks, '  ...

num2str(rbSize) ' Subcarriers each'])

% 生成数据符号

bitsIn = randi([0 1], bitsPerSubCarrier*numDataCarriers, 1);

symbolsIn = qamMapper(bitsIn);

% 将数据打包成OFDM符号

offset = (numFFT-numDataCarriers)/2; % 为通带中心

symbolsInOFDM = [zeros(offset,1); symbolsIn; ...

zeros(numFFT-offset-numDataCarriers,1)];

ifftOut = ifft(ifftshift(symbolsInOFDM));

%预先加上循环前缀

txSigOFDM = [ifftOut(end-cpLen+1:end); ifftOut];

% 使用零填充过滤器来刷新尾部。 获取发射信号

txSigFOFDM = filtTx([txSigOFDM; zeros(L-1,1)]);

%绘制功率谱密度(PSD)

[psd,f] = periodogram(txSigFOFDM, rectwin(length(txSigFOFDM)), ...

numFFT*2, 1, 'centered');

plot(f,10*log10(psd));

% 计算峰均功率比(PAPR)

PAPR = comm.CCDF('PAPROutputPort', true, 'PowerUnits', 'dBW');

[~,~,paprFOFDM] = PAPR(txSigFOFDM);

disp(['Peak-to-Average-Power-Ratio for F-OFDM = ' num2str(paprFOFDM) ' dB']);

% 绘制OFDM信号的功率谱密度(PSD)

[psd,f] = periodogram(txSigOFDM, rectwin(length(txSigOFDM)), numFFT*2, ...

1, 'centered');

hFig1 = figure('Position', figposition([46 15 30 30]));

plot(f,10*log10(psd));

grid on

axis([-0.5 0.5 -200 -20]);

xlabel('Normalized frequency');

ylabel('PSD (dBW/Hz)')

title(['OFDM, ' num2str(numRBs*rbSize) ' Subcarriers'])

%计算峰均功率比(PAPR)

PAPR2 = comm.CCDF('PAPROutputPort', true, 'PowerUnits', 'dBW');

[~,~,paprOFDM] = PAPR2(txSigOFDM);

disp(['Peak-to-Average-Power-Ratio for OFDM = ' num2str(paprOFDM) ' dB']);

%添加AWGN和Rayleigh

chan = rayleighchan(0.00001,100);

rxSig = filter(chan,txSigFOFDM);

rxtSig = awgn(rxSig, snrdB, 'measured');

% 接收匹配的过滤器,

rxSigFilt = filtRx(rxtSig);

% 考虑过滤器延迟

rxSigFiltSync = rxSigFilt(L:end);

% 删除循环前缀

rxSymbol = rxSigFiltSync(cpLen+1:end);

%执行FFT

RxSymbols = fftshift(fft(rxSymbol));

%选择数据副载波

dataRxSymbols = RxSymbols(offset+(1:numDataCarriers));

% 绘制接收符号星座图

switch bitsPerSubCarrier

case 2  % QPSK

refConst = qammod((0:3).', 4, 'UnitAveragePower', true);

case 4  % 16QAM

refConst = qammod((0:15).', 16,'UnitAveragePower', true);

case 6  % 64QAM

refConst = qammod((0:63).', 64,'UnitAveragePower', true);

case 8  % 256QAM

refConst = qammod((0:255).', 256,'UnitAveragePower', true);

end

constDiagRx = comm.ConstellationDiagram( ...

'ShowReferenceConstellation', true, ...

'ReferenceConstellation', refConst, ...

'Position', figposition([20 15 30 40]), ...

'MeasurementInterval', length(dataRxSymbols), ...

'Title', 'F-OFDM Demodulated Symbols', ...

'Name', 'F-OFDM Reception', ...

'XLimits', [-1.5 1.5], 'YLimits', [-1.5 1.5]);

constDiagRx(dataRxSymbols);

% 由于没有通道建模，因此通道均衡不是必要的

%解映射和BER计算

qamDemod = comm.RectangularQAMDemodulator('ModulationOrder', ...

2^bitsPerSubCarrier, 'BitOutput', true, ...

'NormalizationMethod', 'Average power');

BER = comm.ErrorRate;

% 执行艰难的决策并测量错误

rxBits = qamDemod(dataRxSymbols);

ber = BER(bitsIn, rxBits);

disp(['F-OFDM Reception, BER = ' num2str(ber(1)) ' at SNR = ' ...

num2str(snrdB) ' dB']);

%恢复RNG状态

rng(s);