2D Capon and APES MATLAB examples from JMR:152 57-69 (2001)

Petre Stoica and Tomas Sundin, "Nonparametric NMR Spectroscopy", Journal of Magnetic Resonance", vol 152. pp 57-69, 2001.


2D Capon

% capon_2dtest.m
% Fred J. Frigo
% Marquette University,   Milwaukee, WI  USA
% Copyright 2002 - All rights reserved.
%
% May 25, 2002 - implements 2D Capon example by Stoica and Sundin
%                Journal of Magnetic Resonance (JMR)  #152, Fig 2  p62
%
% See notes from Dr. James A. Heinen  4/25/2002 - 5/16/2002
%
%

function capon_2dtest()

clear all;
close all;
 

% Get start time
ctime = clock;
start_time = sprintf('%d:%2.2d:%2.2d',ctime(4), ctime(5),round(ctime(6)))

i = sqrt(-1.0);

% Parameters given on page JMR: #152  p61
N= 768;
%w1=0.1885;   % original value from JMR: 152 p61 - did not work for N=768, M=256
w1=0.184077;  % Alternate value for N=768, M=256
alpha1=0.016;
a1=1.0;
%w2=0.2136;   % original value from JMR: 152 p61 - did not work for N=768, M=256
w2=0.220892;  % Alternate value for N=768 and M=256
alpha2=0.026;
a2=2.5;
noise_weight= 0.001;
x=[1:N];

% Capon filter length
M = 256;

% Output Length
L=N-M;

% Ouput plot parameters
omega_start = 1;
omega_stop = 40;

% create input signal and plot
for n = 1:N
    noise(n) = rand(1)*noise_weight + i*rand(1)*noise_weight;
    sig(n)= a1*exp( (-alpha1 + (i*w1))*(n-1) ) + ...
            a2*exp( (-alpha2 + (i*w2))*(n-1) ) + noise(n);
end
 
figure(1);
plot(x,real(sig), 'b', x, imag(sig), 'r--');
title('Simulated Input signal');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform 2D Capon on the simulated signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Form the sample covariance matrix; JMR:152, p60 [22]
% See also Intro to Spectral Analysis, Stoica, p199, 1997 (5.4.18)
R=zeros(M,M);
for index = 1:(L+1)
   R = R + (sig(index:index+M-1).'*conj(sig(index:index+M-1)));
end
%R=R/L;  % The divide by L is not in JMR: 152

% compute the inverse of R
IR=inv(R);

% Ramp for exponentials used to determine Y(alpha, omega) Stoica JMR:152, p59 [8]
ramp=linspace(1,L,L);
 
% Set up y(t) defined in JMR:152, p 58 [4]
for m=1:M
   yt(m,:)=sig(m:m+L-1);  % Stoica calls this y(t)
end

% Compute Y(alpha,omega) for different values of alpha
P=35;
alpha_init = 0.0;
alpha_step = 0.001;
for p=1:P
   % compute alpha for this loop
   alpha = alpha_init + (p-1)*alpha_step;
 
   % Compute scale factor: L(alpha) JMR:152 p59 [9]
   if alpha == 0.0
       % to prevent divide by zero, take limit to get L
       Lp = (L*1.0);
   else
       Lp=exp(-2.0*alpha)*((exp(-2.0*alpha*L) - 1.0)/(exp(-2.0*alpha)- 1.0));
   end

 
   % Compute Y(alpha, omega) given by Stoica JMR:152 p 59 [8]
   for L_index=omega_start:omega_stop
 
      % Compute s(alpha,omega) given by Stoica JMR:152 p 60, [24]
      L_omega = (2.0*pi)*L_index/(L*1.0);
      for m=1:M
          s(m)=exp( (-alpha + (i*L_omega))*(m-1) );
      end
 
      % Compute hcapon given by Stoica JMR:152 p 60, [23]
      hcapon=(IR*s.')/(conj(s)*IR*s.');
 
      % compute exponential for [8]
      myexp = exp( (-alpha +(-i*L_omega))*ramp);

      % Compute Y(alpha,omega) from JMR:152  p 59 [8]
      for m=1:M
          ytexp=yt(m,:).*myexp;
          yalpha_omega(m) = sum(ytexp)/Lp;
      end
 
      capon_results(L_index) = hcapon'*yalpha_omega.';
   end
 
   %Capon Energy Spectrum - Stoica  JMR:152 p59 [10]
   ces(p,:) = sqrt(abs(capon_results)*Lp);
 
end

% Surface and Contour plot
figure(2);
contour(ces);
figure(3);
surf(ces);

% Get stop time
ctime = clock;
stop_time = sprintf('%d:%2.2d:%2.2d',ctime(4), ctime(5),round(ctime(6)))


2D APES

% apes_2dtest.m
% Fred J. Frigo
% Marquette University,   Milwaukee, WI  USA
% Copyright 2002 - All rights reserved.
%
% May 25, 2002 - implements 2D APES example by Stoica and Sundin
%                Journal of Magnetic Resonance (JMR)  #152, Fig 2  p62
%
% See notes from Dr. James A. Heinen  4/25/2002 - 5/16/2002
%

function apes_2dtest()

clear all;
close all;

% Get start time
ctime = clock;
start_time = sprintf('%d:%2.2d:%2.2d',ctime(4), ctime(5),round(ctime(6)))

i = sqrt(-1.0);

% Parameters given on page JMR: #152  p61
N= 768;
%w1=0.1885;   % original value from JMR: 152 p61 - did not work for N=768, M=256
w1=0.184077;  % Alternate value for N=768, M=256
alpha1=0.016;
a1=1.0;
%w2=0.2136;   % original value from JMR: 152 p61 - did not work for N=768, M=256
w2=0.220892;  % Alternate value for N=768 and M=256
alpha2=0.026;
a2=2.5;
noise_weight= 0.001;
x=[1:N];

% APES filter length
M = 256;

% Output Length
L=N-M;

% Ouput plot parameters
omega_start = 1;
omega_stop = 40;

% create input signal and plot
for n = 1:N
    noise(n) = rand(1)*noise_weight + i*rand(1)*noise_weight;
    sig(n)= a1*exp( (-alpha1 + (i*w1))*(n-1) ) + ...
            a2*exp( (-alpha2 + (i*w2))*(n-1) ) + noise(n);
end
 
figure(1);
plot(x,real(sig), 'b', x, imag(sig), 'r--');
title('Simulated Input signal');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform 2D APES on the simulated signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Form the sample covariance matrix; JMR:152, p60 [22]
% See also Intro to Spectral Analysis, Stoica, p199, 1997 (5.4.18)
R=zeros(M,M);
for index = 1:(L+1)
   R = R + (sig(index:index+M-1).'*conj(sig(index:index+M-1)));
end
%R=R/L;  % The divide by L is not in JMR: 152

% compute the inverse of R
IR=inv(R);

% Ramp for exponentials used to determine Y(alpha, omega) Stoica JMR:152, p59 [8]
ramp=linspace(1,L,L);
 
% Set up y(t) defined in JMR:152, p 58 [4]
for m=1:M
   yt(m,:)=sig(m:m+L-1);  % Stoica calls this y(t)
end

% Compute Y(alpha,omega) for different values of alpha
P=35;
alpha_init = 0.0;
alpha_step = 0.001;
for p=1:P
   % compute alpha for this loop
   alpha = alpha_init + (p-1)*alpha_step;
 
   % Compute scale factor: L(alpha) JMR:152 p59 [9]
   if alpha == 0.0
       % to prevent divide by zero, take limit to get L
       Lp = (L*1.0);
   else
       Lp=exp(-2.0*alpha)*((exp(-2.0*alpha*L) - 1.0)/(exp(-2.0*alpha)- 1.0));
   end

 
   % Compute Y(alpha, omega) given by Stoica JMR:152 p 59 [8]
   for L_index=omega_start:omega_stop
 
      % Compute s(alpha,omega) given by Stoica JMR:152 p 60, [24]
      L_omega = ((2.0*pi)*L_index)/(L*1.0);
      for m=1:M
          s(m)=exp( (-alpha + (i*L_omega))*(m-1) );
      end
 
      % compute exponential for [8]
      myexp = exp( (-alpha +(-i*L_omega))*ramp);

      % Compute Y(alpha,omega) from JMR:152  p 59 [8]
      for m=1:M
          ytexp=yt(m,:).*myexp;
          yalpha_omega(m) = sum(ytexp)/Lp;
      end
 
      % Compute Q as given by Stoica JMR:152  p60  [27]
      % Q = R + (yalpha_omega.'*conj(yalpha_omega));
      % IQ = inv(Q);
 
      % Compute inv(Q) using matrix inversion lemma  JMR:152,  p61 [30]
      IQ = IR + ((IR*Lp*yalpha_omega.'*conj(yalpha_omega)*IR)/...
                 (Lp*conj(yalpha_omega)*IR*yalpha_omega.'));
 
      % Compute hapes given by Stoica JMR:152 p 61, [29]
      hapes = (IQ*s.')/(conj(s)*IQ*s.');
 
      % compute APES results by Stoica JMR:152, p59  [7]
      apes_results(L_index) = hapes'*yalpha_omega.';
   end
 
   %APES Energy Spectrum - Stoica  JMR:152 p59 [10]
   apes_es(p,:) = sqrt(abs(apes_results)*Lp);
 
end

% Surface and Contour plot
figure(2);
contour(apes_es);
figure(3);
surf(apes_es);

% Get stop time
ctime = clock;
stop_time = sprintf('%d:%2.2d:%2.2d',ctime(4), ctime(5),round(ctime(6)))



last updated:  August 4, 2002