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);