Code Matlab

Orthogonal Matching Pursuit Algorithm (OMP)

clear allclose allclc
M = 8;J = 10;Q = 2;L = 4;P = M + L;N=M+Q+L-1;Mq=M+Q-1;Pq=N;%%%%%%%%%%%%%              VpT = [];%%%%%%%%%%%%%%%%%%55    UU = [];
iter0 = 100;iter1 = 10;ber1 = zeros(1,iter1);ber0 = zeros(1,iter0);snr = 10:5:45;ber = zeros(1,length(snr));for m = 1:length(snr)    m     for k = 1:iter0        h1 = randn(L+1,iter1)+1j*randn(L+1,iter1);%ok        U1 = 2*(randn(M,iter1*J)>0) - 1 + 1j*(2*(randn(M,iter1*J)>0) - 1);        for ell = 1:iter1            h = h1(:,ell);            U = U1(:,(ell - 1)*J+1:ell*J);            eb = sum(sum(abs(U).^2))/(2*M*J);            att = sqrt(eb/2)*10^(-snr(m)/20);%eb/2           H = FBToeplitz(h,M);            Hf = FBToeplitz(h,Mq);            YJ = H*U + att*(randn(P,J) + 1j*randn(P,J));            YJQ = zeros(P+Q-1,J*Q);            for n = 1:J                YJQ(:,(n-1)*Q+1:n*Q) = FBToeplitz(YJ(:,n),Q);            end            Yio=[1 2 2 2 ]            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     [hf,UoT,Uof] = ZPblindf(YJQ,Hf,L,N,Mq,Pq);          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     r=3;    %%%% UoT es U primera    %%% Uo son todas las Ui    % [x] = OMP (r,Y,A)     % [x] = OMP (r,UoT,Uof)                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           % [AVf,hf]             Ufft = fft(YJ,Pq);             hfft=hf*ones(1,J);              Ur = ifft(Ufft./hfft);            Ur1 = Ur(1:M,:);            Ur1r = 2*(real(Ur1)>0) - 1;            Ur1i = 2*(imag(Ur1)>0) - 1;            Uur = real(U);            Uui = imag(U);            ber1(ell) = (sum(sum(abs(Ur1r-Uur))) + sum(sum(abs(Ur1i-Uui))))/(4*M*J);            %pause        end        ber0(k) = sum(ber1)/iter1;    end    ber(m) = sum(ber0)/iter0;endsemilogy(snr,ber)

clear allclose allclc
M = 8;J = 10;Q = 2;L = 4;P = M + L;N=M+Q+L-1;Mq=M+Q-1;Pq=N;%%%%%%%%%%%%%              VpT = [];%%%%%%%%%%%%%%%%%%55    UU = [];
iter0 = 100;iter1 = 10;ber1 = zeros(1,iter1);ber0 = zeros(1,iter0);snr = 10:5:45;ber = zeros(1,length(snr));for m = 1:length(snr)    m     for k = 1:iter0        h1 = randn(L+1,iter1)+1j*randn(L+1,iter1);%ok        U1 = 2*(randn(M,iter1*J)>0) - 1 + 1j*(2*(randn(M,iter1*J)>0) - 1);        for ell = 1:iter1            h = h1(:,ell);            U = U1(:,(ell - 1)*J+1:ell*J);            eb = sum(sum(abs(U).^2))/(2*M*J);            att = sqrt(eb/2)*10^(-snr(m)/20);%eb/2           H = FBToeplitz(h,M);            Hf = FBToeplitz(h,Mq);            YJ = H*U + att*(randn(P,J) + 1j*randn(P,J));            YJQ = zeros(P+Q-1,J*Q);            for n = 1:J                YJQ(:,(n-1)*Q+1:n*Q) = FBToeplitz(YJ(:,n),Q);            end            Yio=[1 2 2 2 ]            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     [hf,UoT,Uof,AVf] = ZPblindf(YJQ,Hf,L,N,Mq,Pq);     %as=AVF;     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    % ff=[as(1,1); 0; 0; 0; as(2,1) 0 0 0 h(3,1) 0 0 0 h(4,1) 0 0 0 h(5,1)];    % r=5;    %%%% UoT es U primera    %%% Uo son todas las Ui    % [x] = OMP (r,Y,A)    % [x] = OMP (r,ff,Uof)                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           % [AVf,hf]             Ufft = fft(YJ,Pq);             hfft=hf*ones(1,J);              Ur = ifft(Ufft./hfft);            Ur1 = Ur(1:M,:);            Ur1r = 2*(real(Ur1)>0) - 1;            Ur1i = 2*(imag(Ur1)>0) - 1;            Uur = real(U);            Uui = imag(U);            ber1(ell) = (sum(sum(abs(Ur1r-Uur))) + sum(sum(abs(Ur1i-Uui))))/(4*M*J);            %pause        end        ber0(k) = sum(ber1)/iter1;    end    ber(m) = sum(ber0)/iter0;endsemilogy(snr,ber)

function U = CS_Hankel(R)%UNTITLED2 Summary of this function goes here% Detailed explanation goes hereL = size(R);m = L(1) - L(2);U1 = zeros(L(2)+1,m); for n = 1:L(2)    r = R(:,n).';    for ell = 1:L(2)+1        U1(ell,:) = r(ell:m+ell-1);    end    U(:,(n-1)*m+1:n*m) = U1;endend

clear allclose allclc
M = 8;J = 2;Q = 8;L = 4;P = M + L;N=M+Q+L-1;Pq=N;Mq=M+Q-1;
h = randn(L+1,1)+1j*randn(L+1,1);
H = FBToeplitz(h,M);
U = 2*(randn(M,J)>0) - 1 + 1j*(2*(randn(M,J)>0) - 1);%%%%%%%%%%%%%%%%%%%%%dominio de tiempo%%%%%%%%%%%%%%%%%YJ = H*U;YJQ = zeros(P+Q-1,J*Q);for n = 1:J    YJQ(:,(n-1)*Q+1:n*Q) = FBToeplitz(YJ(:,n),Q);end
[ht,Rt,Ut] = ZPblind(YJQ,L);
[hm,I] = max(abs(h));ht = ht*h(I)/ht(I);[h,ht]%%%%%%%%%%%%%%%%%%dominio de frecuencia%%%%%%%%%%%%%%%%%%%

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           % En dominio de frecuencia                                 VNxPQ = zeros(N,Pq);                      VpT = [];                      for k=0:N-1               coe = exp((2*pi*1i*k)/N);               VpT = [VpT coe];           end                      VpTf = VpT;                      for i=0:Pq-1               VNxPQ(:,i+1) = VpTf .^-i;           end;                      VNxMQ = VNxPQ(:,1:Mq);                      AV = zeros(N,N);                     Hf = FBToeplitz(h,Mq);            for i=1:N              AV(i,i) = VNxPQ(i,:) * Hf(:,1);             end                        YJf = YJQ               Zf = VNxPQ * YJf;                              lf = size(Zf);               mf = lf(1) - L;                           [S1 V1 D] = svd(Zf);                       S = S1(:,length(S1)-(N-Mq)+1:length(S1));
                      UoT = S';
           [x,y] = size(UoT);
           UU = [];
           VNxMQTranspuesta = VNxMQ.';
           for i=1:x               Ufila = UoT(i,:);
               for j=1:Mq
                   Ui(j,:) = Ufila;
               end
               Ui = Ui .* VNxMQTranspuesta;
               UU = [UU Ui.'];%%%%%%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
           end                       U = UU.';                      UAux = U' * U;                     [EigenVectors EigenValues] = eigs(UAux,3,'SM');           hf = EigenVectors(:,1);           [hmf,If] = max(abs(diag(AV)));           AVf=diag(AV);            hf = 1*hf*AVf(If)/hf(If);            [AVf,hf]           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                                                        %diagEigenValues = diag(EigenValues);                      %minValue = min(diagEigenValues);                     % indexMinValue = -1;                      %for h=1:length(EigenValues)            % if diagEigenValues(h) == minValue             % indexMinValue = h;              % end           %end                      %%%%%%%%%%%%%%%%%%%%%%% AV max Value and max Index          % diagAV = diag(AV);                      %maxValue = max(diagAV);                      %indexMaxValue = -1;                      %for h=1:length(AV)            % if diagAV(h) == maxValue             % indexMaxValue = h;              % end          % end                                 %colEigenVectors = EigenVectors(:,indexMinValue);                      %coeficienteEigenVec = colEigenVectors(indexMaxValue);                      %divCoe = maxValue / coeficienteEigenVec;                      %estimationH = colEigenVectors * divCoe;           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                       %[diagAV, hf]            a=2;

clear allclose allclc
M = 8;J = 2;Q = 8;L = 4;P = M + L;
iter0 = 100;iter1 = 10;ber1 = zeros(1,iter1);ber0 = zeros(1,iter0);snr = 10:5:45;ber = zeros(1,length(snr));for m = 1:length(snr)    m    for k = 1:iter0        h1 = randn(L+1,iter1)+1j*randn(L+1,iter1);        U1 = 2*(randn(M,iter1*J)>0) - 1 + 1j*(2*(randn(M,iter1*J)>0) - 1);        for ell = 1:iter1            h = h1(:,ell);            U = U1(:,(ell - 1)*J+1:ell*J);            eb = sum(sum(abs(U).^2))/(2*M*J);            att = sqrt(eb/2)*10^(-snr(m)/20);            H = FBToeplitz(h,M);            YJ = H*U + att*(randn(P,J) + 1j*randn(P,J));            YJQ = zeros(P+Q-1,J*Q);            for n = 1:J                YJQ(:,(n-1)*Q+1:n*Q) = FBToeplitz(YJ(:,n),Q);            end            [ht,Rt,Ut] = ZPblind(YJQ,L);            [hm,I] = max(abs(h));            ht = ht*h(I)/ht(I);            %[h ht]            Ufft = fft(YJ);            hfft = fft(h,P+Q-1)*ones(1,J);            Ur = ifft(Ufft./hfft);            Ur1 = Ur(1:M,:);            Ur1r = 2*(real(Ur1)>0) - 1;            Ur1i = 2*(imag(Ur1)>0) - 1;            Uur = real(U);            Uui = imag(U);            ber1(ell) = (sum(sum(abs(Ur1r-Uur))) + sum(sum(abs(Ur1i-Uui))))/(4*M*J);            %pause        end        ber0(k) = sum(ber1)/iter1;    end    ber(m) = sum(ber0)/iter0;endsemilogy(snr,ber)

clear allclose allclc
M = 8;J = 10;Q = 2;L = 4;P = M + L;
iter0 = 100;iter1 = 40;ber1 = zeros(1,iter1);ber0 = zeros(1,iter0);snr = 10:5:45;ber = zeros(1,length(snr));for m = 1:length(snr)    m    for k = 1:iter0        h1 = randn(L+1,iter1)+1j*randn(L+1,iter1);        U1 = 2*(randn(M,iter1*J)>0) - 1 + 1j*(2*(randn(M,iter1*J)>0) - 1);        for ell = 1:iter1            h = h1(:,ell);            U = U1(:,(ell - 1)*J+1:ell*J);            eb = sum(sum(abs(U).^2))/(2*M*J);            att = sqrt(eb/2)*10^(-snr(m)/10);%eb/2            H = FBToeplitz(h,M);            YJ = H*U + att*(randn(P,J) + 1j*randn(P,J));            YJQ = zeros(P+Q-1,J*Q);            for n = 1:J                YJQ(:,(n-1)*Q+1:n*Q) = FBToeplitz(YJ(:,n),Q);            end            [ht,Rt,Ut] = ZPblind(YJQ,L);            [hm,I] = max(abs(h));            ht = ht*h(I)/ht(I);            %[h ht]            Ufft = fft(YJ,P+Q-1);            hfft = fft(h,P+Q-1)*ones(1,J);            Ur = ifft(Ufft./hfft);            Ur1 = Ur(1:M,:);            Ur1r = 2*(real(Ur1)>0) - 1;            Ur1i = 2*(imag(Ur1)>0) - 1;            Uur = real(U);            Uui = imag(U);            ber1(ell) = (sum(sum(abs(Ur1r-Uur))) + sum(sum(abs(Ur1i-Uui))))/(4*M*J);            %pause        end        ber0(k) = sum(ber1)/iter1;    end    ber(m) = sum(ber0)/iter0;endsemilogy(snr,ber)

function T = FBToeplitz(v,k)%UNTITLED Summary of this function goes here% Detailed explanation goes hereT = zeros(length(v)+k-1,k);v = v(:);for ell = 1:k    T(:,ell) = [zeros(ell-1,1);v;zeros(k-ell,1)];endend

function [x] = OMP (K,y,A)%Orthogonal Matching Pursuit Algorithm (OMP) is a greedy compressed sensing%recovery algorithm which selects the best fitting column of the sensing%matrix in each iteration. A least squares (LS) optimization is then %performed in the subspace spanned by all previously picked columns. Res = y.' ;[m,n] = size (A) ;
Q = zeros (m,K) ;R = zeros (K,K) ;Rinv = zeros (K,K) ;w = zeros (m,K) ;x = zeros (1,n) ;
for J = 1 : K        %Index Search    [V ,kkk] = max(abs(A'*Res)) ;    kk (J) = kkk ;        %Residual Update    w (:,J) = A (:,kk (J)) ;    for I = 1 : J-1        if (J-1 ~= 0)            R (I,J) = Q (:,I)' * w (:,J) ;            w (:,J) = w (:,J) - R (I,J) * Q (:,I) ;        end    end    R (J,J) = norm (w (:,J)) ;    Q (:,J) = w (:,J) / R (J,J) ;    Res = Res - (Q (:,J) * Q (:,J)' * Res) ;  end
%Least Squaresfor J = 1 : K    Rinv (J,J) = 1 / R (J,J) ;    if (J-1 ~= 0)        for I = 1 : J-1            Rinv (I,J) = -Rinv (J,J) * (Rinv (I,1:J-1) * R (1:J-1,J)) ;        end    endend
xx = Rinv * Q' * y.' ;
for I = 1 : K    x (kk (I)) = xx (I) ;end

function [h,R,U] = ZPblind(Z,L)%UNTITLED3 Summary of this function goes here% Detailed explanation goes herel = size(Z);m = l(1) - L;[R1,S,V] = svd(Z);R = R1(:,m+1:l(1));U = CS_Hankel(R);U1 = U*U';[h1,D] = eigs(U1,3,'SM');h = h1(:,1);end

function [hf,UoT,Uo,AVf] = ZPblindf(YJf,Hf,L,N,Mq,Pq)             VpT = [];            UU=[];           for k=0:N-1               coe = exp((2*pi*1i*k)/N);               VpT = [VpT coe];           end                      VpTf = VpT;           VNxPQ = zeros(N,Pq);           for i=0:Pq-1               VNxPQ(:,i+1) = VpTf .^(-i);           end;                      VNxMQ = VNxPQ(:,1:Mq);             AV = zeros(N,N);                    % Hf = FBToeplitz(h,Mq);            for i=1:N              AV(i,i) = VNxPQ(i,:) * Hf(:,1);             end                                       Zf = VNxPQ * YJf;                              lf = size(Zf);               mf = lf(1) - L;                           [S1 V1 D] = svd(Zf);                       S = S1(:,length(S1)-(N-Mq)+1:length(S1));            UoT = S';
           [x,y] = size(UoT);              VNxMQTranspuesta = VNxMQ.';
           for i=1:x               Ufila = UoT(i,:);
               for j=1:Mq
                   Ui(j,:) = Ufila;
               end
               Ui = Ui .* VNxMQTranspuesta;
               UU = [UU Ui.'];%%%%%%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
           end                       Uo = UU.';                      UAux = Uo' * Uo;                     [EigenVectors EigenValues] = eigs(UAux,3,'SM');           hf = EigenVectors(:,1);           [hmf,If] = max(abs(diag(AV)));           AVf=diag(AV);            hf = 1*hf*AVf(If)/hf(If);

Google Sites

Report abuse