function [ ] = RFB_DNA(fil_nam)
% Detecting DNA microsatellites using the Ramanujan Filter Banks
% RFB_DNA(fil_nam) produces a period vs location plot for the DNA sequnce
% contained in the file fil_nam. Here, fil_nam must be the name of the 
% text file containing the DNA sequence in Fasta format.
% For example, fil_nam could be 'AC010136.txt'.
%
% Copyright (c) 2015, California Institute of Technology. [Based on
% research sponsored by ONR.] All rights reserved.
% Based on work from the lab of P P Vaidyanathan.
% Code developed by Srikanth V. Tenneti.

% No-endorsement clause: Neither the name of the California Institute 
% of Technology (Caltech) nor the names of its contributors may be used to 
% endorse or promote products derived from this software without specific 
% prior written permission.

% Disclaimer: THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 
% CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, 
% BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
% FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 
% HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 
% SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 
% TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
% EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.



%% Selecting The Filter Bank Parameters

Pmax = 50; % The largest period that needes to be detected
Rcq = 7; % Number of repeats in each Ramanujan filter
Rav = 3; % Number of repeats in each averaging filter

Penalty = []; % Peanlty function. Can be (optionally) used to set preference 
% to a certain set of periods in the time vs period plane 
for i=1:Pmax
    k = 1:i;
    k=k(gcd(k,i)==1);
    s = length(k);
    Penalty = cat(1,Penalty,(s)^0);
end

%% Constructing the Filter Bank

FR = cell(Pmax,1); %The set of Ramanujan Filters
FA = cell(Pmax,1); %The set of Averaging Filters

for i = 1:Pmax
    
    cq = zeros(i,1); % cq shall be the ith ramanujan sum sequnece.
k_orig = 1:i; k = k_orig(gcd(k_orig,i)==1);

for n = 0:(i-1)
    for a = k
    cq(n+1) = cq(n+1) + exp(1j*2*pi*a*(n)/i);
    end
end

cq=real(cq);

    FR{i,1} = repmat(cq,Rcq,1);
    FR{i,1} = FR{i,1}./norm(FR{i,1});
    FA{i,1} = repmat(ones(i,1),Rav,1);
    FA{i,1} = FA{i,1}./length(FA{i,1});
end


%% Reading the input signal

[xheader, xdata] = fastaread(fil_nam);
x = zeros(length(xdata),1);
y_final  = zeros(length(xdata),Pmax);
N_iter = 24;

all_assign = (perms([1 2 3 4]))';

for i_iter = 1:N_iter
fspec = 'Iteration %d out of %d.';
f_str = sprintf(fspec,i_iter,24);
disp(f_str);

nA = all_assign(1,i_iter);
nC = all_assign(2,i_iter);
nT = all_assign(3,i_iter);
nG = all_assign(4,i_iter);

for i = 1:length(xdata)
    switch xdata(i)
        case 'A' 
            x(i) = nA;
        case 'C' 
            x(i) = nC;
        case 'T' 
            x(i) = nT;
        case 'G' 
            x(i) = nG;
    end
end
            
            


% Computing the Outputs of the Filter Bank

y = zeros(length(x),Pmax);

for i = 1:Pmax
    
    y_temp = conv(x,FR{i,1},'same'); 
    y_temp = (abs(y_temp)).^2;
    y_temp = y_temp./Penalty(i);
    y_temp = conv(y_temp,FA{i,1},'same');
    y(:,i) = y_temp;
end

% Inner Thresholding 


y(:,1:2) = 0; % Periods 1 and 2 give strong features on the time vs period planes. Hence, zeroing them out to view the other periods better. 
y = y./max(y(:));
y(y<0.2)=0;

y_final = y_final + y;

end

% Outer Thresholding

y_final = y_final./max(y_final(:));
y_final(y_final<0.4)=0;

figure;imagesc(y_final');
colormap(flipud(gray));
disp('Please use the zoom tool in the figure window to view the detected micro-satellites'); 
end










