💥💥💥💞💞💞欢迎来到本博客❤️❤️❤️💥💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
目录
💥1 概述
📚2 运行结果
🎉3 参考文献
👨💻4 Matlab代码实现
💥1 概述
2006年, Donoho等人[4提出压缩感知(Com-pressive Sensing , CS)概念框架,并用数学模型为理论提供支撑。压缩感知理论突破了奈奎斯特采样定理对信号维度的限制,避免稀疏信号在奈奎斯特采样时会产生的大量冗余信息,并缓解硬件设备和算法负担。
鲸鱼优化算法（Whale Optimization Algorithm，简称WOA）是一种启发式算法，受到鲸鱼在寻找猎物时的行为启发而提出的。它主要用于解决优化问题，包括信道估计问题。
5G信道估计是为了获取无线信道的状态信息，以便进行资源分配、功率控制等操作。在使用鲸鱼优化算法进行5G信道估计时，按照以下步骤进行：
1. 定义问题：明确信道估计的目标和约束条件。例如，可以将问题定义为最小化信道误差的均方根（Root Mean Square Error，RMSE）。
2. 初始化种群：随机生成一些鲸鱼个体作为初始种群。
3. 迭代更新：通过模拟鲸鱼的搜索行为进行迭代更新。具体来说，可以按照以下步骤进行：
- 计算适应度：根据当前鲸鱼个体的位置，计算其对应的适应度值，即信道误差。
- 更新最优个体：记录当前种群中适应度最好的个体作为当前的最优解。
- 更新位置：根据当前个体的适应度值和最优个体的位置，更新每个鲸鱼个体的位置。
- 更新搜索半径：根据当前迭代次数和最大迭代次数的比值，更新每个鲸鱼个体的搜索半径。
- 达到停止条件：判断是否达到停止条件，如果是，则停止迭代；否则，返回到第一步进行下一轮迭代。
4. 输出结果：根据迭代结束后的最优个体，得到信道估计结果。
需要注意的是，鲸鱼优化算法的效果与具体的问题、算法参数设置以及迭代次数等因素有关，因此在实际应用中，可以根据具体情况进行调整和优化。此外，还可以结合其它算法和技术，如机器学习方法和统计方法，来进一步提高5G信道估计的性能。
📚2 运行结果
部分代码：
% Channel estimation using LS, WOA and MMSE
% Number of OFDM Symbol = 1e1
% Channel model: TDLC-300
clc, clear; close all;
methods = {'LS ', 'WOA', 'MMSE'} % Channel estimation methods
snrRange = 0:5:30; % Signal to noise ratio in dB
numSymbol = 1e1; % Number of symbols
numFft = 4096; % Size of DFT
numCp = numFft/4; % Number of CP
subCarrierSpacing = 30e3; % Subcarrier Spacing
numBitPerSym = 4; % Number of bits per (modulated) symbol
numSymPerPilot = 12; % Number of (modulated) symbol per pilots
numBitBerSecond = 1e3; % Number of bits per second
signalEnergy = 100; % Energy of signal
% Propagation Channel Model Configuration
% Create a TDL channel model object and specify its propagation characteristics.
numTapEst = 400; % Number of est. channel taps
numTap = 320; % Number of true channel taps
% TDLC300-100
tapDelay = [0 65 70 190 195 200 245 325 520 1045 1510 2595]; % in ns
tapPower = [-6.9 0 -7.7 -2.5 -2.4 -9.9 -8.0 -6.6 -7.1 -13.0 -14.2 -16.0]; % in dB
% WOA Alg
maxIter = 8; % maximum number of generations
numAgent = 8; % Number of search agents
ub = [50 100 400]; % [ub1,ub2,...,ubn] where ubn is the upper bound of variable n
lb = [0 20 0]; % [lb1,lb2,...,lbn] where lbn is the lower bound of variable n
dim = 3; % Number of variables
positions = rand(numAgent, dim).*(ub-lb) + lb;
visualization = 0;
saveOrNot = 0; % = 1 for save
sampRate = numFft*subCarrierSpacing; % Sample rate
numPilot = ceil(numFft/numSymPerPilot); % Number of pilots per OFDM symbol
pilotLoc = zeros(numPilot, 1); % Pilot's Location
pathLoss = zeros(numTap, 1);
tapSample = round(tapDelay*1e-9*sampRate);
pathLoss(tapSample+1) = 10.^(tapPower/10); % Path loss of channel
M = 2^numBitPerSym; % M - QAM
A = sqrt(3/2/(M-1)*signalEnergy); % QAM normalization factor
Nofdm = numFft + numCp; % Number of OFDM
numData = numFft - numPilot; % Number of data
MSEs_snr = zeros(length(snrRange),length(methods));
ber_snr = zeros(length(snrRange),length(methods));
fileIdx = getFileId(saveOrNot);
tic
for snrIdx = 1:length(snrRange)
SNR = snrRange(snrIdx);
er = zeros(1,length(methods));
MSE = zeros(1,length(methods));
for nsym=1:numSymbol
msgint = randi([0 M-1],numFft-numPilot,1); % Symbol generation
data = qammod(msgint, M);
% Add pilot
p = randi([0, M-1], numPilot, 1); % Pilot sequence generation
pilot = qammod(p, M);
ip = 0;
X = zeros(numFft, 1);
for k=1:numFft
if rem(k,numSymPerPilot)== 1
ip = ip+1;
X(k)=pilot(floor(k/numSymPerPilot)+1); % For pilot
pilotLoc(ip) = k; % For pilot location
else
X(k) = data(k-ip); % For data
end
end
% OFDM
x = ifft(X,numFft); % IFFT
xt = [x(numFft-numCp+1:numFft); x]; % add CP
% PA
tx = A*xt;
signalPowerdB = 10*log10(cov(tx));
% Channal gain
h = (randn(numTap, 1)+1j*randn(numTap, 1))...
.*sqrt(pathLoss/2); % Channel gain
H = fft(h,numFft); % True channel frequency respond
H_power_dB = 10*log10(abs(H.*conj(H))); % True channel power in dB
y_channel = conv(tx,h); % Channel path (convolution)
% Add noise
rx = awgn(y_channel, SNR, 'measured');
% rx = y_channel + 1/(sqrt(2.0)*10^(SNR/20))*complex(randn(size(y_channel)),randn(size(y_channel)));
% sto = sto_est(rx, numFft, numCp);
% Receiver
y = rx(numCp+1:Nofdm); % Remove CP
Y = fft(y); % FFT
% Channel estimation
for methodIdx = 1:length(methods)
method = methods{methodIdx};
if method(1) == 'L'
% LS estimation with linear interpolation
H_est = LS_CE(Y,pilot,pilotLoc,numFft, 'linear');
elseif method(1) == 'W'
% WOA estimation
[H_est, positions] = WOA_CE(Y,pilot,pilotLoc,numFft, numSymPerPilot, numBitPerSym, ...
positions, numAgent, maxIter, lb, ub, dim);
elseif method(1) == 'M'
% MMSE estimation
H_est = MMSE_CE(Y,pilot,pilotLoc,numFft,numSymPerPilot,h,SNR);
end
if method(end) == 'T'
h_est = ifft(H_est); % Esti channel gain
h_est = h_est(1:numTapEst); % N-tap channel gain
H_est = fft(h_est,numFft); % DFT-based channel estimation
end
H_est_power_dB = ...
10*log10(abs(H_est.*conj(H_est))); % Esti channel power in dB
Y_eq = Y./H_est;
Data_extracted = zeros(numFft- numPilot, 1);
ip = 0;
for k=1:numFft
if mod(k,numSymPerPilot)==1
ip=ip+1;
else
Data_extracted(k-ip)=Y_eq(k);
end
end
msg_detected = qamdemod(Data_extracted, M);
bitDetected = de2bi(msg_detected, numBitPerSym);
bitTrans = de2bi(msgint, numBitPerSym);
er(methodIdx) = er(methodIdx) + sum(sum(bitDetected~=bitTrans));
MSE(methodIdx) = MSE(methodIdx) + (H-H_est/A)'*(H-H_est/A);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if visualization
figure(snrIdx)
subplot(1, length(methods), methodIdx)
hold on
scatter(real(Data_extracted), imag(Data_extracted), 'b.')
title(methods{methodIdx})
pause(0)
end
end
end
MSEs = MSE/(numFft*numSymbol);
MSEs_snr(snrIdx, :) = MSEs;
ber_snr(snrIdx, :) = er/(numSymbol*numData*numBitPerSym);
for methodIdx = 1:length(methods)
str = sprintf('SNR = %2.0f dB: BER of %s \t= %6.5f\n', SNR, methods{methodIdx}, ber_snr(snrIdx, methodIdx));
fprintf(str)
if saveOrNot
if methodIdx == 1
fprintf(fileIdx, '\n-----------------------------------------\n');
end
fprintf(fileIdx, str);
end
end
end
if saveOrNot
fclose(fileIdx);
end
figure
subplot(121)
semilogy(snrRange, ber_snr)
legend(methods{:})
xlabel('SNR')
ylabel('BER')
grid on
subplot(122)
semilogy(snrRange, MSEs_snr)
legend(methods{:})
xlabel('SNR')
ylabel('MSE')
grid on
toc
% Local Functions
function [H_interpolated] = interpolate(H,pilot_loc,Nfft,method)
% Input: H = Channel estimate using pilot sequence
% pilot_loc = Location of pilot sequence
% Nfft = FFT size
% method = 鈥檒inear鈥�/鈥檚pline鈥�
% Output: H_interpolated = interpolated channel
if pilot_loc(1)>1
slope = (H(2)-H_est(1))/(pilot_loc(2)-pilot_loc(1));
H = [H(1)-slope*(pilot_loc(1)-1); H]; pilot_loc = [1; pilot_loc];
end
if pilot_loc(end) <Nfft
slope = (H(end)-H(end-1))/(pilot_loc(end)-pilot_loc(end-1));
H = [H; H(end)+slope*(Nfft-pilot_loc(end))];
pilot_loc = [pilot_loc; Nfft];
end
if lower(method(1))=='l'
H_interpolated = interp1(pilot_loc,H,[1:Nfft]', 'linear');
else
H_interpolated = interp1(pilot_loc,H,[1:Nfft]', 'spline');
end
end
function H_LS = LS_CE(Y,Xp,pilot_loc,Nfft,int_opt)
% LS channel estimation function
% Inputs:
% Y = Frequency-domain received signal
% Xp = Pilot signal
% pilot_loc = Pilot location
% N = FFT size
% Nps = Pilot spacing
% int_opt = 鈥檒inear鈥� or 鈥檚pline鈥�
% output:
% H_LS = LS Channel estimate
LS_est = Y(pilot_loc)./Xp; % LS channel estimation
if lower(int_opt(1))=='l'
method='linear';
else
method='spline';
end
% Linear/Spline interpolation
H_LS = interpolate(LS_est,pilot_loc,Nfft,method);
end
function H_MMSE = MMSE_CE(Y,Xp,pilot_loc,Nfft,Nps,h,SNR)
% MMSE channel estimation function
% Inputs:
% Y = Frequency-domain received signal
% Xp = Pilot signal
% pilot_loc = Pilot location
% Nfft = FFT size
% Nps = Pilot spacing
% h = Channel impulse response
% SNR = Signal-to-Noise Ratio[dB]
% output:
% H_MMSE = MMSE channel estimate
% Calculate RMS delay spread
Ph = h.*conj(h);
Ptotal = h'*h;
t_sym = 1*(0:length(h)-1)';
t_mean = sum(t_sym.*Ph/Ptotal);
t_cov = sum(t_sym.^2.*Ph/Ptotal);
t_rms = sqrt(t_cov-t_mean^2);
f_max = 100;
H_MMSE = MMSE_ideal(Y,Xp,pilot_loc,Nfft,Nps,SNR, t_rms, f_max);
end
function [H_WOA, Positions] = WOA_CE(Y,Xp,pilot_loc,Nfft,Nps, Nbs, ...
Positions, NumAgent, Max_iter, lb, ub, dim)
% fobj = @CostFunction
% dim = number of your variables
% Max_iteration = maximum number of generations
% SearchAgents_no = number of search agents
% lb = [lb1,lb2,...,lbn] where lbn is the lower bound of variable n
% ub = [ub1,ub2,...,ubn] where ubn is the upper bound of variable n
fobj = @ (x) MMSE_loss(Y, Xp, pilot_loc, Nfft, Nps, Nbs, x(1), x(2), x(3));
x = WhaleOptAlg(NumAgent,Max_iter,lb,ub,dim,fobj);
SNR = x(1);
t_rms = x(2);
f_max = x(3);
H_WOA = MMSE_ideal(Y,Xp,pilot_loc,Nfft,Nps,SNR, t_rms, f_max);
end
function fileId = getFileId(enable)
if enable == 0
fileId = 0;
return
end
ctime = clock;
cmonth = ctime(2); smonth = num2str(cmonth);
cday = ctime(3); sday = num2str(cday);
chour = ctime(4); shour = num2str(chour);
cminute = ctime(5); sminute = num2str(cminute);
fileName = ['CE', smonth, sday, shour, sminute, '.txt'];
fileId = fopen(fileName, 'w');
end
🎉3 参考文献
