%%% copyright 2012, Emilio Kropff
%%% output:  first dim = number of frequencies
%%%         second dim = number of trials
%%%          third dim = number of timestamps



function h=bandpass_and_hilbert(x,flow,fhigh,sampling)
    l=size(x,1); 

    NFFT= 2^nextpow2(l); % power of 2 is a faster fft 

    x_fft=fft(x,NFFT);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Hilbert Transform, taken from the MATLAB function hilbert(). h is
    % the 'filter' applied on the fft of the signal. The filter
    % consists in multiplying the fft by i*sign(freq). But the output
    % is hilbert(x)=x-i*transform(x), which results in a total filter that multyplies
    % by 2 positive frequencies and by 0 negative frequencies (and by 1 the
    % two central frequencies).
    % http://en.wikipedia.org/wiki/Hilbert_transform#Relationship_with_
    % the_Fourier_transform
    hilbert_filter  = zeros(NFFT,1); 
    % NFFT is even (it is a power of 2) so there are two central frequencies
    hilbert_filter([1 NFFT/2+1]) = 1; 
    hilbert_filter(2:NFFT/2) = 2;   %positive freqs (upper part of the circle)        

    
    fft_freqs=sampling/2*(linspace(0,1,NFFT/2+1))'; %This includes only the upper part of the circle
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Bandpass Filter
    bandpass_filter=bsxfun(@gt,fft_freqs,flow) & bsxfun(@lt,fft_freqs,fhigh);
    bandpass_filter=[bandpass_filter; bandpass_filter(end-1:-1:2,:)]; 

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Finally, perform the bandpass and hilbert transform in one step
    % (blocked)
    
%     nfreqmax=floor(.5*10^8/numel(x_fft));
    nfreqmax=floor(.1*10^8/numel(x_fft));
    fstart=1; fend=min(numel(flow),fstart-1+nfreqmax);
    h=NaN(l,numel(flow));
    while fstart<=fend
        haux=ifft(bsxfun(@times,x_fft,bsxfun(@times,bandpass_filter(:,fstart:fend),hilbert_filter)));%combined_filter(:,ones(1,car.nlaps)));
        h(:,fstart:fend)=haux(1:l,:);
        fstart=fend+1;
        fend=min(numel(flow),fstart-1+nfreqmax);
    end

%     h=squeeze(h);
%     if nargout>1
%         bpx=zeros(nf,nx,l);
%         for f=1:nf
%             bandpass_filter=(fft_freqs>flow).*(fft_freqs<fhigh);
%             bandpass_filter=[bandpass_filter; bandpass_filter(end-1:-1:2)];             
%             hilbert_transf=ifft(bsxfun(@times,x_fft,bandpass_filter));
%             bpx(f,:,:)=hilbert_transf(1:l,:);
%         end
%         bpx=squeeze(bpx);
%     end