subcarrierSpacing = 15e3;
nFFT = ceil(totalBandwidth / subcarrierSpacing);
noisePower_W = 10^(noisePower_dBW / 10);
impulseNoisePower_dBW = -113;
impulseNoisePower_W = 10^(impulseNoisePower_dBW / 10);
berEstAWGN = zeros(size(EbNoVec));
berEstImp = zeros(size(EbNoVec));
snr_dB = convertSNR(EbNoVec, "ebno", "snr", BitsPerSymbol=bps);
for n = 1:length(EbNoVec)
while numErrsAWGN < 200 && numBits < 1e6
txSymbols = randi([0 M-1], nFFT, 1);
txBits = int2bit(txSymbols, bps);
txGrid = qammod(txSymbols, M, 'UnitAveragePower', true);
txOut = ifft(txGrid, nFFT);
signalPower_dBW = snr_dB(n) + noisePower_dBW;
signalPower_W = 10^(signalPower_dBW / 10);
scalingFactor = sqrt(signalPower_W / mean(abs(txOut).^2));
txOutScaled = txOut * scalingFactor;
noiseSignal_notscaled = randn(size(txOutScaled)) + 1j * randn(size(txOutScaled));
noiseAWGN = sqrt(noisePower_W / 2) * noiseSignal_notscaled;
impulseNoise = zeros(size(txOutScaled));
burstPositions = rand(size(txOutScaled)) < 1;
impulseNoise(burstPositions) = sqrt(impulseNoisePower_W / 2) * ...
(randn(sum(burstPositions), 1) + 1j * randn(sum(burstPositions), 1));
rxInAWGN = txOutScaled + noiseAWGN;
rxInImp = txOutScaled + noiseAWGN + impulseNoise;
rxGridAWGN = fft(rxInAWGN, nFFT);
rxGridImp = fft(rxInImp, nFFT);
rxSymbolsAWGN = qamdemod(rxGridAWGN, M, 'UnitAveragePower', true);
rxSymbolsImp = qamdemod(rxGridImp, M, 'UnitAveragePower', true);
rxBitsAWGN = int2bit(rxSymbolsAWGN, bps);
rxBitsImp = int2bit(rxSymbolsImp, bps);
[bitErrorsAWGN, ~] = biterr(txBits, rxBitsAWGN);
[bitErrorsImp, ~] = biterr(txBits, rxBitsImp);
numErrsAWGN = numErrsAWGN + bitErrorsAWGN;
numErrsImp = numErrsImp + bitErrorsImp;
numBits = numBits + nFFT * bps;
berEstAWGN(n) = numErrsAWGN / numBits;
berEstImp(n) = numErrsImp / numBits;
semilogy(EbNoVec, berEstAWGN, 'b-o', 'LineWidth', 2);
semilogy(EbNoVec, berEstImp, 'r--o', 'LineWidth', 2);
title('OFDM BER for BPSK with AWGN and Impulse Noise');
ylabel('Bit Error Rate (BER)');
legend('AWGN Only', 'AWGN + Impulse Noise');