%
% Demonstrate discrete irregular reconstruction using irregrecon.m
%
% each run uses a different random sampling depending on case
%
% written by D.G. Long at BYU 30 Dec 2006
% updated by D.G. Long at BYU 30 Dec 2014
%
clear

% set parameters to generate synthetic signal 
M1=10;            % low pass coefficients
d1=3;             % set perfect value of d
N1=floor((2*M1+1)*d1);   % signal period (for exact uniform spacing) (N>2*M+1)
R1=(2*M1+1);      % number of irregular samples to generate (R>=2*M+1)

% override values for experimenting
%N1=N1+1;   % override period for non-integer regular sampling
%R1=floor(1.5*R1);  % number of irregular samples to generate (R>=2*M+1)
%R1=R1+1;   % number of irregular samples to generate (R>=2*M+1)


% set sampling location option flag
rflag=0; % uniformly distributed random distribution
%rflag=1; % regularly spaced locations
%rflag=2; % constrained random sampling


%
% generate synthetic 1d signal
% Fy is the "true" signal we want to reconstruct 
%

% generate ideal lowpass filter
lpf=zeros([1 N1]); lpf(1:M1+1)=1; lpf(N1-M1+1:N1)=1; 

% generate bandlimited random signal
if 0
  y=randn([1 N1]);          % random signal, not bandlimited
else
  tmp=(rand([1 N1])-0.5).*lpf; % random phase over lpf bandwidth
  tmp(1)=0; tmp(2:M1+1)=-tmp(N1:-1:N1-M1+1); % make phase odd
  y=real(ifftshift(ifft(exp(j*pi*tmp))))*N1/M1; % bandlimited signal
end
y=y+0.25*sin(M1*pi*([1:N1]-1)/N1); % optionally add single frequency signal
Fy=real(ifft(fft(y).*lpf));% apply lp filter to generate bandlimited signal

% generate R1>2*M1+1 random samples, ensuring that they are unique
% and within one period (1:N1).  Without loss of generality,
% the first sample is set to be at 1.
sloc=randsample1d(R1,N1,M1,rflag);

sl=[1:floor(d1):floor(d1)*(2*M1+1)];  % index to uniform samples
sl1=[sl sl(1)+N1];       % array used for computing max location spacing
sloc1=[sloc sloc(1)+N1]; % array used for computing max location spacing

% compute the actual irregular samples
x=[0:N1-1];
fr=Fy(sloc)';   % column vector of the irregularly sampled signal
xr=x(sloc)';    % column vector of the sample locations

if 0 % optionally print sample locations
  mdelta=max(diff(sloc1));
  fprintf('Sample locations  max_delta=%d\n',mdelta);
  fprintf('%4d %4d %4d\n',[(1:R1); sl; sloc]);
end

if 1 % optionally plot sample location summary
  myfigure(3);clf
  subplot(2,1,1)
  plot([1 N1],[0 0],'k');
  hold on;
  plot(sl,sl*0,'b.');
  plot(sloc,sloc*0,'rx')
  hold off;
  xlabel('index');
  title('Sample locations b=uniform r=irregular')
  subplot(2,1,2)
  plot([1 R1],[d1 d1],'k');
  hold on;
  plot(diff(sl1),'b');
  plot(diff(sloc1),'r');
  hold off
  mdelta=max(diff(sloc1));
  title(sprintf('Sample location spacing, max delta=%d of %d',mdelta,d1))
  %print -dpng fig6.png
end

%%
%% reconstruct signal using irregrecon, return optional parameters for plotting
%%
[Frecon fu xu Do D K]=irregrecon(fr,xr,N1,M1);

% reconstruction parameters used in irregrecon
N=length(Frecon);
M=(length(fu)-1)/2;
R=length(fr);
d=N/(2*M+1);    % uniform sample spacing (in samples)

% plot the signal and its spectrum
myfigure(1);clf

% signal plot
subplot(2,1,1)
% first plot original signal
plot([x(1),x(end)],[0,0],'b');
hold on
for k=1:length(Fy)
  plot([x(k) x(k)],[0 Fy(k)],'b');
end
% plot equivalent uniform samples		    
plot(xu,0*fu,'k.');   % plot equiv uniform sample locations on axis
for k=1:length(fu)
  plot([xu(k) xu(k)],[0 fu(k)],'c'); % equiv uniform sample values
end
plot(xu,fu,'c.');     % plot equiv uniform sample values
% plot irregular sample values and locations
plot(xr,fr,'r.');
% plot reconstructed signal
plot(x,Frecon,'g'); 
hold off
axis([0 N-1 -1.5 1.5])
xlabel('sample index 0:N-1 over one period')
ylabel('value')
title(sprintf('Period N=%d  Bandlimit M=%d  Samples R=%d  d=%d',N,M,R,d));


% spectrum plot
subplot(2,1,2)
F=fft(Fy);    xfn=(0:N-1)-floor(N/2);
Ffu=fft(fu);  xfn2=(0:2*M)-M;
% original signal spectrum
plot(xfn,fftshift(20*log10(abs(F)+1.e-25)));
hold on
plot(xfn,fftshift(20*log10(abs(F)+1.e-25)),'b.');
% spectrum of equivalent uniform samples
plot((0:2*M)-M,fftshift(20*log10(abs(d*Ffu)+1.e-25)),'k.');
hold off
axis([-floor(N/2) floor(N/2) -20 30])
xlabel('frequency index (-N/2:N/2) -> (-\pi:\pi)')
title('Signal magnitude spectrum');
ylabel('dB');

%print -dpng fig1.png

% print spectral coefficients
fi=[-M:M]; ind=[N-M+1:N 1:M+1];
ind2=[M+2:2*M+1 1:M+1];
Ffu=fft(fu);

if 0 % optionally print spectral coefficients
  fprintf('Signal spectrum:\n');
  fprintf('  j           a_j(N)                  d a_j(2*M+1)     Percent err\n');
  fprintf('%4d %9.3f + %9.3f i   %9.3f + %9.3f i  (%9.3f)\n',...
      [fi' real(F(ind))' imag(F(ind))' real(d*Ffu(ind2)) imag(d*Ffu(ind2)) ...
	100*abs(F(ind).'-d*Ffu(ind2))./abs(F(ind))']');
end

if 0 % optionally print equivalent uniform samples
  % note that equivalent uniform samples are only "samples" of
  % Fy if d is an integer
  fi=1:2*M+1;
  fprintf('Equivalent signal samples: d=%d  N=%d  M=%d  R=%d  minR=%d\n',d,N,M,R,2*M+1);
  fprintf('   j    x(j d)    fu(j) [x(j d)] Fy([x])  Fr([x])   err*100\n');
  fprintf('%4d %9.2f %9.2f %4d %9.2f %9.2f %9.4f\n',[fi' xu fu ...
	floor(xu) Fy(floor(xu+1)).' Frecon(floor(xu+1)) ...
	(Fy(floor(xu+1))'-Frecon(floor(xu+1)))*100 ]'); 
end

% compute and print summary errors
ndelta=min(diff(sloc));
mdelta=max(diff(sloc));
RMS_recon=100*sqrt(sum((Fy'-Frecon).^2))/sqrt(sum(Fy.^2));
fprintf('Equivalent signal samples: d=%d  N=%d  M=%d  R=%d  minR=%d\n',d,N,M,R,2*M+1);
fprintf('Sample location spacing d=%d  min_spacing=%d  max_spacing=%d\n',d,ndelta,mdelta);
fprintf('RSS_reconstruction error in percent: %f\n',RMS_recon);
fprintf('Reconstruction Do matrix condition number: %f\n',K);


if 0 % optinally use SVD to compute Do (pseudo) inverse and eigenvalues
     % and plot matrix eigenvalues
  [U S V]=svd(Do);
  thres=1.e-10;
  Sinv=S';
  ind=find(Sinv>thres);
  Sinv(ind)=1./Sinv(ind);
  Doinv=V*Sinv*U';
  
  myfigure(4);cf
  k=min(size(S));
  plot(log10(S(sub2ind(size(S),[1:k],[1:k]))));
  xlabel('order')
  ylabel('Log10 Eigenvalue of Do')
  title('Eigenvalues of irregular reconstruction matrix');
end

if 0 % optionally plot reconstruction response function for given samples
  myfigure(6);clf
  plot(Doinv');
  title('Response for each irregular measurement');
  xlabel('Equivalent uniform sample index 1:2*M+1');
  ylabel('Response');

  myfigure(7);clf
  surf(Doinv);
  ylabel('Irregular sample index 1:R');
  xlabel('Equivalent uniform sample index 1:2*M+1');
  title('Response for each irregular measurement');
end
