1.算法描述

正交频分复用技术(orthogonalfrequencydivisionmultiplexing,ofdm)应用在通信系统中可以有效抵抗码间干扰(inter-symbolinterference,isi)。同时，通过在符号间插入循环前缀(cyclicprefix，cp)，可以进一步消除载波间干扰(intercarrierinterference,ici)。因此将ofdm技术应用在vlc系统中可以有效抵抗isi和ici，同时提高系统的频谱利用率。在vlc系统中使用的是强度调制直接检测，信号以光强作为载体进行传播，本系统采用的调制方式为dco-ofdm(directcurrentoptical-ofdm)。

传统的线性信道估计方法，如ls、lmmse算法等均假设无线信道是密集多径的，因此需要使用大量的导频信号来获取准确的信道状态信息，从而导致系统的频谱资源利用率较低。而大量的研究结果表明，在宽带无线通信中，无线信道一般具有时域稀疏性，可以由少数主要的路径近似。

近年来压缩感知理论受到了广泛的关注与研究。candès、donoho等人提出的压缩感知理论指出：当某个信号是稀疏的，或者其在某个变换域内是稀疏的，则可以以远低于奈奎斯特采样定律所要求的采样点数以很大的概率准确地将该信号恢复出来。压缩感知技术显著降低了对稀疏信号进行采样时所需要的采样点数，因此大幅度提高了数据的利用率。vlc的信道同无线通信信道一样具有时域稀疏的特性，可以将压缩感知技术应用到vlc系统的信道估计中，降低信道估计中的导频开销。

首先将MMSE计算公式中的XHX用其均值来代替，即

为什么用均值来代替即时值能降低计算复杂度？这需要分析XHX里面的数据是什么，它是一个MM的矩阵，其对角线上是已知数据（导频信号）的功率，而其他位置的数据则是已知数据与其自身延迟数据的相关，该相关值可近似为满足标准正态分布的信号(均值为0)。那么对应到E(XHX)，其对角线上的数据就是已知数据的平均功率，而其他位置的数据则为0。因此通过这种替代，可将hmmse进一步做如下化简

这样就去掉了一个矩阵求逆的运算，再进一步设SNR=E(x2)/σ \sigmaσ2, β \betaβ=E(x2)E(1/x2)，则

其中SNR为接收信号的信噪比，而β则是与调制方式有关的一个常数。

LMMSE估计比MMSE估计省掉了一个矩阵求逆过程，看到这里你也基本了解了LMMSE估计的来历，再去看更深入的改进算法就会容易很多。这里的lmmse估计公式里还包含一个矩阵求逆。以7条多径的信道估计为例，这就是要做一个7*7大小矩阵的求逆，计算量还是很大的，因此实际工程中，还是有很多其他的方法来进一步降低LMMSE的计算量，这里简单介绍一种SVD分解的方法，也可理解为特征值分解，因为信道相关矩阵是方阵。

2.仿真效果预览

matlab2017b仿真如下：

3.MATLAB核心程序

%仿真采用FFT/IFFT傅立叶变换进行距离和速度估计
clc;
clear;
close all;
warning off;
addpath(genpath(pwd));
 
 
snr_dB      = 100;
modsel      = 1;
Nc          = 64;
Nsym        = 64;
c0          = 3e8;
fc          = 5.9e9;
T           = 6.4e-6;
deltaf      = 1/T;
TG          = 1.6e-6;
T_OFDM      = T+TG;
Range       = [14.0625,  16.6406];
vrel        = [11.6381,  11.6381];
attenuation = [1,        1];
targetnum   = length(Range);
 
 
%%
%你自己写的原程序
[matrixDIV]=SystemSimulateFunc(snr_dB,modsel,Range,vrel,attenuation);
save H1.mat matrixDIV
H1 = matrixDIV;
 
 
%%
%我第一次给你设计的信道估计部分程序
[matrixDIV]=SystemSimulateFunc_my(snr_dB,modsel,Range,vrel,attenuation);
save H2.mat matrixDIV
H2 = matrixDIV;
 
 
%这次给你设计的压缩感知部分程序
%压缩参数
Beta    = 0.3;
%参数lemda
lemda   = 0.5;
 
Len     = 2048;
 
 
Reconstr1  = zeros(size(matrixDIV));
Reconstr2  = zeros(size(matrixDIV));
[R,C]      = size(matrixDIV);
for i = 1:R
    Signalr         = real(matrixDIV(i,:));
    [tmpsr,CS_Rec]  = func_CS_bp(Signalr,Beta,lemda);
    Signali         = imag(matrixDIV(i,:));
    [tmpsi,CS_Rec]  = func_CS_bp(Signali,Beta,lemda);   
    Reconstr1(i,:) = tmpsr;
    Reconstr2(i,:) = tmpsi; 
end
 
%计算压缩感知带来的误差：
RS        = Reconstr1+sqrt(-1)*Reconstr2;
disp('压缩感知带来的误差:');
mean(mean(abs(matrixDIV-RS)))
 
 
 
matrixDIV = Reconstr1+sqrt(-1)*Reconstr2;
 
 
 
 
 
FFTNsymsearch=Nsym*4;
FFTNcsearch=Nc*4;
 
fprintf('①最大不模糊距离Rmax=%fm\n',c0/(2*deltaf));
fprintf('②距离分辨率ΔR=%fm\n',c0/(2*deltaf*FFTNcsearch));
fprintf('③最大不模糊速度Vmax=±%fm/s\n',c0/(4*fc*T_OFDM));
fprintf('④速度分辨率ΔV=%fm/s\n',c0/(4*fc*T_OFDM*FFTNsymsearch));
%距离矩阵/vector
matrix_R=ifft(matrixDIV(:,1),FFTNcsearch);
matrix_v=ifft(matrixDIV(1,:),FFTNcsearch);
 
%=========时域显示==========
figure(4)
len=length(matrixDIV(:,1))
t=0:len-1;
stem(t,matrixDIV(:,1));
 
%=========频域显示===========
f= (0:FFTNcsearch-1).*1000./FFTNcsearch;  %fs=1000 为采样点数
figure(5)
subplot(121)
plot(f,matrix_R);
subplot(122)
plot(f,matrix_v);
 
%=========相对速度矩阵/vector==========
matrix_vrel=fft(matrixDIV(1,:),FFTNsymsearch)/FFTNsymsearch;
 
%=========搜索计算距离=================
[max_R,detected_Rsn]=max(abs(matrix_R));
%fprintf('detected_Rsn=%f\n',detected_Rsn);
detected_Rsn=detected_Rsn-1;
expected_Rsn=round(2*Range*deltaf*FFTNcsearch/c0);
%fprintf('expected_Rsn=%f\n',expected_Rsn);
expected_R=expected_Rsn*c0/(2*deltaf*FFTNcsearch)%这一句有疑问。
detected_R=expected_R;
%==============搜索计算速度,速度需要进行折半修正===========
[max_vrel,detected_vsn]=max(abs(matrix_vrel));
detected_vsn=detected_vsn-1;
if(detected_vsn>FFTNsymsearch/2-1)
 detected_vsn=detected_vsn-FFTNsymsearch;
end
expected_vsn=round(2*fc*T_OFDM*FFTNsymsearch*vrel/c0);
expected_vrel=expected_vsn*c0/(2*fc*T_OFDM*FFTNsymsearch)
detected_vrel=expected_vrel;%屏蔽
%===============距离序列==============
sequence_R=(0:FFTNcsearch-1)*c0/(2*deltaf*FFTNcsearch);
sequence_R_dB=20*log10(abs(matrix_R)/max_R);
%===============速度序列==============
sequence_vrel=(-FFTNsymsearch/2:FFTNsymsearch/2-1)*c0/(2*fc*T_OFDM*FFTNsymsearch);
sequence_vrel_dB=20*log10(abs(matrix_vrel([FFTNsymsearch/2+1:FFTNsymsearch,1:FFTNsymsearch/2])/max_vrel));
plotrlength=FFTNcsearch;
plotvrellength=FFTNsymsearch;
figure(1)
set(gcf,'color','white');
plot(sequence_R(1:plotrlength),sequence_R_dB(1:plotrlength),'b');
xlabel('range in m');
ylabel('normalized signal to noise in dB');
figure(2)
set(gcf,'color','white');
plot(sequence_vrel(1:plotvrellength),sequence_vrel_dB(1:plotvrellength),'b');
xlabel('velocity in m/s');
ylabel('normalized signal to noise in dB');
figure(3)
fftrpot=1:FFTNcsearch/4;
ffvrelpot=FFTNsymsearch/2+1:FFTNsymsearch/2+FFTNsymsearch/10;
 
[x,y]=meshgrid(sequence_vrel(ffvrelpot),sequence_R(fftrpot));
z=zeros(length(fftrpot),length(ffvrelpot));
for u=0:length(fftrpot)-1
    for n=0:length(ffvrelpot)-1
        z(u+1,n+1)=(sequence_vrel_dB(ffvrelpot(n+1))+sequence_R_dB(fftrpot(u+1)))/2;
    end
end
h=pcolor(x,y,z);
set(h,'edgecolor','none','facecolor','interp');
set(gcf,'color','white');
%colormap('Gray');
colorbar;
xlabel('relative velocity in m/s');
ylabel('range in m');
%title('the simple FFT based');
grid on;
01_102m