--- /dev/null
+clc;
+
+%SNR_vals = 20:3:60; N_SNR = length(SNR_vals);
+SNR_vals = -30:0.5:10; N_SNR = length(SNR_vals);
+
+N = 100;
+M = 10;
+
+BWNN = 2*asind(2/M);
+BWR = 2*asind(0.891/M);
+
+synth_noise = 0;
+use_symm = 1;
+
+DOARMSE_idec = zeros(2, N_SNR);
+DOARMSE_sa = zeros(2, N_SNR);
+DOARMSE_toep = zeros(2, N_SNR);
+
+avgDOARMSE_idec = zeros(1, N_SNR);
+avgDOARMSE_sa = zeros(1, N_SNR);
+avgDOARMSE_toep = zeros(1, N_SNR);
+
+POR_idec = zeros(1, N_SNR);
+POR_sa = zeros(1, N_SNR);
+POR_toep = zeros(1, N_SNR);
+
+DOA_idec = zeros(2, N_SNR);
+DOA_sa = zeros(2, N_SNR);
+DOA_toep = zeros(2, N_SNR);
+
+N_runs = 10;
+
+f1 = figure;
+
+for n_run=1:N_runs
+for n_SNR=1:N_SNR
+ if n_run > 1
+ plotix_max = N_SNR;
+ else
+ plotix_max = n_SNR;
+ end
+
+SNR1 = SNR_vals(n_SNR);
+
+disp(['SNR=', num2str(SNR1), 'dB']);
+
+SNR2 = SNR1;
+SNR3 = SNR2;
+theta1 = BWR;
+theta2 = 0;
+theta3 = -31.0;
+gamma = 45*0;
+d = 0.5;
+k = 2*pi;
+T = 0.0;
+omega1 = 2*pi;
+omega2 = 2*pi;
+omega3 = 2*pi;
+
+randn('seed',0);
+
+noise = randn(M,N) + j*randn(M,N);
+
+amp1 = sqrt(2)*10^(SNR1/20);
+amp2 = sqrt(2)*10^(SNR2/20);
+amp3 = 0*sqrt(2)*10^(SNR3/20);
+
+phi1 = -k*d*sin(theta1*pi/180);
+phi2 = -k*d*sin(theta2*pi/180);
+phi3 = -k*d*sin(theta3*pi/180);
+
+if use_symm
+ a1 = exp(j*phi1).^[-(M-1)/2:(M-1)/2]';
+ a2 = exp(j*phi2).^[-(M-1)/2:(M-1)/2]';
+ a3 = exp(j*phi3).^[-(M-1)/2:(M-1)/2]';
+else
+ a1 = exp(j*phi1).^[0:M-1]';
+ a2 = exp(j*phi2).^[0:M-1]';
+ a3 = exp(j*phi3).^[0:M-1]';
+end
+
+A = [a1 a2 a3];
+
+phase1 = j*2*pi*rand(N,1);
+phase2 = j*2*pi*rand(N,1);
+phase3 = j*2*pi*rand(N,1);
+
+n_corr = 1;
+if n_corr
+ %phase1 = zeros(N,1);
+ phase2 = phase1;
+ phase3 = phase2;
+end
+
+signal1 = amp1*exp(phase1);
+signal2 = amp2*exp(phase2);
+signal3 = amp3*exp(phase3);
+
+s = [signal1 signal2 signal3]';
+
+x = A*s + (1-synth_noise)*noise;
+
+delta = 0;
+
+N_mean = 100*M;
+L = M;
+K = M-L+1;
+
+L_sa = M-1;
+%L_sa=round(M/2);
+%L_sa=M;
+K_sa = M-L_sa+1;
+
+x_ = x;
+N_it = 20;
+
+int_f = @arrayest;
+
+R_scm = x*x'/N;
+R_toep = zeros(1,L);
+
+for lix=0:L-1
+ R_toep(lix+1) = mean(diag(R_scm,lix));
+end
+R_toep = toeplitz(R_toep);
+
+R_sa = zeros(L_sa);
+for kix=1:K_sa
+ R_sa = R_sa + x(kix:kix+L_sa-1,:)*x(kix:kix+L_sa-1,:)'/K_sa/N;
+end
+
+R = getDecorrelatedMatrix(x,L,N_it,int_f);
+
+[V,D] = eig(R);
+[V_toep,D_toep] = eig(R_toep);
+[V_scm,D_scm] = eig(R_scm);
+[V_sa,D_sa] = eig(R_sa);
+
+N_sig = 2;
+
+D = diag(D); [D_vals,D_ix]=sort(D); D_ix = D_ix(1:end-N_sig);
+D_toep = diag(D_toep); [D_toep_vals,D_toep_ix]=sort(D_toep);D_toep_ix = D_toep_ix(1:end-N_sig);
+D_sa = diag(D_sa); [D_sa_vals,D_sa_ix]=sort(D_sa); D_sa_ix = D_sa_ix(1:end-N_sig);
+D_scm = diag(D_scm); [D_scm_vals,D_scm_ix]=sort(D_scm);D_scm_ix = D_scm_ix(1:end-N_sig);
+
+V = V(:,D_ix);
+R_n = V*V';
+
+V_scm = V_scm(:,D_scm_ix);
+R_scm_n = V_scm*V_scm';
+
+V_toep = V_toep(:,D_toep_ix);
+R_toep_n = V_toep*V_toep';
+
+V_sa = V_sa(:,D_sa_ix);
+R_sa_n = V_sa*V_sa';
+
+% Calculate DOAs using root algorithms
+N_roots = 2*L-1;
+polcoeffs_idec = zeros(1,N_roots);
+polcoeffs_scm = zeros(1,N_roots);
+polcoeffs_toep = zeros(1,N_roots);
+polcoeffs_sa = zeros(1,N_roots);
+n_diag = -(L-1);
+
+for n_roots=1:N_roots
+ polcoeffs_idec(n_roots) = sum(diag(R_n,n_diag));
+ polcoeffs_scm(n_roots) = sum(diag(R_scm_n,n_diag));
+ polcoeffs_toep(n_roots) = sum(diag(R_toep_n,n_diag));
+ polcoeffs_sa(n_roots) = sum(diag(R_sa_n,n_diag));
+
+ n_diag = n_diag+1;
+end
+
+roots_idec = getMUSICangles(polcoeffs_idec, N_sig);
+roots_scm = getMUSICangles(polcoeffs_scm, N_sig);
+roots_toep = getMUSICangles(polcoeffs_toep, N_sig);
+roots_sa = getMUSICangles(polcoeffs_sa, N_sig);
+
+if (abs(roots_idec(1) + theta1).^2/BWNN < 1) && (abs(roots_idec(2) + theta2).^2/BWNN < 1)
+ POR_idec(n_SNR) = POR_idec(n_SNR) + 1;
+end
+
+if (abs(roots_sa(1) + theta1).^2/BWNN < 1) && (abs(roots_sa(2) + theta2).^2/BWNN < 1)
+ POR_sa(n_SNR) = POR_sa(n_SNR) + 1;
+end
+
+if (abs(roots_toep(1) + theta1).^2/BWNN < 1) && (abs(roots_toep(2) + theta2).^2/BWNN < 1)
+ POR_toep(n_SNR) = POR_toep(n_SNR) + 1;
+end
+
+DOARMSE_idec(1,n_SNR) = DOARMSE_idec(1,n_SNR) + abs(roots_idec(1) + theta1).^2/BWNN/N_runs;
+DOARMSE_sa(1,n_SNR) = DOARMSE_sa(1,n_SNR) + abs(roots_sa(1) + theta1).^2/BWNN/N_runs;
+DOARMSE_toep(1,n_SNR) = DOARMSE_toep(1,n_SNR) + abs(roots_toep(1) + theta1).^2/BWNN/N_runs;
+
+DOARMSE_idec(2,n_SNR) = DOARMSE_idec(2,n_SNR) + abs(roots_idec(2) + theta2).^2/BWNN/N_runs;
+DOARMSE_sa(2,n_SNR) = DOARMSE_sa(2,n_SNR) + abs(roots_sa(2) + theta2).^2/BWNN/N_runs;
+DOARMSE_toep(2,n_SNR) = DOARMSE_toep(2,n_SNR) + abs(roots_toep(2) + theta2).^2/BWNN/N_runs;
+
+avgDOARMSE_idec(n_SNR) = mean(DOARMSE_idec(:,n_SNR));
+avgDOARMSE_sa(n_SNR) = mean(DOARMSE_sa(:,n_SNR));
+avgDOARMSE_toep(n_SNR) = mean(DOARMSE_toep(:,n_SNR));
+
+DOA_idec(1,n_SNR) = roots_idec(1);
+DOA_sa(1,n_SNR) = roots_sa(1);
+DOA_toep(1,n_SNR) = roots_toep(1);
+
+DOA_idec(2,n_SNR) = roots_idec(2);
+DOA_sa(2,n_SNR) = roots_sa(2);
+DOA_toep(2,n_SNR) = roots_toep(2);
+
+figure(f1);
+subplot(2,2,1);
+hold off;
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_idec(1,1:plotix_max))), 'r');
+hold on;
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_idec(2,1:plotix_max))), 'r:');
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_sa(1,1:plotix_max))), 'm');
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_sa(2,1:plotix_max))), 'm:');
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_toep(1,1:plotix_max))), 'b');
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_toep(2,1:plotix_max))), 'b:');
+title('Individual DOA RMSE');
+
+subplot(2,2,2);
+hold off;
+plot(SNR_vals(1:plotix_max), db(sqrt(avgDOARMSE_idec(1:plotix_max))), 'r');
+hold on;
+plot(SNR_vals(1:plotix_max), db(sqrt(avgDOARMSE_sa(1:plotix_max))), 'm');
+plot(SNR_vals(1:plotix_max), db(sqrt(avgDOARMSE_toep(1:plotix_max))), 'b');
+title('Average DOA RMSE');
+
+subplot(2,2,3);
+hold off;
+plot(SNR_vals(1:plotix_max), DOA_idec(1,1:plotix_max), 'r');
+hold on;
+plot(SNR_vals(1:plotix_max), DOA_idec(2,1:plotix_max), 'r:');
+plot(SNR_vals(1:plotix_max), DOA_sa(1,1:plotix_max), 'm');
+plot(SNR_vals(1:plotix_max), DOA_sa(2,1:plotix_max), 'm:');
+plot(SNR_vals(1:plotix_max), DOA_toep(1,1:plotix_max), 'b');
+plot(SNR_vals(1:plotix_max), DOA_toep(2,1:plotix_max), 'b:');
+
+title('DOA Estimates (deg)')
+
+subplot(2,2,4);
+hold off;
+plot(SNR_vals(1:plotix_max), POR_idec(1,1:plotix_max)/n_run, 'r-o');
+hold on;
+plot(SNR_vals(1:plotix_max), POR_sa(1,1:plotix_max)/n_run, 'm');
+plot(SNR_vals(1:plotix_max), POR_toep(1,1:plotix_max)/n_run, 'b');
+
+drawnow;
+
+end
+end
+%%
--- /dev/null
+clc;
+
+%SNR_vals = 20:3:60; N_SNR = length(SNR_vals);
+SNR_vals = -30:2:10; N_SNR = length(SNR_vals);
+
+N = 100;
+M = 10;
+
+BWNN = 2*asind(2/M);
+BWR = 2*asind(0.891/M);
+
+synth_noise = 0;
+use_symm = 1;
+
+DOARMSE_idec = zeros(2, N_SNR);
+DOARMSE_sa = zeros(2, N_SNR);
+DOARMSE_toep = zeros(2, N_SNR);
+
+avgDOARMSE_idec = zeros(1, N_SNR);
+avgDOARMSE_sa = zeros(1, N_SNR);
+avgDOARMSE_toep = zeros(1, N_SNR);
+
+POR_idec = zeros(1, N_SNR);
+POR_sa = zeros(1, N_SNR);
+POR_toep = zeros(1, N_SNR);
+
+DOA_idec = zeros(2, N_SNR);
+DOA_sa = zeros(2, N_SNR);
+DOA_toep = zeros(2, N_SNR);
+
+N_runs = 10;
+
+f1 = figure;
+
+for n_run=1:N_runs
+for n_SNR=1:N_SNR
+ if n_run > 1
+ plotix_max = N_SNR;
+ else
+ plotix_max = n_SNR;
+ end
+
+SNR1 = SNR_vals(n_SNR);
+
+disp(['SNR=', num2str(SNR1), 'dB']);
+
+SNR2 = SNR1;
+SNR3 = SNR2;
+theta1 = BWR;
+theta2 = 0;
+theta3 = -31.0;
+gamma = 45*0;
+d = 0.5;
+k = 2*pi;
+T = 0.0;
+omega1 = 2*pi;
+omega2 = 2*pi;
+omega3 = 2*pi;
+
+randn('seed',0);
+
+noise = randn(M,N) + j*randn(M,N);
+
+amp1 = sqrt(2)*10^(SNR1/20);
+amp2 = sqrt(2)*10^(SNR2/20);
+amp3 = 0*sqrt(2)*10^(SNR3/20);
+
+phi1 = -k*d*sin(theta1*pi/180);
+phi2 = -k*d*sin(theta2*pi/180);
+phi3 = -k*d*sin(theta3*pi/180);
+
+if use_symm
+ a1 = exp(j*phi1).^[-(M-1)/2:(M-1)/2]';
+ a2 = exp(j*phi2).^[-(M-1)/2:(M-1)/2]';
+ a3 = exp(j*phi3).^[-(M-1)/2:(M-1)/2]';
+else
+ a1 = exp(j*phi1).^[0:M-1]';
+ a2 = exp(j*phi2).^[0:M-1]';
+ a3 = exp(j*phi3).^[0:M-1]';
+end
+
+A = [a1 a2 a3];
+
+phase1 = j*2*pi*rand(N,1);
+phase2 = j*2*pi*rand(N,1);
+phase3 = j*2*pi*rand(N,1);
+
+n_corr = 1;
+if n_corr
+ %phase1 = zeros(N,1);
+ phase2 = phase1;
+ phase3 = phase2;
+end
+
+signal1 = amp1*exp(phase1);
+signal2 = amp2*exp(phase2);
+signal3 = amp3*exp(phase3);
+
+s = [signal1 signal2 signal3]';
+
+x = A*s + (1-synth_noise)*noise;
+
+delta = 0;
+
+N_mean = 100*M;
+L = M;
+K = M-L+1;
+
+L_sa = M-1;
+%L_sa=round(M/2);
+%L_sa=M;
+K_sa = M-L_sa+1;
+
+x_ = x;
+N_it = 20;
+
+int_f = @arrayest;
+
+R_scm = x*x'/N;
+R_toep = zeros(1,L);
+
+for lix=0:L-1
+ R_toep(lix+1) = mean(diag(R_scm,lix));
+end
+R_toep = toeplitz(R_toep);
+
+R_sa = zeros(L_sa);
+for kix=1:K_sa
+ R_sa = R_sa + x(kix:kix+L_sa-1,:)*x(kix:kix+L_sa-1,:)'/K_sa/N;
+end
+
+R = getDecorrelatedMatrix(x,L,N_it,int_f);
+
+[V,D] = eig(R);
+[V_toep,D_toep] = eig(R_toep);
+[V_scm,D_scm] = eig(R_scm);
+[V_sa,D_sa] = eig(R_sa);
+
+N_sig = 2;
+
+D = diag(D); [D_vals,D_ix]=sort(D); D_ix = D_ix(1:end-N_sig);
+D_toep = diag(D_toep); [D_toep_vals,D_toep_ix]=sort(D_toep);D_toep_ix = D_toep_ix(1:end-N_sig);
+D_sa = diag(D_sa); [D_sa_vals,D_sa_ix]=sort(D_sa); D_sa_ix = D_sa_ix(1:end-N_sig);
+D_scm = diag(D_scm); [D_scm_vals,D_scm_ix]=sort(D_scm);D_scm_ix = D_scm_ix(1:end-N_sig);
+
+V = V(:,D_ix);
+R_n = V*V';
+
+V_scm = V_scm(:,D_scm_ix);
+R_scm_n = V_scm*V_scm';
+
+V_toep = V_toep(:,D_toep_ix);
+R_toep_n = V_toep*V_toep';
+
+V_sa = V_sa(:,D_sa_ix);
+R_sa_n = V_sa*V_sa';
+
+% Calculate DOAs using root algorithms
+N_roots = 2*L-1;
+polcoeffs_idec = zeros(1,N_roots);
+polcoeffs_scm = zeros(1,N_roots);
+polcoeffs_toep = zeros(1,N_roots);
+polcoeffs_sa = zeros(1,N_roots);
+n_diag = -(L-1);
+
+for n_roots=1:N_roots
+ polcoeffs_idec(n_roots) = sum(diag(R_n,n_diag));
+ polcoeffs_scm(n_roots) = sum(diag(R_scm_n,n_diag));
+ polcoeffs_toep(n_roots) = sum(diag(R_toep_n,n_diag));
+ polcoeffs_sa(n_roots) = sum(diag(R_sa_n,n_diag));
+
+ n_diag = n_diag+1;
+end
+
+roots_idec = getMUSICangles(polcoeffs_idec, N_sig);
+roots_scm = getMUSICangles(polcoeffs_scm, N_sig);
+roots_toep = getMUSICangles(polcoeffs_toep, N_sig);
+roots_sa = getMUSICangles(polcoeffs_sa, N_sig);
+
+if (abs(roots_idec(1) + theta1).^2/BWNN < 1) && (abs(roots_idec(2) + theta2).^2/BWNN < 1)
+ POR_idec = POR_idec
+end
+
+DOARMSE_idec(1,n_SNR) = DOARMSE_idec(1,n_SNR) + abs(roots_idec(1) + theta1).^2/BWNN/N_runs;
+DOARMSE_sa(1,n_SNR) = DOARMSE_sa(1,n_SNR) + abs(roots_sa(1) + theta1).^2/BWNN/N_runs;
+DOARMSE_toep(1,n_SNR) = DOARMSE_toep(1,n_SNR) + abs(roots_toep(1) + theta1).^2/BWNN/N_runs;
+
+DOARMSE_idec(2,n_SNR) = DOARMSE_idec(2,n_SNR) + abs(roots_idec(2) + theta2).^2/BWNN/N_runs;
+DOARMSE_sa(2,n_SNR) = DOARMSE_sa(2,n_SNR) + abs(roots_sa(2) + theta2).^2/BWNN/N_runs;
+DOARMSE_toep(2,n_SNR) = DOARMSE_toep(2,n_SNR) + abs(roots_toep(2) + theta2).^2/BWNN/N_runs;
+
+avgDOARMSE_idec(n_SNR) = mean(DOARMSE_idec(:,n_SNR));
+avgDOARMSE_sa(n_SNR) = mean(DOARMSE_sa(:,n_SNR));
+avgDOARMSE_toep(n_SNR) = mean(DOARMSE_toep(:,n_SNR));
+
+DOA_idec(1,n_SNR) = roots_idec(1);
+DOA_sa(1,n_SNR) = roots_sa(1);
+DOA_toep(1,n_SNR) = roots_toep(1);
+
+DOA_idec(2,n_SNR) = roots_idec(2);
+DOA_sa(2,n_SNR) = roots_sa(2);
+DOA_toep(2,n_SNR) = roots_toep(2);
+
+figure(f1);
+subplot(2,2,1);
+hold off;
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_idec(1,1:plotix_max))), 'r');
+hold on;
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_idec(2,1:plotix_max))), 'r:');
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_sa(1,1:plotix_max))), 'm');
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_sa(2,1:plotix_max))), 'm:');
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_toep(1,1:plotix_max))), 'b');
+plot(SNR_vals(1:plotix_max), db(sqrt(DOARMSE_toep(2,1:plotix_max))), 'b:');
+title('Individual DOA RMSE');
+
+subplot(2,2,2);
+hold off;
+plot(SNR_vals(1:plotix_max), db(sqrt(avgDOARMSE_idec(1:plotix_max))), 'r');
+hold on;
+plot(SNR_vals(1:plotix_max), db(sqrt(avgDOARMSE_sa(1:plotix_max))), 'm');
+plot(SNR_vals(1:plotix_max), db(sqrt(avgDOARMSE_toep(1:plotix_max))), 'b');
+title('Average DOA RMSE');
+
+subplot(2,2,3);
+hold off;
+plot(SNR_vals(1:plotix_max), DOA_idec(1,1:plotix_max), 'r');
+hold on;
+plot(SNR_vals(1:plotix_max), DOA_idec(2,1:plotix_max), 'r:');
+plot(SNR_vals(1:plotix_max), DOA_sa(1,1:plotix_max), 'm');
+plot(SNR_vals(1:plotix_max), DOA_sa(2,1:plotix_max), 'm:');
+plot(SNR_vals(1:plotix_max), DOA_toep(1,1:plotix_max), 'b');
+plot(SNR_vals(1:plotix_max), DOA_toep(2,1:plotix_max), 'b:');
+
+title('DOA Estimates (deg)')
+
+subplot(2,2,4);
+
+drawnow;
+
+end
+end
+%%
\@writefile{toc}{\contentsline {section}{\numberline {I}Introduction}{1}}
\@writefile{toc}{\contentsline {section}{\numberline {II}Background}{1}}
\@writefile{toc}{\contentsline {subsection}{\numberline {\unhbox \voidb@x \hbox {II-A}}Array Model, our Objective [the Spatial Covariance Matrix, and its Estimation]}{1}}
-\newlabel{eq:def_toeplitzcm}{{6}{1}}
\citation{vantrees:OAP}
\citation{schmidt:music}
\citation{roy:esprit}
\citation{shan:ss}
\citation{45613}
\citation{124962}
+\newlabel{eq:def_toeplitzcm}{{6}{2}}
\newlabel{eq:def_scm}{{7}{2}}
-\@writefile{toc}{\contentsline {subsection}{\numberline {\unhbox \voidb@x \hbox {II-B}}DOA Estimation: Root-MUSIC (?)}{2}}
-\@writefile{toc}{\contentsline {subsection}{\numberline {\unhbox \voidb@x \hbox {II-C}}Problems/Deviations with the Covariance Matrix}{2}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {\unhbox \voidb@x \hbox {II-B}}DOA Estimation with Root-MUSIC}{2}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {\unhbox \voidb@x \hbox {II-C}}Robust Covariance Matrix Estimation}{2}}
\newlabel{eq:def_Rra}{{12}{2}}
-\@writefile{toc}{\contentsline {section}{\numberline {III}(Iterative?) Covariance Matrix Decorrelation}{2}}
+\@writefile{toc}{\contentsline {section}{\numberline {III}Covariance Matrix Decorrelation Through Optimal Toeplitz Completion}{2}}
\newlabel{eq:defR_s}{{13}{2}}
\newlabel{eq:def_beamforming}{{14}{2}}
\newlabel{eq:def_beamforming_wavenumber}{{15}{2}}
\newlabel{eq:def_mvdr}{{16}{2}}
-\newlabel{eq:def_Sest}{{18}{2}}
\bibstyle{IEEEtran}
\bibdata{IterativeCovarianceMatrixDecorrelation}
\bibcite{1144039}{1}
\bibcite{schmidt:music}{5}
\bibcite{roy:esprit}{6}
\bibcite{shan:ss}{7}
+\newlabel{eq:def_Sest}{{18}{3}}
\newlabel{eq:def_Re}{{20}{3}}
\newlabel{eq:def_S0_const}{{21}{3}}
\newlabel{eq:def_S0}{{22}{3}}
\newlabel{eq:def_S1}{{23}{3}}
-\@writefile{toc}{\contentsline {subsection}{\numberline {\unhbox \voidb@x \hbox {III-A}}Computational Complexity}{3}}
\@writefile{toc}{\contentsline {section}{\numberline {IV}Simulations and Discussion}{3}}
\@writefile{toc}{\contentsline {section}{\numberline {V}Conclusion}{3}}
\@writefile{toc}{\contentsline {section}{References}{3}}
-This is pdfTeX, Version 3.1415926-2.3-1.40.12 (TeX Live 2011) (format=pdflatex 2011.7.3) 4 JUN 2012 11:31
+This is pdfTeX, Version 3.1415926-2.3-1.40.12 (TeX Live 2011) (format=pdflatex 2011.7.3) 4 JUN 2012 13:49
entering extended mode
restricted \write18 enabled.
%&-line parsing enabled.
LaTeX Font Info: ... okay on input line 424.
LaTeX Font Info: Checking defaults for U/cmr/m/n on input line 424.
LaTeX Font Info: ... okay on input line 424.
-LaTeX Font Info: Try loading font information for U+msa on input line 540.
+LaTeX Font Info: Try loading font information for U+msa on input line 541.
(/usr/local/texlive/2011/texmf-dist/tex/latex/amsfonts/umsa.fd
File: umsa.fd 2009/06/22 v3.00 AMS symbols A
)
-LaTeX Font Info: Try loading font information for U+msb on input line 540.
+LaTeX Font Info: Try loading font information for U+msb on input line 541.
(/usr/local/texlive/2011/texmf-dist/tex/latex/amsfonts/umsb.fd
File: umsb.fd 2009/06/22 v3.00 AMS symbols B
)
-Underfull \hbox (badness 1472) in paragraph at lines 528--540
+Underfull \hbox (badness 1472) in paragraph at lines 528--541
\OT1/ptm/m/n/10 and show that a slightly mod-i-fied form of re-dun-dancy
[]
-Underfull \hbox (badness 10000) in paragraph at lines 528--540
+Underfull \hbox (badness 10000) in paragraph at lines 528--541
[]
-Underfull \hbox (badness 1895) in paragraph at lines 540--543
+Underfull \hbox (badness 1895) in paragraph at lines 541--544
\OT1/ptm/m/n/10 iza-tion'' (un-der-ly-ing causes not im-me-di-ately clear)?
[]
-Underfull \hbox (badness 10000) in paragraph at lines 576--582
+Underfull \hbox (badness 10000) in paragraph at lines 577--583
\OT1/ptm/m/n/10 When the sig-nals of in-ter-est are un-cor-re-lated, i.e.
[]
]
-[2] (./IterativeCovarianceMatrixDecorrelation.bbl) [3]
+[2]
+Overfull \hbox (7.58418pt too wide) detected at line 733
+[][] \OT1/cmr/m/n/10 = [] [][] []\OML/cmm/m/it/10 :
+ []
+
+
+Overfull \hbox (18.96274pt too wide) in paragraph at lines 812--812
+[]
+ []
+
+(./IterativeCovarianceMatrixDecorrelation.bbl) [3] [4
+
+]
(./IterativeCovarianceMatrixDecorrelation.aux) )
Here is how much of TeX's memory you used:
1832 strings out of 493633
26420 string characters out of 3143378
- 87895 words of memory out of 3000000
+ 88895 words of memory out of 3000000
5158 multiletter control sequences out of 15000+200000
40090 words of font info for 77 fonts, out of 3000000 for 9000
831 hyphenation exceptions out of 8191
- 27i,11n,24p,279b,263s stack positions out of 5000i,500n,10000p,200000b,50000s
+ 27i,13n,24p,293b,263s stack positions out of 5000i,500n,10000p,200000b,50000s
{/usr/local/texlive/2011/texmf-
dist/fonts/enc/dvips/base/8r.enc}</usr/local/texlive/2011/texmf-dist/fonts/type
1/public/amsfonts/cm/cmbx10.pfb></usr/local/texlive/2011/texmf-dist/fonts/type1
ts/cm/cmsy10.pfb></usr/local/texlive/2011/texmf-dist/fonts/type1/public/amsfont
s/cm/cmsy5.pfb></usr/local/texlive/2011/texmf-dist/fonts/type1/public/amsfonts/
cm/cmsy7.pfb></usr/local/texlive/2011/texmf-dist/fonts/type1/public/amsfonts/sy
-mbols/msbm10.pfb></usr/local/texlive/2011/texmf-dist/fonts/type1/urw/times/utmb
-8a.pfb></usr/local/texlive/2011/texmf-dist/fonts/type1/urw/times/utmbi8a.pfb></
-usr/local/texlive/2011/texmf-dist/fonts/type1/urw/times/utmr8a.pfb></usr/local/
-texlive/2011/texmf-dist/fonts/type1/urw/times/utmri8a.pfb>
-Output written on IterativeCovarianceMatrixDecorrelation.pdf (3 pages, 178376 b
+mbols/msam10.pfb></usr/local/texlive/2011/texmf-dist/fonts/type1/public/amsfont
+s/symbols/msbm10.pfb></usr/local/texlive/2011/texmf-dist/fonts/type1/urw/times/
+utmb8a.pfb></usr/local/texlive/2011/texmf-dist/fonts/type1/urw/times/utmbi8a.pf
+b></usr/local/texlive/2011/texmf-dist/fonts/type1/urw/times/utmr8a.pfb></usr/lo
+cal/texlive/2011/texmf-dist/fonts/type1/urw/times/utmri8a.pfb>
+Output written on IterativeCovarianceMatrixDecorrelation.pdf (4 pages, 189321 b
ytes).
PDF statistics:
- 79 PDF objects out of 1000 (max. 8388607)
- 57 compressed objects within 1 object stream
+ 86 PDF objects out of 1000 (max. 8388607)
+ 62 compressed objects within 1 object stream
0 named destinations out of 1000 (max. 500000)
1 words of extra memory for PDF output out of 10000 (max. 10000000)
reduced for each iteration. The new technique is shown through\r
simulations to have better treshold and resolution properties than\r
both redundancy averaging avd spatial averaging, as well as a DOA RMSE\r
-that is superior to redundancy averaging and approaches spatial averaging.\\\r
+that is superior to redundancy averaging and approaches spatial\r
+averaging as SNR increases.\\\r
\\\r
{\bf Notes:}\r
\begin{itemize}\r
\hat{\mat{R}} = \frac{1}{N}\sum_{n=0}^{N-1}{\vec{x}[n]\vec{x}^{H}[n]}.\r
\end{equation}\r
\r
-\subsection{DOA Estimation: Root-MUSIC (?)}\r
-\r
-There are several\r
-non-parametric and parametric techniques for achieving our primary\r
+\subsection{DOA Estimation with Root-MUSIC}\r
+There are several non-parametric and parametric techniques for achieving our primary\r
goal. The non-parametric methods are typically beamscan algorithms\r
(see e.g. \cite{vantrees:OAP}), which form a spatial spectrum over all wavenumbers and locate the $D$\r
largest peaks within. The most famous parametric methods are MUSIC\r
The wavenumbers are the roots of the MUSIC\r
power spectrum:\r
\begin{equation}\r
-Q_{MUSIC}(k) = \vec{v}^{H}(k)\outerp{\mat{U}_{N}}\vec{v}(k)\r
+Q_{MUSIC}(k) = \vec{v}^{H}(k)\outerp{\mat{V}_{N}}\vec{v}(k)\r
\end{equation}\r
\r
-\subsection{Problems/Deviations with the Covariance Matrix}\r
+\subsection{Robust Covariance Matrix Estimation}\r
In practice, the sample covariance matrix of (\ref{eq:def_scm}) will\r
not have a Toeplitz structure. This can be due to different\r
reasons. If $\hat{\mat{R}}$ is based on a finite number of temporal samples, the\r
(\ref{eq:def_toeplitzcm}) can be written in the wavenumber domain:\r
\begin{equation}\r
\label{eq:defR_s}\r
-\left[ \mat{R}_{s} \right]_{m,n} = \int_{-\pi}^{\pi}{\abssq{S(k)}e^{i(m-n)k}\,dk}\r
+\left[ \mat{R}_{s} \right]_{m,n} = \frac{1}{2\pi}\int_{-\pi}^{\pi}{\abssq{S(k)}e^{i(m-n)k}\,dk}\r
\end{equation}\r
The wavefield magnitude at wavenumber $k_{s}$, $\abssq{S(k_{s})}$, is generally unknown, but can be\r
estimated using beamforming:\r
domain using a Fourier transformation:\r
\begin{equation}\r
\label{eq:def_beamforming_wavenumber}\r
-\vec{w}^{H}\mat{R}\vec{w} = \E{\abssq{\int_{-\pi}^{\pi}{W(k)S(k)\,dk}}}.\r
+\vec{w}^{H}\mat{R}\vec{w} = \E{\abssq{\frac{1}{2\pi}\int_{-\pi}^{\pi}{W(k)S(k)\,dk}}}.\r
\end{equation}\r
An optimum estimator of $\abssq{S(k)}$ is the Minimum Variance Distortionless\r
Response (MVDR) beamformer, which is defined as:\r
Under the assumption of an uncorrelated wavefield, (\ref{eq:def_mvdr})\r
can be written in the wavenumber domain as in (\ref{eq:def_beamforming_wavenumber}):\r
\begin{align}\r
-\r
& \vec{w}_{MV} = \argmin_{\vec{w}}\r
-\int_{-\pi}^{\pi}{\abssq{W(k)}\E{\abssq{S(k)}}\,dk} \\\r
+\frac{1}{2\pi}\int_{-\pi}^{\pi}{\abssq{W(k)}\E{\abssq{S(k)}}\,dk} \\\r
& \mbox{subject to } W(k_{s})=1 \mbox{ and } w_{m} = 0 \mbox{ for }\r
m<0 \mbox{ or } m\r
\geq M \nonumber\r
\label{eq:def_Sest}\r
& \abssq{\hat{S}^{(e+1)}(k_{s})} =\r
\vec{w}_{MV}^{(e+1)H}\hat{\mat{R}}\vec{w}_{MV}^{(e+1)} \\\r
-& \mbox{where } \vec{w}_{MV}^{(e+1)} = \argmin_{\vec{w}} \int_{-\pi}^{\pi}{\abssq{W(k)\hat{S}^{(e)}(k)}\,dk}. \nonumber\r
+& \mbox{where } \vec{w}_{MV}^{(e+1)} = \argmin_{\vec{w}} \frac{1}{2\pi}\int_{-\pi}^{\pi}{\abssq{W(k)\hat{S}^{(e)}(k)}\,dk}. \nonumber\r
\end{align}\r
It is well known that the solution of this MVDR optimization problem\r
is:\r
\vec{w}_{MV}^{(e+1)} = \frac{\r
\hat{\mat{R}}^{-(e)}\vec{v}(k)}{\vec{v}^{H}(k)\hat{\mat{R}}^{-(e)}\vec{v}(k)},\r
\end{equation}\r
-where $\hat{\mat{R}}^{-(e)}$ is shorthand for $\left(\hat{\mat{R}}^{(e)}\right)^{-1}$.\r
-\r
-This estimate can be used with (\ref{eq:defR_s}) to find the\r
-Toeplitz covariance matrix estimate of the uncorrelated wavefield:\r
+where $\hat{\mat{R}}^{-(e)}$ is shorthand for\r
+$\left(\hat{\mat{R}}^{(e)}\right)^{-1}$.This estimate can be used with\r
+(\ref{eq:defR_s}) to specify an iterative Toeplitz covariance matrix\r
+estimator of the associated uncorrelated wavefield:\r
\begin{equation}\r
\label{eq:def_Re}\r
-\left[ \hat{\mat{R}}^{(e+1)} \right]_{m,n} = \int_{-\pi}^{\pi} { \frac{\vec{v}^{H}(k)\hat{\mat{R}}^{-(e)}\hat{\mat{R}}\hat{\mat{R}}^{-(e)}\vec{v}(k)}{\abssq{\vec{v}^{H}(k)\hat{\mat{R}}^{-(e)}\vec{v}(k)}} e^{i(m-n)k}\,dk}.\r
+\left[ \hat{\mat{R}}^{(e+1)} \right]_{m,n} = \frac{1}{2\pi}\int_{-\pi}^{\pi} { \frac{\vec{v}^{H}(k)\hat{\mat{R}}^{-(e)}\hat{\mat{R}}\hat{\mat{R}}^{-(e)}\vec{v}(k)}{\abssq{\vec{v}^{H}(k)\hat{\mat{R}}^{-(e)}\vec{v}(k)}} e^{i(m-n)k}\,dk}.\r
\end{equation}\r
As for any other iterative algorithm, there is the question of\r
intialization. A plausible initial estimate $\abssq{\hat{S}^{(0)}(k)}$\r
sense to initialize with the conventional beamformer directly:\r
\begin{equation}\r
\label{eq:def_S0}\r
-\abssq{\hat{S}^{(0)}(k)} = \vec{v}^{H}(k)\hat{\mat{R}}\vec{v}(k)\r
+\abssq{\hat{S}^{(0)}(k)} = \frac{1}{M^{2}}\vec{v}^{H}(k)\hat{\mat{R}}\vec{v}(k)\r
\end{equation}\r
An interesting observation comes when inserting (\ref{eq:def_S0}) into\r
(\ref{eq:def_Sest}) and, subsequently, the result into\r
(\ref{eq:def_Re}) to estimate $\hat{\mat{R}}^{(1)}$, namely that the\r
-first covariance matrix estimate becomes:\r
+first covariance matrix estimate becomes (see Appendix A for a proof):\r
\begin{equation}\r
\label{eq:def_S1}\r
[\hat{\mat{R}}^{(1)}]_{m,n} =\r
-\frac{1}{M}\sum_{p=0}^{P-1}{\left[\hat{\mat{R}}\right]_{p,P+p}}, P = \abs{m-n}.\r
+\frac{1}{M^{2}}\sum_{p=0}^{P-1}{\left[\hat{\mat{R}}\right]_{p,P+p}}, P = \abs{m-n}.\r
\end{equation}\r
This is equal to the redundancy averaged estimator of (\ref{eq:def_Rra}) except for a term\r
$-P$ in the denominator of the leading fraction. As we will\r
ultimately makes more sense to start with (\ref{eq:def_S1}) instead of\r
the conventional beamformer estimate.\r
\r
-{\bf Other points:}\r
-\begin{itemize}\r
-\item Computational complexity.\r
-\item Convergence.\r
-\end{itemize}\r
+%{\bf Other points:}\r
+%\begin{itemize}\r
+%\item Computational complexity.\r
+%\item Convergence.\r
+%\end{itemize}\r
+\r
+%\subsection{Computational Complexity}\r
+\r
+Unfortunately, it is not trivial to prove the convergence of the\r
+iterative estimator in (\ref{eq:def_Re}). In the next section, we will\r
+demonstrate convergence for a range of cases through simulations.\r
\r
-\subsection{Computational Complexity}\r
It is worth noting that the suggested\r
-method is much more computationally demanding than both spatial and\r
+method is more computationally demanding than both spatial and\r
redundancy averaging. Each iteration requires the formation of an $M\r
\times M$ Toeplitz matrix, where the value along each of the $M$\r
diagonals is calculated from the integral in (\ref{eq:def_Re}). As\r
solution, while the subsequent solutions result in improved DOA\r
estimates and resolution. \r
\r
+\section*{Appendix A: Proof of Equation (\ref{eq:def_S1})}\r
+Because the covariance matrix of (\ref{eq:defR_s}) is Toeplitz, we\r
+will investigate the value along the $l^{th}$ diagonal:\r
+\begin{align}\r
+\left[ \hat{\mat{R}}^{(1)} \right]_{m,n} & =\r
+\frac{1}{2\pi}\int_{-\pi}^{\pi}{\abssq{\hat{S}^{(0)}(k)}e^{ilk}\,dk}, \nonumber \\ &\r
+\mbox{for } l\triangleq m-n \nonumber \\\r
+& = \frac{1}{2\pi M^{2}}\int_{-\pi}^{\pi}{\vec{v}^{H}(k)\hat{\mat{R}}\vec{v}(k)e^{ilk}\,dk} \nonumber \\\r
+& =\r
+\frac{1}{2\pi TM^{2}}\sum_{n=0}^{T-1}{\sum_{m'=0}^{M-1}{\sum_{n'=0}^{M-1}{x_{m'}[t]x_{n'}^{*}[t]\int_{-\pi}^{\pi}{e^{i(l-l')k}\,dk}}}}\r
+,\nonumber \\ & \mbox{for } l' \triangleq m'-n'\nonumber \\\r
+& = \frac{1}{M^{2}}\sum_{m'=0}^{M-1}{\sum_{n'=0}^{M-1}{\delta_{l-l'}\frac{1}{T}\sum_{t=0}^{T-1}{x_{m'}[t]x_{n'}^{*}[t]}}} \nonumber \\\r
+& = \frac{1}{M^{2}} \sum_{p=0}^{l-1}{\left[ \hat{\mat{R}} \right]_{p,p+l}}, \nonumber \\\r
+\end{align}\r
+which we confirm is equal to (\ref{eq:def_S1}) and nearly equal to the redundancy averaged estimate of the\r
+covariance matrix as given in (\ref{eq:def_Rra}), except for the denominator in the leading fraction.\r
+\r
\bibliographystyle{IEEEtran}\r
\bibliography{IterativeCovarianceMatrixDecorrelation}\r
\r