Description of SpectralMajorisationDemo

0001 function SpectralMajorisationDemo();
0002 %SpectralMajorisationDemo()
0003 %
0004 %  This demo demonstrated spectral majorisation achieved by a number of
0005 %  iterative PEVD algorithms after 50 and 200 iterations for a situation
0006 %  with known ground truth.
0007 %  The ground truth is provided by a source model generated by the
0008 %  function GenerateMIMOSources(), which here contains 8 sources with
0009 %  spectrally majorised PSDs and a paraunitary convolutive mixing matrix
0010 %  of order 16.
0011 %
0012 %  The results reflect approximately those shown in Figures 4 and 5 of [1].
0013 %
0014 %  Reference:
0015 %
0016 %  [1] S. Redif, S. Weiss, J.G. McWhirter: "Sequential Matrix Diagonalisation
0017 %      Algorithms for Polynomial EVD of Parahermitian Matrices", to appear in
0018 %      IEEE Transactions on Signal Processing.
0019 
0020 % S. Weiss, University of Strathclyde, 2/11/2014
0021 
0022 %-----------------------------------------------------
0023 %  Ground truth and simulation parameters
0024 %-----------------------------------------------------
0025 % set seeds for a specific source model
0026 randn('seed',10); rand('seed',10);
0027 % generate source model with the following parameters:
0028 L = 8;         %  # of sources
0029 P = 16;        %  order of source innovation filter
0030 M = 8;         %  # of sensors
0031 K = 16;        %  order of paraunitary mixing matrix
0032 gamma = 0.08;  %  max radii of zeros
0033 % convolutive mixing model
0034 [H,D,F] = GenerateMIMOSources(L,P,M,K,gamma);
0035 % PSDs of ground truth sources
0036 P = PolyMatDiagSpec(D,1024);    
0037 % space-time covariance matrix at sensors
0038 R = PolyMatConv(H,PolyMatConv(D,ParaHerm(H)));
0039 % some display parameters
0040 LWidth = 1;    % line widths for spectra
0041 LWidth2 = 3;   % line width for ground truth
0042 ShadeColour=[.7 .7 .7];  % grey scale for ground truth
0043 w= (0:1023)/512;   % frequency scale
0044 
0045 %-----------------------------------------------------
0046 %   SMD, ME-SMD, SBR2 and SBR2C after 50 iterations
0047 %-----------------------------------------------------
0048 figure(1);
0049 clf;
0050 [H_est,D_est] = SMD(R,50,0.00000001,0.00000001,'SMD');
0051 P_est = PolyMatDiagSpec(D_est,1024);
0052 subplot(221);
0053 % entries for legend
0054 plot([-2 -1],[1 2],'b-o','LineWidth',LWidth);
0055 hold on;
0056 plot([-2 -1],[1 2],'r-.*','LineWidth',LWidth);
0057 plot([-2 -1],[1 2],'k--s','LineWidth',LWidth);
0058 plot([-2 -1],[1 2],'-+','Color',[0 .5 0],'LineWidth',LWidth);
0059 AlgorithmPSDs50Plot(P_est,P,LWidth,LWidth2,ShadeColour);
0060 % legend and labels
0061 dummy = legend('$l=1$','$l=2$','$l=3$','$l=4$');
0062 set(dummy,'interpreter','latex','location','SouthWest');
0063 text(0.05,1,'(a)');
0064 title('SMD after 50 iterations');
0065 
0066 % ME-SMD
0067 subplot(222);
0068 [H_est,D_est] = SMD(R,50,0.00000001,0.00000001,'MESMD');
0069 P_est = PolyMatDiagSpec(D_est,1024);
0070 AlgorithmPSDs50Plot(P_est,P,LWidth,LWidth2,ShadeColour);
0071 text(0.05,1,'(b)');
0072 title('ME-SMD after 50 iterations');
0073 
0074 % SBR2
0075 subplot(223);
0076 [H_est,D_est] = SBR2(R,50,0.00000001,0.00000001,'SBR2');
0077 P_est = PolyMatDiagSpec(D_est,1024);
0078 AlgorithmPSDs50Plot(P_est,P,LWidth,LWidth2,ShadeColour);
0079 text(0.05,1,'(c)');
0080 title('SBR2 after 50 iterations');
0081 
0082 % SBR2C
0083 subplot(224);
0084 [H_est,D_est] = SBR2(R,50,0.00000001,0.00000001,'SBR2C');
0085 P_est = PolyMatDiagSpec(D_est,1024);
0086 AlgorithmPSDs50Plot(P_est,P,LWidth,LWidth2,ShadeColour);
0087 text(0.05,1,'(d)');
0088 title('SBR2C after 50 iterations');
0089 
0090 %-----------------------------------------------------
0091 %   SMD, ME-SMD, SBR2 and SBR2C after 200 iterations
0092 %-----------------------------------------------------
0093 figure(2);
0094 clf;
0095 [H_est,D_est] = SMD(R,200,0.00000001,0.00000001,'SMD');
0096 P_est = PolyMatDiagSpec(D_est,1024);
0097 subplot(221);
0098 % entries for legend
0099 plot([-2 -1],[1 2],'b-o','LineWidth',LWidth);
0100 hold on;
0101 plot([-2 -1],[1 2],'r-.*','LineWidth',LWidth);
0102 plot([-2 -1],[1 2],'k--s','LineWidth',LWidth);
0103 plot([-2 -1],[1 2],'-+','Color',[0 .5 0],'LineWidth',LWidth);
0104 plot([-2 -1],[1 2],'b-.d','LineWidth',LWidth);
0105 plot([-2 -1],[1 2],'r--x','LineWidth',LWidth);
0106 plot([-2 -1],[1 2],'k-o','LineWidth',LWidth);
0107 plot([-2 -1],[1 2],'-.*','Color',[0 .5 0],'LineWidth',LWidth);
0108 AlgorithmPSDs200Plot(P_est,P,LWidth,LWidth2,ShadeColour);
0109 dummy = legend('$l=1$','$l=2$','$l=3$','$l=4$','$l=5$','$l=6$','$l=7$','$l=8$');
0110 set(dummy,'interpreter','latex','location','SouthWest');
0111 text(0.05,1,'(a)');
0112 title('SMD after 200 iterations');
0113 
0114 % ME-SMD
0115 subplot(222);
0116 [H_est,D_est] = SMD(R,200,0.00000001,0.00000001,'MESMD');
0117 P_est = PolyMatDiagSpec(D_est,1024);
0118 AlgorithmPSDs200Plot(P_est,P,LWidth,LWidth2,ShadeColour);
0119 text(0.05,1,'(b)');
0120 title('ME-SMD after 200 iterations');
0121 
0122 % SBR2
0123 subplot(223);
0124 [H_est,D_est] = SBR2(R,200,0.00000001,0.00000001,'SBR2');
0125 P_est = PolyMatDiagSpec(D_est,1024);
0126 AlgorithmPSDs200Plot(P_est,P,LWidth,LWidth2,ShadeColour);
0127 text(0.05,1,'(c)');
0128 title('SBR2 after 200 iterations');
0129 
0130 % SBR2C
0131 subplot(224);
0132 [H_est,D_est] = SBR2(R,200,0.00000001,0.00000001,'SBR2C');
0133 P_est = PolyMatDiagSpec(D_est,1024);
0134 AlgorithmPSDs200Plot(P_est,P,LWidth,LWidth2,ShadeColour);
0135 text(0.05,1,'(d)');
0136 title('SBR2C after 200 iterations');
0137 
0138 
0139 %=============================================
0140 function AlgorithmPSDs50Plot(P_est,P,LWidth,LWidth2,ShadeColour);
0141 % Plot power spectral densities of first four columns in P and P_est
0142 % background
0143 w = (0:1023)/512;
0144 plot(w,10*log10(abs(P(:,1))),'Color',ShadeColour,'LineStyle','-','LineWidth',LWidth2);
0145 hold on;
0146 plot(w,10*log10(abs(P(:,2))),'Color',ShadeColour,'LineStyle','-.','LineWidth',LWidth2);
0147 plot(w,10*log10(abs(P(:,3))),'Color',ShadeColour,'LineStyle','--','LineWidth',LWidth2);
0148 plot(w,10*log10(abs(P(:,4))),'Color',ShadeColour,'LineStyle','-','LineWidth',LWidth2);
0149 % lines
0150 plot(w,10*log10(abs(P_est(:,1))),'b-','LineWidth',LWidth);
0151 plot(w,10*log10(abs(P_est(:,2))),'r-.','LineWidth',LWidth);
0152 plot(w,10*log10(abs(P_est(:,3))),'k--','LineWidth',LWidth);
0153 plot(w,10*log10(abs(P_est(:,4))),'-','Color',[0 .5 0],'LineWidth',LWidth);
0154 % markers
0155 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,1))),'bo','LineWidth',LWidth);
0156 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,2))),'r*','LineWidth',LWidth);
0157 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,3))),'ks','LineWidth',LWidth);
0158 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,4))),'+','Color',[0 .5 0],'LineWidth',LWidth);
0159 % legend and labels
0160 xlabel('norm. angular frequency $\Omega/\pi$','fontsize',14,...
0161   'Interpreter','Latex');
0162 ylabel('$10\log_{10}S^{(50)}_{s,l}(e^{j\Omega})$ / [dB]','fontsize',...
0163   14,'Interpreter','Latex');
0164 axis([0 2 -8 2]);
0165 grid on;
0166 
0167 %=============================================
0168 function AlgorithmPSDs200Plot(P_est,P,LWidth,LWidth2,ShadeColour);
0169 % Plot power spectral densities of all eight columns in P and P_est
0170 % background
0171 w= (0:1023)/512;
0172 
0173 % background
0174 plot(w,10*log10(abs(P(:,1))),'Color',ShadeColour,'LineStyle','-','LineWidth',LWidth2);
0175 hold on;
0176 plot(w,10*log10(abs(P(:,2))),'Color',ShadeColour,'LineStyle','-.','LineWidth',LWidth2);
0177 plot(w,10*log10(abs(P(:,3))),'Color',ShadeColour,'LineStyle','--','LineWidth',LWidth2);
0178 plot(w,10*log10(abs(P(:,4))),'Color',ShadeColour,'LineStyle','-','LineWidth',LWidth2);
0179 plot(w,10*log10(abs(P(:,5))),'Color',ShadeColour,'LineStyle','-.','LineWidth',LWidth2);
0180 plot(w,10*log10(abs(P(:,6))),'Color',ShadeColour,'LineStyle','--','LineWidth',LWidth2);
0181 plot(w,10*log10(abs(P(:,7))),'Color',ShadeColour,'LineStyle','-','LineWidth',LWidth2);
0182 plot(w,10*log10(abs(P(:,8))),'Color',ShadeColour,'LineStyle','-.','LineWidth',LWidth2);
0183 % lines
0184 plot(w,10*log10(abs(P_est(:,1))),'b-','LineWidth',LWidth);
0185 plot(w,10*log10(abs(P_est(:,2))),'r-.','LineWidth',LWidth);
0186 plot(w,10*log10(abs(P_est(:,3))),'k--','LineWidth',LWidth);
0187 plot(w,10*log10(abs(P_est(:,4))),'-','Color',[0 .5 0],'LineWidth',LWidth);
0188 plot(w,10*log10(abs(P_est(:,5))),'b-.','LineWidth',LWidth);
0189 plot(w,10*log10(abs(P_est(:,6))),'r--','LineWidth',LWidth);
0190 plot(w,10*log10(abs(P_est(:,7))),'k-','LineWidth',LWidth);
0191 plot(w,10*log10(abs(P_est(:,8))),'-.','Color',[0 .5 0],'LineWidth',LWidth);
0192 % markers
0193 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,1))),'bo','LineWidth',LWidth);
0194 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,2))),'r*','LineWidth',LWidth);
0195 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,3))),'ks','LineWidth',LWidth);
0196 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,4))),'+','Color',[0 .5 0],'LineWidth',LWidth);
0197 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,5))),'bd','LineWidth',LWidth);
0198 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,6))),'rx','LineWidth',LWidth);
0199 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,7))),'ko','LineWidth',LWidth);
0200 plot(w(64:128:end),10*log10(abs(P_est(64:128:end,8))),'g*','Color',[0 .5 0],'LineWidth',LWidth);
0201 % axes
0202 xlabel('norm. angular frequency $\Omega/\pi$','fontsize',14,...
0203   'Interpreter','Latex');
0204 ylabel('$10\log_{10}S^{(200)}_{s,l}(e^{j\Omega})$ / [dB]','fontsize',...
0205   14,'Interpreter','Latex');
0206 axis([0 2 -14 2]);
0207 grid on;
0208
SpectralMajorisationDemo

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
