function [u,varargout]=igft(data,B,varargin)
% Computes the reconstruction of input data using a matrix as a
% user-defined orthogonal basis. A warning is produced if the input 
% matrix B is non-orthogonal.  Plots results as an option.  The
% reconstructed function and amplitude spectra (optional) are returned.
% 
% [u,varargout]=igft(data,B,varargin)
% 
% INPUT
% data: an Mx1 vector
%    B: an MxN orthogonal matrix.  If B is non-orthongal, and output will
%       still be produced, but there will be a warning message.  Any value
%       for M may be used, but it is recommended that M > N is used, as
%       this function is intended to be used as a tool for approximating
%       data using 'best bases'.
% plot: (Optional) string 'plot': an option to plot the amplitude spectra.
%       Recommended to visualize how basis reconstructs given data.
% 
% OUTPUT
%  u:  the reconstructed data.
% amp:  (Optional) the generalized Fourier amplitude spectra.
% 
% Example:
%
% B=rand(100,50);
% Q=orth(B);
% data=rand(100,1);
% [u,amp]=igft(data,Q,'plot');
% 
% Also available as part of VSPLAB from site:
% http://www.ess.washington.edu/~joshuadc

opt='none';
Md=length(data);
[M,N]=size(B);
Ndisp=min([N,35]);
spctr=zeros(N,1);
colors={'k','r','b'};

if(Md~=M)
    error('data and basis vectors do not match');
end;

if(nargin==3)
    if(strcmp(varargin{1},'plot'))
        opt=varargin{1};
    end;
end;

ercount=0;
u=zeros(M,1);
for k=1:N
    if(k<N)
        if(B(:,k)'*B(:,k+1)>10^(-5))
            ercount=ercount+1;
            if(ercount==1)
            disp('Your basis is not orthonormal to tolerance of 10^(-5)')
            end;
        end;
    end;
    spctr(k)=B(:,k)'*data;
    u=u+(spctr(k))*B(:,k);
end;

if(strcmp('plot',opt))
    figure; hold on; subplot(2,2,1:2);
    plot(data,'k','linewidth',1.8); hold on;
    plot(u,'r','linewidth',1.8);
    set(gcf,'Color', [1,1,1]);
    set(gca,'linewidth',2);
    legend('Input Data','Reconstruction','Location','NorthWest')
    title('Reconstruction of Data','fontsize',12);
    ylabel('Reconstruction and Original','fontsize',12);
end;

t=0:N-1;
if(strcmp('plot',opt))
    subplot(2,2,3);
    set(gcf,'Color', [1,1,1]);
    for k=1:Ndisp
        s=stem(t(k),abs(spctr(k)),char(colors(mod(k-1,3)+1))); 
        set(s,'linewidth',2);
        hold on;
    end;
    set(gca,'Color', [0.85,0.9,0.95]);
    set(s,'linewidth',2);
    title(sprintf('First %i Fourier Coefficients from Basis',Ndisp),'fontsize',12);
    ylabel('Fourier Spectral Amplitude','fontsize',12);
    xlabel('Coefficient Index','fontsize',12);
    set(gca,'fontsize',12);
    subplot(2,2,4);
    plotXmatrix(B(:,1:Ndisp),'align');
    set(gca,'linewidth',2);
    set(gcf,'Color', [1,1,1]);
    title(sprintf('First %i Basis Functions used in Reconstruction',Ndisp),'fontsize',12);
    ylabel('Basis Index','fontsize',12);
    xlabel('Sample Number','fontsize',12);
end;
if(nargout==2)
    varargout{1}=abs(spctr);
end;