# SMDDemo

## PURPOSE

SMDDemo(SourceSwitch);

## SYNOPSIS

function SMDDemo(SourceSwitch);

## DESCRIPTION

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

## CROSS-REFERENCE INFORMATION

This function calls:
This function is called by:

## SOURCE CODE

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

Generated on Wed 03-Dec-2014 19:38:11 by m2html © 2005