channel_order_est_main.m

Subspace Projection Based Blind Channel Order Estimation of MIMO Systems m file for a classical

%=======================================
% Channel Order Estimation
% Example 1: 2 inputs 6 outputs MIMO
%=======================================
clear all
clc

%---------------------------------
% Parameter initialization
%---------------------------------
Mt = 2;         % inputs
Mr = 6;         % outputs
N_max = 100;    % Monte Carlo runs
L = 10;         % upper bound of the system's orders
w = 10;         % predictor's size
L1 = 3;         % order of subsystem 1
L2 = 5;         % order of subsystem 2
Q = 400;        % sample number
snr = [17 20 22 25 27 30 32 40 50];       % Signal-to-noise ratio

err = zeros(1,length(snr));
for s = 1:length(snr)
    SNR = snr(s)
    h = waitbar(0,'Please wait...');
    for n = 1:N_max
        waitbar(n/N_max)
        %---------------------------------
        % Generate the source
        %---------------------------------
        Smybol_Len = Q+L;      % symbol number
        Src = randint(Mt,2*Smybol_Len);  % binary source
        mod_bits = tx_QPSK(Src);    % QPSK modulation
        
        %---------------------------------
        % Transmit through MIMO Channel
        %---------------------------------
        Rec_Sig = MIMO_channel(mod_bits,Mr,L,Q,SNR);
        
        %--------------------
        % Estimate order
        %--------------------
        E = cell(1,L+1);
        rank = zeros(1,L+1);
        for l = 0:L
            E{l+1} = form_E_mat(Rec_Sig,l,w,Q);
            [U,S,V] = svd(E{l+1},'econ');
            rank(l+1) = effective_rank(S);
        end

        %--------------------
        % Step 4: Find J1
        %--------------------
        delta = zeros(1,L);
        for l = 1:L
            delta(l) = norm(E{l},inf)/norm(E{l+1},inf);
        end
        [C,I] = min(delta);
        J1 = I;
        m1 = rank(J1+1);
        J2 = J1+1;
        while J2 <= L
            if rank(J2+1)-rank(J2) == m1
                J2 = J2 + 1;
            else
                break;
            end
        end
        err(s) = err(s)+((J1 ~= min(L1,L2))||(J2 ~= max(L1,L2)));
    end  % end N
    close(h)
end

err = err/N_max

save('err_newdeno_400.mat','snr','err')