Petre Stoica and Tomas Sundin, "Nonparametric NMR Spectroscopy", Journal of Magnetic Resonance", vol 152. pp 57-69, 2001.
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)))
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)))