RLS代码

chianelxie

function [OutHard, OutSoft] = RLS(SynSignal, TS_Symbol, VectorC)
global gcDataSource_bak;

% 群时延
K0 = 3;
% 前馈滤波器阶数
Lf = 4;
% 反馈滤波器阶数
Lb = 2;
nState = 2^Lb;
% 输入信号倍速
SampleRate = size(SynSignal, 1);
% 滤波器初始化
f_k = zeros(SampleRate, Lf);
b_k = zeros(1, Lb);
% MLSE初始化
BurstLen = 85;  
Path_TB = zeros(BurstLen+1, nState);
GMSK_Metric = zeros(1, nState*2);
old_p = 0;
new_p0 = nState;
new_p1 = nState*1.5;
% RLS相关初始化
Delta = 0.004; % 0.004
P_k = 1/Delta * eye(Lf*SampleRate+Lb);
Lamda = 1; % 0 < Lamda <= 1
TrainingModeLen = 40;

for ii = 0:BurstLen
    r_k = SynSignal(:,62+Lb+K0+ii:-1:63+Lb+K0-Lf+ii);
    if ii <= TrainingModeLen
        a_k = TS_Symbol(Lb+ii:-1:1+ii);
        a_k_k0 = TS_Symbol(Lb+1+ii);   
        for jj = 1:SampleRate
            s_k_Temp(jj) = real(sum(r_k(jj,:) .* f_k(jj,:)) * conj(VectorC));
        end
        Total_s_k_Temp = sum(s_k_Temp);
        v_k_Temp = Total_s_k_Temp - b_k * a_k';
        e_k_Temp = v_k_Temp - a_k_k0;

        % RLS
        if SampleRate == 2
            u_k = conj([r_k(1,:), r_k(2,:), (-1)*a_k]');
        else
            u_k = conj([r_k(1,:), (-1)*a_k]');
        end
        pai_k = P_k * u_k;
        k_k = pai_k / (Lamda + (u_k)' * pai_k);
        k_k_Temp = conj((k_k)');
        for jj = 1:SampleRate
            f_k(jj,:) = conj(conj(f_k(jj,:)) - k_k_Temp((jj-1)*Lf+1:jj*Lf) * e_k_Temp);
        end
     %   b_k = b_k - k_k_Temp(end-Lb+1:end) * conj(e_k_Temp);
        b_k = b_k - real(k_k_Temp(end-Lb+1:end)) * e_k_Temp;
        P_k = P_k / Lamda - k_k * (u_k)' * P_k / Lamda;
     %   s_k = Total_s_k_Temp;
    end
   
    % 应用更新后的前、反馈系数做信号投影
    for jj = 1:SampleRate
        s_k(jj) = real(sum(r_k(jj,:) .* (f_k(jj,:))) * conj(VectorC));
    end
    s_k = sum(s_k);
  
    % 2^Lb状态 MLSE 原始分支度量公式
    if nState == 16
        S_L_AK(8) = b_k(1) + b_k(2) + b_k(3);    % 111
        S_L_AK(7) = b_k(1) + b_k(2) - b_k(3);    % 110
        S_L_AK(6) = b_k(1) - b_k(2) + b_k(3);    % 101
        S_L_AK(5) = b_k(1) - b_k(2) - b_k(3);    % 100
        for mm = 1:4
            S_L_AK(mm) = -S_L_AK(8+1-mm);
        end
    elseif nState == 8
        S_L_AK(4) = b_k(1) + b_k(2);             % 11
        S_L_AK(3) = b_k(1) - b_k(2);             % 10
        S_L_AK(2) = (-1) * S_L_AK(3);            % 01
        S_L_AK(1) = (-1) * S_L_AK(4);            % 00
    elseif nState == 4
        S_L_AK(2) = b_k(1);                      % 1
        S_L_AK(1) = (-1) * S_L_AK(2);            % 0
    else
        disp(strcat('MIC-MLSE Equalizer do not support state',num2str(nState)));
    end
    % ACS
    for jj = 1:nState/2
        meas_up0 = GMSK_Metric(old_p+1) + abs(s_k - (-1) - S_L_AK(jj) - (-1)*b_k(Lb))^2;        % 输入 0   输出 0
        meas_down0 = GMSK_Metric(old_p+1) + abs(s_k - 1 - S_L_AK(jj) - (-1)*b_k(Lb))^2;         % 输入 1   输出 0
        meas_up1 = GMSK_Metric(old_p+2) + abs(s_k - (-1) - S_L_AK(jj) - b_k(Lb))^2;             % 输入 0   输出 1
        meas_down1 = GMSK_Metric(old_p+2) + abs(s_k - 1 - S_L_AK(jj) - b_k(Lb))^2;              % 输入 1   输出 1

        SOVA_GMSK_Info(ii+1,jj*2-1) = meas_up1 - meas_up0;
        if (meas_up0 > meas_up1)
            GMSK_Metric(new_p0+1) = meas_up1;
            Path_TB(ii+1, 2*jj-1) = 1;
        else
            GMSK_Metric(new_p0+1) = meas_up0;
            Path_TB(ii+1, 2*jj-1) = 0;
        end

        SOVA_GMSK_Info(ii+1,jj*2) = meas_down1 - meas_down0;
        if meas_down0 > meas_down1                        
            GMSK_Metric(new_p1+1) = meas_down1;
            Path_TB(ii+1, 2*jj) = 1;
        else                                             
            GMSK_Metric(new_p1+1) = meas_down0;
            Path_TB(ii+1, 2*jj) = 0;
        end
        old_p = old_p + 2;
        new_p0 = new_p0 + 1;
        new_p1 = new_p1 + 1;
    end
    if old_p == nState*2
        old_p = 0;
        new_p0 = nState;
        new_p1 = 1.5*nState;
    else
        new_p0 = 0;
        new_p1 = nState/2;
    end   
end

% 回溯
[TempValue, MaxMetricState] = min(GMSK_Metric(1:nState));
for ii = BurstLen+1:-1:1
    MSB = floor(MaxMetricState/(nState/2 + 1));                                
    Demod_Out_hard_Temp(ii) = MSB;
    Demod_Out_soft_Temp(ii) = SOVA_GMSK_Info(ii,2*(mod(MaxMetricState-1, nState/2)) + MSB + 1);
    New_LSB = Path_TB(ii, 2*(mod(MaxMetricState-1, nState/2)) + MSB + 1);
    MaxMetricState = 2*(mod(MaxMetricState-1, nState/2)) + New_LSB + 1;
end
% 输出Burst右侧数据软、硬判结果
OutHard = Demod_Out_hard_Temp(1+26-Lb:58+26-Lb);
OutSoft = Demod_Out_soft_Temp(1+26-Lb+Lb:58+26-Lb+Lb);

% 测试变量
result1 = OutHard - gcDataSource_bak(59:116);
%result3 = DetectorOut(1+26-Lb:58+26-Lb) - gcDataSource_bak(59:116);
for ii = 1:length(OutSoft)
    if (OutSoft(ii) >= 0)
        TmpSoft(ii) = 0;
    else
        TmpSoft(ii) = 1;
    end
end
result2 = TmpSoft - gcDataSource_bak(59:116);
result1;