%
% Demonstrate general discrete irregular 1d reconstruction using recon.m
% Includes both aperture and no-aperture cases
%
% different random sampling depending on sampling 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=30;            % 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+3;   % 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

% set aperture function parameters
ApCase=3;         % aperture function test case (-1..9)
L=3;              % aperture length (2*L+1) setting

% define aperture functions for each sample
Ap=zeros([R1 2*L+1]);  % aperture matrix
switch ApCase
case -1 % no aperture
  Ap(:,L+1)=1;
  Na=0;
case 0  % delta function
  Ap(:,L+1)=1;
  Na=1;
case 1  % simple 3 term window
  Ap(:,L)=0.25;
  Ap(:,L+1)=1; % center
  Ap(:,L+2)=0.25;
  Na=1;
case 2 % Fejer (triangle) window
  k=1-abs([-L:L])'/(L+1);
  for n=1:R1
    Ap(n,:)=k'; 
  end
  Na=1;
case 3 % Hamming window
  k=pi*[-L:L]'/L;
  for n=1:R1
    Ap(n,:)=0.56+0.48*cos(k'); 
  end
  Na=1;
case 4 % multiple Fejer window sizes
  k=1-abs([-L:L])'/(L+1);
  for n=1:2:R1
    Ap(n,:)=k'; 
  end
  k=[0 1-abs([1-L:L-1])/L 0]';
  for n=2:2:R1
    Ap(n,:)=k';
  end
  Na=R1;
case 5 % multiple Hamming window sizes
  k=pi*[-L:L]'/L;
  for n=1:2:R1
    Ap(n,:)=0.56+0.48*cos(k'); 
  end
  k=pi*[1-L:L-1]'/(L-1);
  for n=2:2:R1
    Ap(n,1)=0;Ap(n,2*L+1)=0;
    Ap(n,:)=[0 0.56+0.48*cos(k') 0]; 
  end
  Na=R1;
case 6 % mix of cases
  Ap(:,L+1)=1;
  for n=2:6:R1
    Ap(:,L)=0.25;
    Ap(:,L+2)=0.25;
  end
  k=1-abs([-L:L])'/(L+1);
  for n=3:6:R1
    Ap(n,:)=k'; 
  end
  k=[0 1-abs([1-L:L-1])/L 0]';
  for n=4:6:R1
    Ap(n,:)=k';
  end
  k=pi*[-L:L]'/L;
  for n=5:6:R1
    Ap(n,:)=0.56+0.48*cos(k'); 
  end
  k=pi*[1-L:L-1]'/(L-1);
  for n=6:6:R1
    Ap(n,1)=0;Ap(n,2*L+1)=0;
    Ap(n,:)=[0 0.56+0.48*cos(k') 0]; 
  end
  Na=R1;
case 7 % multiple Hamming windows
  k=pi*[-L:L]'/L;
  for n=1:R1
    Ap(n,:)=0.56+.48*cos(k'); 
  end
  k=pi*[1-L:L-1]'/(L-1);
  for n=floor(R1/2)+1:R1
    Ap(n,:)=[0 0.56+0.48*cos(k') 0]; 
  end
  Na=R1;
case 8  % simple 3 term averaging window
  Ap(:,L)=1;
  Ap(:,L+1)=1; % center
  Ap(:,L+2)=1;
  Na=1;
case 9  % dual averaging windows
  Ap(:,L)=1;
  Ap(:,L+1)=1; % center
  Ap(:,L+2)=1;
  for n=1:floor(R1/2)
  %for n=floor(R1/2)+1:R1
    Ap(n,L-1)=1;
    Ap(n,L+3)=1;
  end
  Na=R1;
otherwise % delta function
  Ap(:,L+1)=1;
  Na=1;
end
  
% normalize apertures
for n=1:R1
  Ap(n,:)=Ap(n,:)/sum(Ap(n,:));
end

%
% generate a synthetic signal
% Fy is the "true" signal we want to reconstruct 
%
lpf=zeros([1 N1]); lpf(1:M1+1)=1; lpf(N1-M1+1:N1)=1; % generate ideal lowpass filter

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

% generate R1>2*M1+1 random (depending on rflag) samples, 
% ensuring that they are disjoint (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];
xr=x(sloc)';    % column vector of the sample locations
%xr=sloc(:)-1;    % column vector of the sample locations
fr=Fy(sloc)';   % column vector non-aperture-filtered samples

fr1=fr; % save for testing

if 0 % optionally print sample locations
  mdelta=max(diff(sloc1));
  fprintf('Sample locations  max_spacing=%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))
end
  
% plot apertures
myfigure(2);clf
if Na>0
  k=[-L:L];
  subplot(2,1,1)
  for n=1:R1
    temp=Ap(n,:);
    plot(k,temp,'r');
    if n==1, hold on; end
  end
  hold off
  xlabel('realtive index');
  ylabel('value');
  title(sprintf('Apertures (ApCase=%d)',ApCase));

  subplot(2,1,2)
  for n=1:R1
    temp=zeros([1 N1]);
    ind0=sloc(n)+k;
    ind1=find(ind0>N1);
    ind0(ind1)=ind0(ind1)-N1;
    ind1=find(ind0<1);
    ind0(ind1)=ind0(ind1)+N1;
    temp(ind0)=Ap(n,:); % zero pad apertures
    plot(1:N1,temp,'r');
    if n==1, hold on; end
    plot(sloc(n),temp(sloc(n)),'b*');
  end
  hold off
  xlabel('index');
  ylabel('value');
  title('individual apertures with sample locations')
else
  clf;
end

% generate aperture-filtered samples
k=[-L:L];
for n=1:R1
  ind0=sloc(n)+k;
  ind1=find(ind0>N1);
  ind0(ind1)=ind0(ind1)-N1;
  ind1=find(ind0<1);
  ind0(ind1)=ind0(ind1)+N1;
  fr(n)=sum(Ap(n,:).*Fy(ind0)); % aperture filtered sample
end

% generate aperture pointer array
if Na>1
  xa=1:length(xr);
else
  xa=zeros(size(1:length(xr)))+Na;
end

d1=N1/(2*M1+1);    % uniform sample spacing (in samples)
fprintf('Equivalent signal samples: d=%d  N=%d  M=%d  R=%d  minR=%d\n',d1,N1,M1,R1,2*M1+1);
ndelta=min(diff(sloc));
mdelta=max(diff(sloc));
fprintf('Sample location spacing d=%d  min_spacing=%d  max_spacing=%d\n\n',d1,ndelta,mdelta);


%%
%% reconstruct signal using recon, return optional parameters for plotting
%%

% repeat for each method to compute timing and compare accuracy
for method=0:3

  tstart=tic;
  
  Frecon=recon(fr,xr,N1,M1,method,Na,xa,Ap);

  dtime=toc(tstart);  % compute timing of call

  RMS_recon1=sqrt(sum((Fy.'-Frecon).^2)/length(Fy));
  RMS_recon2=100*RMS_recon1/sqrt(sum(Fy.^2)/length(Fy));
  
  fprintf('Method: %d  time: %f\n',method,dtime);
  fprintf('RSS_reconstruction error: %f %f%%\n\n',RMS_recon1,RMS_recon2);

  if method==0 % generate signal plots for default method

    % reconstruction parameters
    N=N1;
    M=M1;
    R=length(fr);
    d=N/(2*M+1);    % uniform sample spacing (in samples)
    xu=d*[0:2*M]'+1;% locations of estimated uniform samples
    fu=Fy(xu);

    % 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.');
    % spectra of apertures
    for k=1:R
      temp=xfn*0;
      temp(1:2*L+1)=Ap(k,:);  % zero pad apertures
      plot(xfn,fftshift(20*log10(abs(fft(temp))+1.e-25)),'r');
    end
    hold off
    axis([-floor(N/2) floor(N/2) -20 30])
    xlabel('frequency index (-N/2:N/2) -> (-\pi:\pi)')
    title('Signal and apertures magnitude spectra');
    ylabel('dB');


    % compute 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
      disp('Signal spectrum:');
      disp('  j           a_j(N)                  d a_j(2*M+1)     Percent err');
      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 min=%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

  end

end

