这周最紧要任务，写多径效应的CIR估计，故将思路历程收集于一个 blog 中，以供今后参考
首先这个任务题目就看不懂，啥叫 CIR估计？
问了学姐，是信道脉冲响应，突然联想到之前一直听说的信道估计，百度百科对于信道估计的解释是：
所谓信道估计，就是从接收数据中将假定的某个信道模型的模型参数估计出来的过程。如果信道是线性的话，那么信道估计就是对系统冲激响应进行估计。需强调的是信道估计是信道对输入信号影响的一种数学表示，而“好”的信道估计则是使得某种估计误差最小化的估计算法。
2021-9-2更新
什么是CIR估计？
CIR包括两个纬度
CIR就是每个延时点上信号的强度
两个纬度：（1）每路信号的延时时间；（2）每路信号的强度
之前老师让看的两篇文章与CIR估计有什么关系？
【文献精读】【通信】Symphony: Localizing Multiple Acoustic Sources with a Single Microphone Array
【文献笔记】【通信】MAVL: Multiresolution Analysis of Voice Localization
DoA/AoA 的估计就是按信号时延来算的，而且用的也是互相关的的方法，因此在算角度的过程中已经把CIR算出来了
如果是多径信道，令 rx 与 tx 信号作互相关，肯定会在多个延时处有不同信号（这里的延时指tx相对rx的延时，不同路径信号延时不同）
不仅如此，对于多径信道来说，作自相关也能算出CIR，比如之前看的同一个source 的 LOS 和 ECHO两个信道，令 tx 作自相关，就能得到 LOS 和 ECHO 之间的时延
那么也就是说，我的任务可以转述为，对有多径效应的信道作信道估计
这其中有两个问题需要搞清楚：
（1）有多径效应的信道该如何用代码表示？
（2）信道估计算法怎么写？
首先，在开始调研之前，我对信道的理解如下：
一个通信系统有编码、调制、信道传输、解调、解码这几个过程
把要传送的数据编码，为的是增加复杂度，减小出错率
调制，是为了能够将数据嵌入一个carrier中，从而能够实现传输
而通信模型中的信道发挥存在感的就是加入噪声干扰，使解调和解码出现错误
由我当前对通信系统以及信道的理解，可以对题目引出如下几点疑问：
（1）信道脉冲响应的意义是什么？
（2）多径效应的信道和其他效应的信道相比有什么区别？或者说，不同的信道在哪些方面会有差别？
问题（1）隔天再看已经可以解答
学信号与系统的时候就知道，在线性时不变系统中，所有输入对应的输出都可以看成是一堆脉冲响应的线性组合，所以只要知道了这个系统的冲激响应，那么任意给一个输入就能算出输出
第二个问题有点深，在探究这个问题之前，我觉得可以先实现一个普通的信道估计算法来感受一下
8-30更新
通过阅读两篇声源定位的论文，可以初步下一个定义，所谓信道估计，就是得到信道传递函数 H 的过程，这个H就是一个关于衰减与延迟的函数
根据 2 中的描述，为什么要作信道估计？
是为了在接收端进行信道补偿，从而提高整个的系统性能。
《MIMO-OFDM无线通信技术及MATLAB实现》有专门的一章来讲信道估计
信道估计一般使用前导 preamble 或者导频 pilot
导频结构有三种：
（1）块状类型：适用于频率选择性信道
（2）梳状类型：适用于快衰弱信道
（3）格状类型
基于训练符号的信道估计
训练符号估计需要发送前导或者导频信号，会产生额外负荷降低传输效率。
使用训练符号时，一般采用 最小二乘（LS）或者 最小均方误差（MMSE）技术来进行信道估计
LS方法，每个子载波上的LS信道估计可以表示为
其中，Y为输出，X为输入
LS方法会增强噪声，但是因为简单而被广泛使用
接下来的一堆公式推导看不懂
根据与导师的交流，对于多径信道来说 LS 和 MMSE根本不够用，并推荐我看了如下两篇论文：
【文献精读】【通信】Symphony: Localizing Multiple Acoustic Sources with a Single Microphone Array
这两篇论文看完后还是懵逼，因为这两篇paper的主题是声源定位，所以只是简要介绍了多径信道的模型，并没有讲信道估计算法
这篇帖子详细贴了 OFDM系统的信道估计代码
OFDM中信道估计技术分析与实现
我亲手打了一遍这个帖子的代码，有一两个小bug，但是过程非常清晰，亲手 coding 之后感觉酣畅淋漓，我自己写的详细注释版如下：
close all;
clear;
clc;
% 单纯打印在命令行窗口
fprintf('OFDM信道估计仿真\n');
%% 参数设置
carrier_count = 64; % 单个OFDM符号中的载波数目
num_symbol = 50; % OFDM符号数量，不包括导频
guard = 8; % 循环前缀 CP, 即每个OFDM的最后8个子载波作为CP加在开头
pilot_inter = 8; % 导频间隔，即每隔几个OFDM符号插入一个导频
modulation_mode = 16;% 调制方式, QAM调制的order
SNR = 0:2:20; % 信噪比取值
num_loop = 15; % 循环次数
% 每种信噪比做15次实验
num_bit_err = zeros(length(SNR), num_loop);
num_bit_err_dft = zeros(length(SNR), num_loop);
num_bit_err_ls = zeros(length(SNR), num_loop);
MSE = zeros(length(SNR), num_loop);
MSE_dft = zeros(length(SNR), num_loop);
MSE_ls = zeros(length(SNR), num_loop);
%% 主程序
for c1 = 1:length(SNR)
fprintf('\n\n\n仿真信噪比为%f\n\n', SNR(c1));
for num1 = 1:num_loop
%---------------产生发送的随机序列-----------------------
bits_len = carrier_count * num_symbol;
bits_tx = randi([0 1], 1, bits_len); % 1行bits_len列矩阵
%---------------符号调制-----------------------
modulated_sequence = qammod(bits_tx, modulation_mode);
%---------------导频格式-----------------------
pilot_len = carrier_count; % 导频与OFDM符号同级别
% rand 默认区间（0，1）
% round 四舍五入
% 结束后 pilot_symbols只有0和1，是一个1行pilot_len列矩阵
pilot_symbols = round(rand(1, pilot_len));
% 将导频信号中的0变为-1，所以目前导频有两种，-1和1
for i = 1:pilot_len
if pilot_symbols(1, i) == 0
pilot_symbols(1, i) = pilot_symbols(1, i) - 1;
else
pilot_symbols(1, i) = pilot_symbols(1, i);
end
end
% 行向量变列向量
pilot_symbols = pilot_symbols';
%---------------计算导频和数据数目-----------------------
% ceil 向上取整
num_pilot = ceil(num_symbol / pilot_inter);
% rem 取余，返回 num_symbol除以pilot_inter后的余数
% 每组符号开头一个导频，结尾还要再加一个导频
if rem(num_symbol, pilot_inter)==0
num_pilot = num_pilot + 1;
end
num_data = num_symbol + num_pilot;
%---------------导频位置计算-----------------------
pilot_index = zeros(1, num_pilot);
% +1 是算上了导频的位置，此时data_index容量大于 num_symbol
data_index = zeros(1, num_pilot * (pilot_inter + 1));
for i = 1:num_pilot-1
pilot_index(1, i) = (i-1) * (pilot_inter + 1) + 1;
end
% 最后一个导频的位置从开头数起要看具体情况
% 但反正是在最后一个
pilot_index(1, num_pilot) = num_data;
% 为了插入导频，第一个、最后一个以及每隔8个位的索引都空了出来
% 缺少的索引值在 pilot_index
for j = 0:num_pilot
data_index(1, (1 + j * pilot_inter) : (j + 1) * pilot_inter) = (2 + j * (pilot_inter + 1)) : ((j + 1) * (pilot_inter + 1));
end
% 将多余的容量去掉
data_index = data_index(1, 1:num_symbol);
%---------------导频插入-----------------------
% 和之前不一样，每一列为一个OFDM调制符号
piloted_ofdm_syms = zeros(carrier_count, num_data);
% piloted_ofdm_syms 与 modulated_sequence的列数不一样
% piloted_ofdm_syms 的列数增加了 导频数
% 而modulated_sequence的列数只有 OFDM符号数
% 先将 modulated_sequence 中应该是 OFDM符号 的列数依次插入 modulated_sequence 的列
% data_index 和 modulated_sequence 的列数都是 num_symbol
piloted_ofdm_syms(:, data_index) = reshape(modulated_sequence, carrier_count, num_symbol);
% 再将 pilot 插入其索引位
% pilot_symbols 只是一个列向量，横向扩充 num_pilot 倍成为一个数组
piloted_ofdm_syms(:, pilot_index) = repmat(pilot_symbols, 1, num_pilot);
%---------------IFFT变换-----------------------
% OFDM调制后默认是频域信号
% 为什么要乘以 sqrt(carrier_count)
time_signal = sqrt(carrier_count) * ifft(piloted_ofdm_syms);
%---------------加循环前缀-----------------------
% CP 与子载波同级别
add_cyclic_signal = [time_signal((carrier_count - guard + 1 : carrier_count), :); time_signal];
% 将矩阵变为一个行向量再传输
tx_data_trans = reshape(add_cyclic_signal, 1, (carrier_count + guard) * num_data);
%---------------信道处理-----------------------
% AWGN 信道
% 计算信号功率
tx_signal_power = sum(abs(tx_data_trans(:)) .^ 2) / length(tx_data_trans(:));
% 为啥这么算？
noise_var = tx_signal_power / (10 ^ (SNR(c1) / 10));
rx_data = awgn(tx_data_trans, SNR(c1), 'measured');
%---------------信号接收、去循环前缀、FFT变换-----------------------
% 信号接收
rx_signal = reshape(rx_data, (carrier_count + guard), num_data);
% 去除CP
rx_signal_matrix = zeros(carrier_count, num_data);
rx_signal_matrix = rx_signal(guard + 1:end, :);
% FFT变换
% 因为接下来要提取导频，导频自然得回到频域才能提取
rx_carriers = fft(rx_signal_matrix) / sqrt(carrier_count);
%---------------导频和数据提取-----------------------
% 提取导频
rx_pilot = rx_carriers(:, pilot_index);
% 提取数据(频域)
rx_fre_data = rx_carriers(:, data_index);
%---------------导频位置信道响应LS估计-----------------------
% pilot_patt 就是发送端插入的 pilot
pilot_patt = repmat(pilot_symbols, 1, num_pilot);
% LS算法的 y/x
pilot_esti = rx_pilot ./ pilot_patt;
%---------------LS估计的线性插值-----------------------
% 使用信道估计，其实就已经默认了rx的data是不可信的
% 所以需要使用导频来恢复 data
int_len = pilot_index;
len = 1:num_data;
channel_H_ls = zeros(carrier_count, num_data);
% 每个符号/导频的子载波数量是固定的
% 现在是要通过 num_pilot 个 导频来估算 num_data 个导频+data
% 因此子载波一个一个来
for ii = 1:carrier_count
% 在 int_len(pilot_index)这些点（点的数量为 num_pilot）
% 函数取值分别为 pilot_esti(ii, 1:(num_pilot))
% 那么同样的函数，将点扩充为 len(num_data)个，各点的估计值是多少？
channel_H_ls(ii, :) = interp1(int_len, pilot_esti(ii, 1:(num_pilot)), len, 'linear');
end
channel_H_data_ls = channel_H_ls(:, data_index);
%---------------LS估计中发送数据的估计值-----------------------
% 这个公式从何而来还不懂
tx_data_estimate_ls = rx_fre_data .* conj(channel_H_data_ls) ./ (abs(channel_H_data_ls) .^ 2);
%---------------DFT估计------------
% DFT估计是为了再ls基础上消除时域噪声
% 但是以下代码没看懂为什么能达到消除时域噪声的效果
% 或许是在时域作超过子载波长度点数的FFT能够消除噪声？
% 先将 pilot_esti 换为时域
% 补充子载波padding后（symbol纬度不变）
% 再换回频域
% 加padding 其实就是做了 carrier_count+1024 点 FFT 变换
tx_pilot_estimate_ifft = ifft(pilot_esti);
padding_zero = zeros(1024, num_pilot);
tx_pilot_estimate_ifft_padding_zero = [tx_pilot_estimate_ifft; padding_zero];
% 对于矩阵来说，fft分别对各列作fft
tx_pilot_estimate_dft = fft(tx_pilot_estimate_ifft_padding_zero);
%---------------DFT估计的线性插值------------
% 插值操作与 ls 一样
int_len = pilot_index;
len = 1:num_data;
channel_H_dft = zeros(carrier_count, num_data);
for ii = 1:carrier_count
channel_H_dft(ii, :) = interp1(int_len, tx_pilot_estimate_dft(ii, 1:(num_pilot)), len, 'linear');
end
channel_H_data_dft = channel_H_dft(:, data_index);
%---------------DFT估计中发送数据的估计值------------
tx_data_estimate_dft = rx_fre_data .* conj(channel_H_data_dft) ./ (abs(channel_H_data_dft) .^ 2);
%---------------DFT符号解调------------
% 原矩阵每一列向量变为行向量
% 并且按顺序铺成 sequence
demod_in_dft = tx_data_estimate_dft(:).';
demod_out_dft = qamdemod(demod_in_dft, modulation_mode);
%---------------LS符号解调------------
demod_in_ls = tx_data_estimate_ls(:).';
demod_out_ls = qamdemod(demod_in_ls, modulation_mode);
%---------------误码率计算------------
for i = 1:length(bits_tx)
if demod_out_dft(i) ~= bits_tx(i)
num_bit_err_dft(c1, num1) = num_bit_err_dft(c1, num1) + 1;
end
if demod_out_ls(i) ~= bits_tx(i)
num_bit_err_ls(c1, num1) = num_bit_err_ls(c1, num1) + 1;
end
end
end
end
% num_bit_err_dft.'表示原来的行变列
% 原先是每行代表一个 SNR 的所有 loop 结果
% 现在变成每列
% 与num_bit_err_dft(:).'展开不同，num_bit_err_dft.'仍是矩阵
% 因为 mean 是按列来算平均值的
BER_dft = mean(num_bit_err_dft.') / length(bits_tx);
BER_ls = mean(num_bit_err_ls.') / length(bits_tx);
%% 绘图
figure
semilogy(SNR, BER_dft, '-mp', SNR, BER_ls, '-k+');
title('OFDM系统的LS和DFT信道估计');
xlabel('SNR');
ylabel('BER');
legend('DFT信道估计','LS信道估计');
通过这次coding 主要明白了以下的几个重要知识点：
（1）pilot, CP 的真面目
导频 pilot 与 符号 symbol 同级别 ，每隔几个 symbol 插入一个 pilot
更新，pilot也可以与子载波同级别，导频有很多种插入法
循环前缀 CP 与 子载波 同级别，每个symbol的最后n个子载波复制了然后插入到开头作为CP
比如一个 symbol [1, 2, 3, 4, 5], 假设 CP 为2，那么加入CP后变为 [4, 5, 1, 2, 3, 4, 5]
（2）明白了pilot的真正用途
使用信道估计，其实就已经默认了rx的data是不可信的，所以需要使用不易被干扰的导频来获取信道对信号的影响，即得到H，然后再通过 rx与 H 的运算来获得可信 data
（3）DFT估计，就是在LS基础上将导频ifft换回时域然后再作N点fft换回频域（N大于载波数量）
并且，通过此次coding 明确了我搞不懂的问题到底在哪里
要做XXX信道估计，要解决两个问题
（1）建立这个信道的模型
加多少噪音、衰弱、延时？以什么方式加？
（2）选择合适的估计算法
如果 y = hx + n，那就 ls, mmse, dft
卡壳点：
（1）微信论文中有给出多径信道模型，代码如何实现？
看起来还是 y = hx + n的形式，为什么 ls 不行了？
（2）不知道该上哪找合适的估计算法
[1] 移动无线信道脉冲响应的计算机模拟
[2] Matlab关于信道估计的一个描述