SMDDemo(SourceSwitch);
SMDDemo(SourceSwitch) demonstrates the polynomial EVD approximated by
the function SMD() for different type of parahermitian matrices. The
input parameter SourceSwitch can be selected as:
SMDDemo('Random') produces an abritrary complex valued 5x5 polynomial
parahermitian matrix of order 10.
SMDDemo('Simple1') decomposes a 3x3 real valued polynomial parahermitian
matrix R(z) of order 4:
| 1 -.4*z 0 |
R(z) = | -.4*z 1 .5*z^{-2} |
| 0 .5*z^2 3 |
SMDDemo('Simple2') decomposes a 5x5 complex valued polynomial
parahermitian matrix R(z) of order 4:
| 1 .5*z^2 -.4j*z^{-1} 0 .2j*z^2 |
| .5*z^{-2} 1 0 0 0 |
R(z) = | .4j*z 0 3 .3z^{-1} 0 |
| 0 0 .3*z .5 0 |
|-.2j*z^{-2} 0 0 0 .25 |
This matrix admits a PEVD with a non-polynomial Gamma(z).
The function generates a number of outputs to highlight diagonalisation
and spectral majorisation.
Input parameter:
SourceSwitch selects parahermitian matrx to be decomposed
default: 'Random'
Output parameter: none
0001 function SMDDemo(SourceSwitch); 0002 %SMDDemo(SourceSwitch); 0003 % 0004 % SMDDemo(SourceSwitch) demonstrates the polynomial EVD approximated by 0005 % the function SMD() for different type of parahermitian matrices. The 0006 % input parameter SourceSwitch can be selected as: 0007 % 0008 % SMDDemo('Random') produces an abritrary complex valued 5x5 polynomial 0009 % parahermitian matrix of order 10. 0010 % 0011 % SMDDemo('Simple1') decomposes a 3x3 real valued polynomial parahermitian 0012 % matrix R(z) of order 4: 0013 % | 1 -.4*z 0 | 0014 % R(z) = | -.4*z 1 .5*z^{-2} | 0015 % | 0 .5*z^2 3 | 0016 % 0017 % SMDDemo('Simple2') decomposes a 5x5 complex valued polynomial 0018 % parahermitian matrix R(z) of order 4: 0019 % | 1 .5*z^2 -.4j*z^{-1} 0 .2j*z^2 | 0020 % | .5*z^{-2} 1 0 0 0 | 0021 % R(z) = | .4j*z 0 3 .3z^{-1} 0 | 0022 % | 0 0 .3*z .5 0 | 0023 % |-.2j*z^{-2} 0 0 0 .25 | 0024 % This matrix admits a PEVD with a non-polynomial Gamma(z). 0025 % 0026 % The function generates a number of outputs to highlight diagonalisation 0027 % and spectral majorisation. 0028 % 0029 % Input parameter: 0030 % SourceSwitch selects parahermitian matrx to be decomposed 0031 % default: 'Random' 0032 % Output parameter: none 0033 0034 % S. Weiss, University of Southampton, 3/10/20 0035 0036 if nargin==0, 0037 SourceSwitch='Random'; 0038 end; 0039 0040 disp('Sequential Matrix Diagonalisation (SMD) Demo'); 0041 disp('--------------------------------------------'); 0042 0043 %--------------------------------------- 0044 % Define Scenario 0045 %--------------------------------------- 0046 % create parahermitian matrix to be decomposed 0047 if strcmp(SourceSwitch,'Random')==1, 0048 A = randn(5,5,6) + sqrt(-1)*randn(5,5,6);; 0049 R = PolyMatConv(A,ParaHerm(A)); 0050 elseif strcmp(SourceSwitch,'Simple1')==1, 0051 R = zeros(3,3,5); 0052 R(:,:,3) = diag([1 1 3]); 0053 R(1,2,1) = .5; 0054 R(2,1,5) = .5; 0055 R(1,2,2) = -.4; 0056 R(2,1,4) = -.4; 0057 elseif strcmp(SourceSwitch,'Simple2')==1, 0058 R = zeros(5,5,5); 0059 R(:,:,3) = diag([1 1 3 .5 .25]); 0060 R(1,2,1) = .5; 0061 R(2,1,5) = .5; 0062 R(1,3,4) = -.4*sqrt(-1); 0063 R(3,1,2) = .4*sqrt(-1); 0064 R(1,5,1) = .2*sqrt(-1); 0065 R(5,1,5) = -.2*sqrt(-1); 0066 R(4,3,2) = .3; 0067 R(3,4,4) = .3; 0068 else 0069 error('option for input parameter SourceSwitch not implemented'); 0070 end; 0071 0072 %--------------------------------------- 0073 % Calculate PEVD via SBR2 0074 %--------------------------------------- 0075 % parameters 0076 maxiter = 100; % max number of iterations 0077 epsilon = 0.001; % alternative stopping criterion 0078 mu = 0; % truncate true zeroes 0079 % decomposition 0080 [H,Gamma] = SMD(R,maxiter,epsilon,mu,'SMD'); 0081 0082 %--------------------------------------- 0083 % Display resulting matrices 0084 %--------------------------------------- 0085 % display original matrix 0086 figure(1); clf; 0087 L=size(R,3); 0088 t=(-(L-1)/2:(L-1)/2); 0089 PolyMatDisplay(abs(R),t); 0090 disp('Figure 1: original parahermitian matrix R(z)'); 0091 0092 % display paraunitary matrix 0093 figure(2); clf; 0094 L=size(H,3); 0095 PolyMatDisplay(abs(H)); 0096 disp('Figure 2: paraunitary matrix H(z)'); 0097 0098 % display diagonalised matrix 0099 figure(3); clf; 0100 L=size(Gamma,3); 0101 t=(-(L-1)/2:(L-1)/2); 0102 PolyMatDisplay(abs(Gamma),t); 0103 disp('Figure 3: diagonalised matrix Gamma(z)'); 0104 0105 % demonstrate approximate spectral majorisation 0106 figure(4); clf; 0107 Ndft = max([256,size(Gamma,3)]); 0108 P = PolyMatDiagSpec(Gamma,Ndft); 0109 plot((0:Ndft-1)/Ndft,10*log10(abs(P))); 0110 xlabel('normalised angular frequency \Omega/(2\pi)'); 0111 ylabel('power spectral densities / [dB]'); 0112 disp('Figure 4: PSDs of diagonal elements of Gamma(z)'); 0113 0114 %--------------------------------------- 0115 % Check result quantitatively 0116 %--------------------------------------- 0117 % remaining off-diagonal energy 0118 N1 = PolyMatNorm(R); 0119 N2 = PolyMatNorm(Gamma,'OffDiag'); 0120 disp(sprintf('ratio off-diagonal/total energy: %f (%f dB)',... 0121 [N2/N1, 10*log10(N2/N1)])); 0122 0123 % accuracy of decomposition and reconstruction 0124 R2 = PolyMatConv(ParaHerm(H),PolyMatConv(Gamma,H)); 0125 L = size(R,3); 0126 L2 = size(R2,3); 0127 Indices = (L2+1)/2 + ( -(L-1)/2:(L-1)/2 ); 0128 N3 = PolyMatNorm(R-R2(:,:,Indices)); 0129 disp(sprintf('error |R(z) - H~(z)Gamma(z)H(z)| = %f (%f dB)',... 0130 [N3, 10*log10(N3)]));